# Universality of striped morphologies

E. Edlund<EMAIL_ADDRESS>M. Nilsson Jacobi<EMAIL_ADDRESS>Complex Systems Group, Department of Energy and Environment, Chalmers
University of Technology, SE-41296 Göteborg, Sweden

###### Abstract

We present a method for predicting the low-temperature behavior of spherical
and Ising spin models with isotropic potentials. For the spherical model the
characteristic length scales of the ground states are exactly determined but
the morphology is shown to be degenerate with checkerboard patterns, stripes
and more complex morphologies having identical energy. For the Ising models we
show that the discretization breaks the degeneracy causing striped
morphologies to be energetically favored and therefore they arise universally
as ground states to potentials whose Hankel transforms have nontrivial minima.

The study of pattern formation in simple systems has received much attention
the last 20 years. Not only are the pictures visually arresting, producing
reviews with considerable artistic qualities Seul and Andelman (1995); Bowman
and Newell (1998), but spatially inhomogeneous phases present difficulties for
standard theoretical methods, calling for new principles to describe the
physics of systems exhibiting them Emery  _et al._ (1999). Whatever weight can
be assigned to the former as an explanation for the interest, the latter is
justification enough, especially as the experimental evidence of stripes,
spots and checkerboards in practically important materials abound. Examples
include lipid monolayers Keller and McConnell (1999), adsorbates on metals
Kern  _et al._ (1991), and various magnetic fluids Rosensweig  _et al._
(1983); Seul and Wolfe (1992). Striped phases are also hypothesized to play a
role in the high-temperature superconductivity of transition metal oxides
Tranquada  _et al._ (1995); Orenstein and Millis (2000).

Many experimental systems displaying heterogeneous patterns involve a
competition between long- and short-range interactions De’Bell  _et al._
(2000) and most theoretical work concentrate on specific examples of such
interactions Garel and Doniach (1982); Löw  _et al._ (1994); Grousson  _et
al._ (2000); Stoycheva and Singer (2000); Giuliani  _et al._ (2006), e.g. spin
models with Hamiltonians on the form

$H=K\sum_{j}s_{j}^{2}-L\sum_{\langle
i,j\rangle}s_{i}s_{j}+\frac{Q}{2}\sum_{i\neq
j}\frac{s_{i}s_{j}}{r_{ij}^{\alpha}},$ (1)

where the spins typically represent some coarse-grained feature of the system
of interest, for example local charge density in a Mott insulator Emery  _et
al._ (1999) or phases in a Langmuir film Weis and McConnell (1984). However,
the qualitative success of such models may have little to do with the
underlying physics as noted by Zaanen in the context of Mott insulators Zaanen
(1998). Indeed, the same general behavior can be observed in models with for
example only short-ranged, purely repulsive forces Malescio and Pellicane
(2003). An explanation for the universality of striped morphologies must
therefore be independent of specific details of the involved forces. The aim
of this Letter is to present such a general treatment. As expected our method
shows that stripes appear naturally for large classes of models, but the added
generality also leads to new tools allowing us to design potentials with
desired properties.

Here we study a generic Hamiltonian with isotropic pairwise interactions

$H=\sum_{ij}^{N}V_{ij}s_{i}s_{j}$ (2)

where $V_{ij}=V(\left|i-j\right|)$ is a matrix representation of the potential
that only depends on the distance between spins $i$ and $j$, here denoted
$\left|i-j\right|$, with (1) as a special case. Depending on considerations
regarding experimental fit or theoretical ease, one may take the spins in (2)
to assume continuous values with the restriction $\sum_{i}s_{i}^{2}=N$,
corresponding to a spherical model Baxter (1982), or take values from some
finite set, where $s\in\\{\pm 1\\}$ and $s\in\\{0,\pm 1\\}$ are the most
common choices, equivalent to different Ising models Löw  _et al._ (1994);
Giuliani  _et al._ (2006).

Consider now the spherical model. Due to the quadratic form of the Hamiltonian
(2) and the constraint, ${\bf s}^{T}V{\bf s}$ and ${\bf s}^{T}{\bf s}=N$ in
matrix notation, the ground state is given by an eigenvector of the
interaction matrix $V$ corresponding to the lowest (energy) eigenvalue Lay
(2003). The central observation for our analysis is the existence of a common
basis of eigenvectors for all radial potentials, namely the Fourier basis.

Figure 1:  Reflection symmetry of triangles on a lattice causes the
interaction matrices from any two radial potentials to commute.

To prove this we start by recalling the fact that if two matrices commute,
then it is possible to find a set of eigenvectors that simultaneously
diagonalize them Horn and Johnson (1990). Consider the commutator for two
interaction matrices $V$ and $W$:

$\sum_{k^{\prime}}V_{ik^{\prime}}W_{k^{\prime}j}-\sum_{k}W_{ik}V_{kj}.$ (3)

Each term, $V_{ik^{\prime}}W_{k^{\prime}j}$, in the first sum can be
represented by a triangle, A in Fig. 1. Assuming that the lattice is periodic
or infinite, there will for each such triangle exist a unique triangle B,
constructed as a reflection of A in the line equidistant from point $i$ and
$j$ (dashed in the figure), corresponding to the term $W_{ik}V_{kj}$ in the
second sum. From the reflection symmetry and the pure radial dependence of the
potential it follows that $V_{ik^{\prime}}W_{k^{\prime}j}-W_{ik}V_{kj}=0$,
which proves that $V$ and $W$ commute.

It now suffices to find a set of eigenvectors for a particular potential.
Perhaps the simplest choice is a nearest neighbor interaction, $V_{ij}=1$ if
$\left|i-j\right|=1$ and $V_{ij}=0$ if $\left|i-j\right|>1$. If we in addition
make an appropriate choice of self-interaction, which only shifts the
eigenvalues and do not affect the eigenvectors, $V_{ii}=-2d$ where $d$ is the
lattice dimension, $V_{ij}$ becomes a discrete finite difference Laplacian on
the lattice. It is well known that both the discrete and continuous Laplacian
have harmonic eigenfunctions, e.g.
$f_{\vec{k}}(\vec{x})=C_{\vec{k}}\prod_{i}^{d}\cos{(2\pi
k_{i}x_{i}/L+\phi_{i})}$ which is an orthogonal eigenbasis in $d$ dimensions
when $\vec{k}$ goes over all distances on the reciprocal lattice, $L$ is the
linear size of the lattice, $\phi_{i}=\pm\pi/4$ and $C_{\vec{k}}$ is an
appropriate normalization constant. We have thus shown that all interaction
matrices have a Fourier eigenbasis. An alternative, more direct but for our
purposes less illustrative, argument for the common Fourier basis is to note
that the structure of $V$ implies that it is a so called circulant matrix
Davis (1994), for which the result is known in the signal processing
literature. That the Fourier base effectively diagonalize the Hamiltonian in
the spherical model with translationally invariant interactions has also been
pointed out by Nussinov Nussinov .

This result has two important consequences. First, it helps us to understand
why systems with different interactions are expected to have similar ground
states. Second, knowledge of the universal eigenbasis allows us to compute the
energy spectrum for any particular system using a linear transform of the
potential. From (2) it follows that the energy per spin for a harmonic
eigenfunction with wave vector $\vec{k}$ are given by
$E(\vec{k})=\frac{1}{N}\sum_{\vec{x},\vec{y}}V(\left|\vec{x}-\vec{y}\right|)f_{\vec{k}}(\vec{x})f_{\vec{k}}(\vec{y})$.
Using various trigonometric identities and the radial structure of $V$ this
expression can be reduced to a radial Fourier transform

$E(\vec{k})=\sum_{\vec{r}}V(\left|\vec{r}\right|)\prod_{i=1}^{d}\cos{(2\pi
k_{i}r_{i}/L)}$ (4)

where the sum goes over all distances $\vec{r}$ on the lattice. Note that the
energy of a configuration can be computed through a fast Fourier transform
over the lattice.

Figure 2:  Examples of eigenmodes of the two dimensional spherical model. a-c,
Any phase shift of a ground state is a new ground state, so checkerboards,
stripes and everything between can be produced by the same model. Linear
combinations of $f_{(2,2)}$ with different phase shifts are shown, all having
the same energy. d, Exchanging the elements of $\vec{k}$ gives a new ground
state and linear combinations of them give rise to complex morphologies. Shown
is $\frac{1}{2}f_{(3,4)}+\frac{1}{2}f_{(4,3)}$.

The ground state of the spherical model is the eigenvector $f_{\vec{k}}$
corresponding to the minimum of $E(\vec{k})$. The simplest ground state
patterns in two dimensions are checkerboards and stripes with the
corresponding wavelength, exemplified in Fig. 2a and c. Further, the subspace
of the eigenbasis corresponding to the minimum can contain two kinds of
degeneracies. First, any change of the phases $\phi_{i}$ leaves the energy
invariant. In two dimensions this means that anything between checkerboards
and stripes can be produced, as illustrated in Fig. 2a-c. Second, the energy
is similarly unaffected by arbitrary permutations of the elements of
$\vec{k}$, reflecting that the energy only depends on the magnitude of the
wave vector (seen most clearly in the continuous limit (5)). Linear
combinations of vectors with different permutations give rise to complex
morphologies, exemplified in Fig. 2d.

Figure 3:  Predicting ground states in different two dimensional Ising models.
(Top) Potentials (inset) and their energy spectra in $|\vec{k}|$-space from
the transform (4). a, a purely repulsive potential Malescio and Pellicane
(2003), b, two competing interactions Löw  _et al._ (1994) (see equation (1)),
and c, an RKKY-like interaction Fischer and Hertz (1993). For small $k$ (long
wavelengths) the spectra only depend on the magnitude of $\vec{k}$, but for
large $k$ (short wavelengths) lattice effects breaks the independence on the
direction of $\vec{k}$ and $E(|\vec{k}|)$ becomes multivalued. To illustrate
this, we link series with constant $k_{y}$ with lines and in a examples of
such states are shown. (Below) Local minima as arrived by through Monte Carlo
annealing as well as ground states of Ising spin-1/2 models with corresponding
potentials.

The eigenmode analysis is exact for the spherical model but also has
implications for the discrete Ising models. It is not directly applicable as
in general an eigenvector of the interaction matrix cannot be constructed in
the restricted discrete space of the Ising spins. However, continuous
eigenvectors are often used to approximate solutions to discrete optimization
problems, for example graph coloring Aspvall and Gilbert (1984) and
partitioning networks into modules with minimal intra-connectivity Fiedler
(1973); Newman (2006). Here we use the same strategy to predict ground states
for Ising spin-$1/2$ models with corresponding potentials by mapping the spins
in the spherical model to $-1$ or $+1$ depending on their sign:
$\hat{f}_{\vec{k}}(\vec{x})=\mbox{sign}(f_{\vec{k}}(\vec{x}))$. 111This
discretization of the Fourier modes to approximate the ground states of Ising
models was also discussed in Nussinov  where it is argued that striped ground
states should be energetically favored if the minimum in the energy spectrum
is sharp. The argument we present does not have this assumption which is
important as many commonly used potentials have broad energy minima and
striped ground states, see e.g. Fig. 3 a and b. The discretization breaks the
energy degeneracy and stripes become energetically favorable compared to
checkerboards and more complex patterns. To see why we note that in each group
of degenerate eigenmodes, with wavelength $|\vec{k}|$, there exist linear
combinations that produce stripes, for example $\cos(\vec{k}\cdot\vec{x})$.
The error introduced by the discretization,
$\|\hat{f}_{\vec{k}}(\vec{x})-f_{\vec{k}}(\vec{x})\|_{2}$ with the standard
$L^{2}$ norm, always increases the energy in the discrete configurations when
compared to the continuous ground state. Due to the $\pm$-symmetry of the
harmonic functions, the difference between $\hat{f}_{\vec{k}}(\vec{x})$ and
$f_{\vec{k}}(\vec{x})$ (with appropriate scaling) is largest in regions where
the continuous function is close to $0$, i.e. at the interface between $+$ and
$-$ regions. From this argument it follows that the errortends to be smallest
for the striped eigenmode since the interface is minimized (assuming that the
width of the stripes is large compared to the lattice spacing).

There are two exceptions when the ground state does not have stripes. Energy
spectra with minimum at the boundary produce ground states that are either a
uniform ferromagnet (the zero frequency mode) or a checkerboard pattern (the
highest frequency mode allowed on the lattice) associated with an anti-
ferromagnet. These two cases can be viewed as degenerate cases of stripes with
infinite respective infinitesimal width.

In Fig. 3 some examples of Ising spin-$1/2$ models with different potentials
are shown together with their energy spectra in $|\vec{k}|$-space, examples of
local minima 222Local minima of the spin systems are found by a standard
Metropolis algorithm on lattices with periodic boundaries. The results shown
in Fig. 3 where generated by simulating 50$\times$50-lattices (in a and b) or
a 150$\times$150-lattice (in c) with the potentials shown at $T=0.1$ for
$2.5\times 10^{5}$ trial flips and then with $T=0$ until convergence. In Fig.
4 a 300$\times$200-lattice with $6\times 10^{5}$ trial flips at $T=0.2$ was
used. and their ground states. First is a short-ranged, purely repulsive
potential related to the model studied in Malescio and Pellicane (2003).
Second is a nearest neighbor ferromagnet with long-range repulsive Coulomb
interaction on the form (1) from Löw  _et al._ (1994). Last is an attenuated
Bessel function, $J_{0}(r)/(r+1)$, chosen for its similarity to the RKKY
interaction in spin glasses Fischer and Hertz (1993). We see that, while the
potentials are qualitatively very different, the ground states are defined
only by the minima in the energy spectrum, i.e. by a single length scale.
Through rescaling, the potentials can be adjusted to have identical ground
states. This illustrates how little observing striped behavior tells us about
the interactions in a system. The local minima do however show a qualitative
difference between the potentials in a and b and the RKKY-like potential in c,
probably related to the difference in localization in energy space.

Equation (4) also has implications for molecular self-assembly. The Fourier
basis in the transformation is orthogonal and can be inverted to find the
potential corresponding to a given energy spectrum. This allows us to design,
from an observed striped state, families of potentials that generate similar
patterns at low-temperature by identifying the dominant wavelength and invert
an energy spectrum with a minimum at this wavelength. A demonstration of the
procedure is shown in Fig. 4: the Fourier power spectrum of a pixelised image
of a metastable state in an experimental system Rosensweig  _et al._ (1983)
was calculated; an energy spectrum was constructed with gaussian minimum at
the same wavelength as the experimental system; and finally the corresponding
potential was found using the inverse transform of (4). The constructed system
has striped metastable states similar to those found in the experimental
system. We conclude that it is relatively easy to construct families of
potentials with desired metastable striped morphologies.

Figure 4:  Designing interactions to imitate observed striped patterns. a,
Stripes in a ferrofluid confined between two glass plates in a magnetic field
[from Rosensweig  _et al._ (1983); Seul and Andelman (1995), reprinted with
permission from AAAS and Elsevier] b, the (negative) radial power spectrum
(blue) of previous picture together with a gaussian (red), (inset) a pair
potential $V(r)$ constructed as the (inverse) transform (4) of said gaussian
and c, metastable state of an Ising spin-$1/2$ model with potential $V(r)$.
Note that the chosen energy spectrum is not unique. Many potentials having a
spectrum minimized at the same wavelength will show similar low-temperature
behavior.

In the continuous limit the transformation (4) becomes a Hankel transform, in
two dimensions defined as

$E(\vec{k})=2\pi\int_{0}^{\infty}rdrV(r)J_{0}(2\pi|\vec{k}|r)$ (5)

where $J_{0}$ is a Bessel function of the first kind. For the general
expression in higher dimensions, see Folland (1992). Note that in the
continuum limit the energy only depends on the magnitude of the wavevector
since the effects of the principal lattice directions disappear. As noted in
Fig. 3, this independence holds true on the lattice as well for small
wavevectors and sufficiently long-range interactions. Equation (5) allows us
to use the analytical properties of the Hankel transform to qualitatively
understand for example why the Bessel function of Fig. 3c has such a sharp
spectrum: the Hankel transform of a Bessel function is a Dirac delta function.

In summary we have shown that the energy spectrum of spherical spin systems
with isotropic interactions can be derived directly from the Fourier transform
of the potential. Due to a degeneracy in the energy eigenstates the spherical
model has ground states with various patterns such as stripes, checkerboards,
and more complicated morphologies. In discrete spin models the degeneracy is
broken leading to striped ground states being energetically favored. We
suggest that this can offer a generic explanation to why striped patterns are
so frequently observed in various experimental and natural systems.

###### Acknowledgements.

The authors would like to thank Olle Häggström for pointing out how purely
repulsive potentials can give rise to striped ground states.

## References

  * Seul and Andelman (1995) M. Seul and D. Andelman, Science 267, 476 (1995).
  * Bowman and Newell (1998) C. Bowman and A. C. Newell, Rev. Mod. Phys. 70, 289 (1998).
  * Emery  _et al._ (1999) V. Emery, S. Kivelson, and J. Tranquada, Proc. Natl. Acad. Sci. USA 96, 8814 (1999).
  * Keller and McConnell (1999) S. L. Keller and H. M. McConnell, Phys. Rev. Lett. 82, 1602 (1999).
  * Kern  _et al._ (1991) K. Kern, H. Niehus, A. Schatz, P. Zeppenfeld, J. Goerge, and G. Comsa, Phys. Rev. Lett. 67, 855 (1991).
  * Rosensweig  _et al._ (1983) R. E. Rosensweig, M. Zahn, and R. Shumovich, J. Magn. Magn. Mater. 39, 127 (1983).
  * Seul and Wolfe (1992) M. Seul and R. Wolfe, Phys. Rev. A 46, 7519 (1992).
  * Tranquada  _et al._ (1995) J. M. Tranquada, B. J. Sternlieb, J. D. Axe, Y. Nakamura, and S. Uchida, Nature 375, 561 (1995).
  * Orenstein and Millis (2000) J. Orenstein and A. Millis, Science 288, 468 (2000).
  * De’Bell  _et al._ (2000) K. De’Bell, A. B. MacIsaac, and J. P. Whitehead, Rev. Mod. Phys. 72, 225 (2000).
  * Garel and Doniach (1982) T. Garel and S. Doniach, Phys. Rev. B 26, 325 (1982).
  * Löw  _et al._ (1994) U. Löw, V. J. Emery, K. Fabricius, and S. A. Kivelson, Phys. Rev. Lett. 72, 1918 (1994).
  * Grousson  _et al._ (2000) M. Grousson, G. Tarjus, and P. Viot, Phys. Rev. E 62, 7781 (2000).
  * Stoycheva and Singer (2000) A. D. Stoycheva and S. J. Singer, Phys. Rev. Lett. 84, 4657 (2000).
  * Giuliani  _et al._ (2006) A. Giuliani, J. L. Lebowitz, and E. H. Lieb, Phys. Rev. B 74, 064420 (2006).
  * Weis and McConnell (1984) R. M. Weis and H. M. McConnell, Nature 310, 47 (1984).
  * Zaanen (1998) J. Zaanen, J. Phys. Chem. Solids 59, 1769 (1998).
  * Malescio and Pellicane (2003) G. Malescio and G. Pellicane, Nat. Mater. 2, 97 (2003).
  * Baxter (1982) R. J. Baxter,  _Exactly Solved Models in Statistical Mechanics_ (Academic Press, London, 1982).
  * Lay (2003) D. C. Lay, “Linear algebra and its applications,” (Pearson Education, Boston, 2003) Chap. 7.2.
  * Horn and Johnson (1990) R. A. Horn and C. R. Johnson,  _Matrix Analysis_ (Cambridge University Press, Cambridge, 1990).
  * Davis (1994) P. J. Davis,  _Circulant Matrices_ (Chelsea Pub Co, Providence, RI., 1994).
  * (23) Z. Nussinov, arxiv:cond-mat/0105253 .
  * Fischer and Hertz (1993) K. H. Fischer and J. A. Hertz,  _Spin Glasses_ (Cambridge University Press, Cambridge, UK, 1993).
  * Aspvall and Gilbert (1984) B. Aspvall and J. Gilbert, SIAM J. Algebra. Disc. 5, 526 (1984).
  * Fiedler (1973) M. Fiedler, Czech. Math. J. 23, 298 (1973).
  * Newman (2006) M. Newman, Proc. Natl. Acad. Sci. USA 103, 8577 (2006).
  * Note (1) This discretization of the Fourier modes to approximate the ground states of Ising models was also discussed in Nussinov  where it is argued that striped ground states should be energetically favored if the minimum in the energy spectrum is sharp. The argument we present does not have this assumption which is important as many commonly used potentials have broad energy minima and striped ground states, see e.g. Fig. 3 a and b.
  * Note (2) Local minima of the spin systems are found by a standard Metropolis algorithm on lattices with periodic boundaries. The results shown in Fig. 3 where generated by simulating 50$\times$50-lattices (in a and b) or a 150$\times$150-lattice (in c) with the potentials shown at $T=0.1$ for $2.5\times 10^{5}$ trial flips and then with $T=0$ until convergence. In Fig. 4 a 300$\times$200-lattice with $6\times 10^{5}$ trial flips at $T=0.2$ was used.
  * Folland (1992) G. B. Folland, “Fourier Analysis and Its Applications,” (Brooks/Cole Pub Co, Pacific Grove, 1992) Chap. 7.5.


# Gradient Descent Optimizes Infinite-Depth ReLU Implicit Networks with Linear
Widths

Tianxiang Gao
Department of Computer Science
Iowa State University
<EMAIL_ADDRESS>
&Hongyang Gao
Department of Computer Science
Iowa State University
<EMAIL_ADDRESS>

###### Abstract

Implicit deep learning has recently become popular in the machine learning
community since these implicit models can achieve competitive performance with
state-of-the-art deep networks while using significantly less memory and
computational resources. However, our theoretical understanding of when and
how first-order methods such as gradient descent (GD) converge on nonlinear
implicit networks is limited. Although this type of problem has been studied
in standard feed-forward networks, the case of implicit models is still
intriguing because implicit networks have infinitely many layers. The
corresponding equilibrium equation probably admits no or multiple solutions
during training. This paper studies the convergence of both gradient flow (GF)
and gradient descent for nonlinear ReLU activated implicit networks. To deal
with the well-posedness problem, we introduce a fixed scalar to scale the
weight matrix of the implicit layer and show that there exists a small enough
scaling constant, keeping the equilibrium equation well-posed throughout
training. As a result, we prove that both GF and GD converge to a global
minimum at a linear rate if the width $m$ of the implicit network is linear in
the sample size $N$, i.e., $m=\Omega(N)$.

##  1 Introduction

Recently, implicit neural networks attracts increasing attention in the
machine learning community, which have achieved competitive or dominated
performances of traditional neural networks in various domains such as
sequence modeling [4] with significantly less usage of computational resources
[8, 9, 3]. In implicit neural networks, the feature vectors are not created
recursively as traditional neural networks but provided implicitly through a
solution of an equilibrium equation. Implicit neural networks generalize the
recursive rules of many commonly used neural network architectures such as
feed-forward, convolution, residual, and recurrent networks [4, 12, 5].
However, the theoretical understanding of when and how a simple first-order
method such as gradient descent (GD) works for these implicit models is
limited, though this type of convergence problem has been well studied in
standard feed-forward networks [10, 1, 27, 21, 2, 23, 20]. Since implicit
networks can have infinitely many layers, the equilibrium equation is not
necessarily well-posed as it may admit zero or multiple solutions during
training. For example, [7, 4, 6, 18] all observe instability of forward
propagations in implicit models. Specifically, the number of iterations
required for forward propagation to find equilibrium points grows with
training epochs. Forward propagation likely becomes divergent as training goes
longer. A line of recent works make efforts to handle this well-posedness
challenge. For example, [12] reformulates the training problem in a so-called
Fenchel divergence formulation and use the projected gradient descent method
to solve the relaxed optimization problem; [25] formulates a splitting problem
for the forward propagation of implicit neural networks and use a proximal
operation to find the fixed point; [6] proposes a special regularization to
ensure the well-posedness. However, none of these works can theoretically
ensure the convergence of a gradient-based method.

Under some simplified setups, a line of recent works have tried to study this
convergence problem from the mathematical theory perspective. For instance,
[18] studies this problem for implicit models with linear activation function.
By applying an extra softmax layer on the shared weight matrix, the well-
posedness challenge is resolved. As a result, they are able to establish the
global linear convergence for gradient flow (GF). Unfortunately, their results
cannot be extended to nonlinear activation, especially for the nonsmooth ReLU
activation, which are critical to the learnability of deep neural networks.
Recently, [13] facility the training process by using a skip connection to the
output and introduces a scaling factor to scale the shared weight matrix of
the network, which can show there exists a small enough scaling constant which
keeps the forward propagation well-posed throughout training. To show the
convergence of ReLU-activated implicit neural networks, they take a neural
tangent kernel (NTK) method [17], where the dynamic of the network prediction
is governed by a Gram matrix that remains positive definite during training.
As a result, they successfully establish the global linear convergence for
both GF and GD as long as the width $m$ of the network is quadratic in the
sample size $N$, i.e., $m=\Omega(N^{2})$. However, their results can only be
applied to restricted range of implicit networks due to the special choice of
the output layer.

Main contribution. In this paper, we propose to establish the global
convergence results for implicit neural networks with the nonlinear ReLU
activation function and regular output layer. Specifically, we provide
sufficient conditions for the initialization under which GF and GD are
guaranteed to converge to a global minimum at a linear rate. Then, we show
that all these initial assumptions can be satisfied by a subset of
initialization where the network has linear width in the sample size, i.e.,
$m=\Omega(N)$. For popular random initialization, we show that these initial
conditions are satisfied with a high probability (w.h.p.) if the width of a
network is quadratic of the sample size, i.e., $m=\Omega(N^{2})$. Although
these results with similar order of overparameterization have been obtained
for finite-depth feed-forward neural networks [16, 21, 20], it is worth noting
that this is the first time, to our best knowledge, that such results are
provided for ReLU-activated implicit networks, which could have infinitely
many layers.

##  2 Preliminaries of Implicit Deep Learning

Notation: For a vector $x$, we use $\|x\|$ to denote its Euclidean norm. For a
matrix $A$, $\|A\|$ is its operator norm, and $\sigma_{\min}(A)$ and
$\sigma_{\max}(A)$ denote its smallest and largest singular values,
respectively. If $A$ is a square matrix, then $\lambda_{\min}(A)$ and
$\lambda_{\max}(A)$ denote the smallest and largest eigenvalue of $A$,
respectively. We use $\text{vec}\left(A\right)$ to denote the vectorization
operation applied on the matrix $A$. Given a function $Y=f(X)$, the derivative
$\partial f/\partial X$ is defined by
$\text{vec}\left(dY\right)=\left(\partial f/\partial X\right)^{T}dX$, where
$X$ and $Y$ can be scalars, vectors, and matrices. We also denote
$[N]:=\\{1,2,\cdots,N\\}$.

Let $X\in{\mbox{$\mathbb{R}$}}^{N\times d}$ and
$y\in{\mbox{$\mathbb{R}$}}^{N}$ be the training data. The implicit neural
network we consider in this paper has the transition at the $\ell$-th layer in
the following form

$\displaystyle Z^{\ell}=\sigma(\gamma Z^{\ell-1}A+\Phi),\quad\ell\geq 1,$ (1)

where $Z^{\ell}\in{\mbox{$\mathbb{R}$}}^{N\times m}$ is the output of the
$\ell$-layer with $Z^{0}=0$, the $i$-th row of
$\Phi\in{\mbox{$\mathbb{R}$}}^{N\times m}$ is the feature vector $\phi(x)$
under the feature map
$\phi:{\mbox{$\mathbb{R}$}}^{d}\rightarrow{\mbox{$\mathbb{R}$}}^{m}$,
$A\in{\mbox{$\mathbb{R}$}}^{m\times m}$ is the weight matrix shared among all
implicit layers, $\sigma(u)=\max\\{0,u\\}$ is the ReLU activation function,
and $\gamma\in(0,1)$ is a fixed scalar. As will be shown in Section 3, the
choice of $\gamma$ is essential to ensure the existence of the limit
$Z=\lim_{\ell\rightarrow\infty}Z^{\ell}$. As $\ell\rightarrow\infty$, an
implicit neural network can be considered as a neural network with infinitely
many layers. Consequently, $Z$ is not only the limit of the sequence
$\\{Z^{\ell}\\}_{\ell=1}^{\infty}$ but also an equilibrium point or fixed
point of the following equilibrium equation:

$\displaystyle Z=\sigma(\gamma ZA+\Phi).$ (2)

In general, the feature map $\phi$ is a nonlinear function, which extracts the
hidden features from the low-dimensional input. In this paper, we consider the
feature map $\phi$ as a simple one layer neural network that is activated by
ReLU activation function, i.e.,

$\displaystyle\Phi:=\phi(X)=\sigma(XW),$ (3)

where $W\in{\mbox{$\mathbb{R}$}}^{d\times m}$ is also a trainable weight
matrix. The training loss is given by

$\displaystyle L(\theta)=\frac{1}{2}\|\hat{y}-y\|^{2},$ (4)

where $\hat{y}=Zb$ is the prediction with weight vector
$b\in{\mbox{$\mathbb{R}$}}^{m}$ and $\theta:=\text{vec}\left(W,A,b\right)$ is
the collection of all trainable parameters.

##  3 Well-Posedness of the Forward and Backward Propagation

In this section, we establish sufficient conditions for the equilibrium
equation Eq.(2) to be well-posed in the sense that the existence of the
equilibrium point $Z$ is uniquely determined. Provided the existence of the
equilibrium point $Z$, we can derive the gradients of the parameters by using
the implicit function theorem instead of back-propagating all intermediate
layers. Previous work [13] has shown the existence of the unique equilibrium
point in the vector-based equilibrium equation Eq. (2) if the scalar
$\gamma>0$ is chosen small enough. We extend this result to a general matrix-
based mapping.

######  Lemma 3.1.

Suppose $\|A\|\leq M$ for some constant $M>0$ and choose the scalar $\gamma>0$
such that $\gamma_{0}:=\gamma M<1$. Then the existence of the fixed point $Z$
is uniquely determined. Moreover, we have
$\|Z^{\ell}\|_{F}\leq\frac{1}{1-\gamma_{0}}\|\Phi\|_{F}$ for all $\ell$, hence
$\|Z\|_{F}\leq\frac{1}{1-\gamma_{0}}\|\Phi\|_{F}$.

Lemma 3.1 shows that the transition Eq.(1) is a contraction mapping if a small
enough scalar $\gamma$ is selected for which $\gamma\|A\|<1$. The proof is
provided in Appendix A.1. However, the operator norm of $A(k)$ changes
throughout the training process. In general, a fixed scalar $\gamma$ cannot
guarantee the well-posedness during the entire training. That is the main
reason [7, 4, 6, 18] observe the number of forward iterations required to find
a fixed point gradually grows with training epochs. Thus, a simple strategy is
to find a sequence $\\{\gamma_{k}\\}_{k=1}^{\infty}$ of appropriate scalars
for each iteration that ensure the equilibrium equation are kept well-posed
over all iterations. However, finding an appropriate scalars at each iteration
is computationally expensive, since it needs to compute quantities related to
the operator norm of the weight matrix $A(k)$ (e.g., [12, 6]). Fortunately, we
can show $\|A(k)\|$ is uniformly upper bounded by some constant. As a result,
there indeed is a small constant for which the forward pass is kept well-posed
throughout training.

A finite-depth neural network needs to store all intermediate parameters and
apply backpropagation to compute the gradients of each weight matrix or
vector. Surprisingly, this expensive computation is not necessary for implicit
neural networks. Instead, the implicit function theorem provides an efficient
way to derive the gradients since the fixed point $Z$ is a root of the
function $f$ given by

$\displaystyle f(Z,A,W):=Z-\sigma(\gamma ZA+\Phi).$ (5)

We can easily show that the partial derivative $\partial f/\partial Z$ is
invertible, provided $\|A\|<\gamma^{-1}$. As a result, we obtain the partial
derivatives $\partial Z/\partial A$ and $\partial Z/\partial W$ by applying
the implicit function theorem. Then the gradients are obtained by using the
chain rule. Appendix A.2 includes the detailed derivation.

######  Lemma 3.2.

Suppose $\|A\|\leq M$ for some constant $M>0$ and choose the scalar $\gamma>0$
such that $\gamma_{0}:=\gamma M<1$. Then

$\displaystyle\lambda_{\min}\\{I_{Nm}-\gamma D(A^{T}\otimes
I_{N})\\}>1-\gamma_{0}>0.$ (6)

Hence the matrix $Q:=I_{Nm}-\gamma D(A^{T}\otimes I_{N})$ is invertible, and

$\displaystyle\frac{\partial L}{\partial W}=$
$\displaystyle\left[DE\left(I_{m}\otimes
X\right)\right]^{T}Q^{-T}\left(b^{T}\otimes I_{N}\right)^{T}(\hat{y}-y),$ (7)
$\displaystyle\frac{\partial L}{\partial A}=$
$\displaystyle\gamma\left[D(I_{m}\otimes Z)\right]^{T}Q^{-T}\left(b^{T}\otimes
I_{N}\right)^{T}(\hat{y}-y),$ (8) $\displaystyle\frac{\partial L}{\partial
b}=$ $\displaystyle Z^{T}(\hat{y}-y),$ (9)

where $D:=\mathop{\bf diag}[\text{vec}\left(\sigma^{\prime}(\gamma
ZA+\Phi)\right)]$, and $E:=\mathop{\bf
diag}\left[\text{vec}\left(\sigma^{\prime}(XW)\right)\right]$.

##  4 Main results

In this section, we first study the dynamics of the prediction $\hat{y}(t)$
that is induced by the gradient flow. It can be shown that a time-variant Gram
matrix controls the dynamics of the prediction. At the same time, the spectral
property of the Gram matrix is consistent throughout the training as long as
the network is over-parameterized. Based on the findings in gradient flow
analysis, we show that gradient descent with a fixed step size converges to a
global minimum of the implicit network at a linear rate.

###  4.1 Continuous time analysis: convergence of gradient flow

The gradient flow is given by $\frac{d\theta}{dt}=-\frac{\partial
L(t)}{\partial\theta}$, where $L(t):=L(\theta(t))$ is the corresponding loss
function for $\theta(t)$ at time $t$. By using the chain rule, we derive the
dynamics of the prediction $\hat{y}(t)$ in the following lemma, and the
derivation is deferred in Appendix A.3.

######  Lemma 4.1.

Assume $\|A(t)\|\leq M$ for all $t\geq 0$ and choose the scalar $\gamma>0$
small enough such that $\gamma_{0}:=\gamma M<1$. Then the dynamics of the
prediction $\hat{y}(t)$ is given by

$\displaystyle\frac{d\hat{y}}{dt}=-H(t)(\hat{y}(t)-y),$ (10)

where

$\displaystyle H(t):=$ $\displaystyle
Z(t)Z(t)^{T}+M(t)M(t)^{T}+\Pi(t)\Pi(t)^{T},$ (11) $\displaystyle M(t):=$
$\displaystyle\gamma[b(t)^{T}\otimes I_{N}]Q(t)^{-1}D(t)[I_{m}\otimes Z(t)],$
(12) $\displaystyle\Pi(t):=$ $\displaystyle\left[b(t)^{T}\otimes
I_{N}\right]Q(t)^{-1}D(t)E(t)\left[I_{m}\otimes X\right].$ (13)

Clearly, the matrix $H(t)$ is positive semidefinite. If there exists a
strictly positive scalar $\lambda_{0}>0$ for which
$\lambda_{\min}(H(t))\geq\lambda_{0}$ for all $t\geq 0$, then $L(t)$
consistently decreases to zero at a linear rate, i.e.,
$L(t)\leq\exp\left\\{-\lambda_{0}t\right\\}L(0)$. Thus, the problem is reduced
to show that the smallest singular value of at least one of the matrices
$Z(t)$, $M(t)$, and $\Pi(t)$ is lower bounded throughout the training. By
using simple matrix analysis results, however, we obtain the following
inequalities:

$\displaystyle\sigma_{\min}\left[M(t)\right]\geq$
$\displaystyle\frac{\gamma}{1-\gamma_{0}}\|b(t)\|\min_{i}\\{D_{ii}\\}\sigma_{\min}\left[Z(t)\right],$
$\displaystyle\sigma_{\min}\left[\Pi(t)\right]\geq$
$\displaystyle\frac{1}{1-\gamma_{0}}\|b(t)\|\min_{i}\\{D_{ii}\\}\min_{i}\\{E_{ii}\\}\sigma_{\min}\left(X\right).$

Thus, to lower bound the singular values of matrices $M(t)$ and $\Pi(t)$, one
must make extra assumptions on the data sample $X$ and activation $\sigma$. An
example method is introduced in [21], where the neural network has to follow a
pyramidal structure and the activation function has to be sufficiently smooth.
This method does not hold for ReLU activation due to the non-smoothness. In
contrast, we adopt the method proposed by [20] for finite-depth ReLU neural
network, where their analysis focuses on the evolution of the last layer of
the network, that is, $Z(t)$ in the implicit neural network. With appropriate
assumptions on the initial conditions, we can establish the global convergence
result for the gradient flow in Theorem 4.2, and the entire proof is provided
in Appendix A.4.

######  Theorem 4.2.

Let $C_{1},C_{2},C_{3}>0$ be given positive numbers. Denote
$\alpha_{0}:=\sigma_{\min}(Z(0))$, $\lambda_{1}:=\|W(0)\|+C_{1}$,
$\lambda_{2}:=\|A(0)\|+C_{2}$, and $\lambda_{3}:=\|b(0)\|+C_{3}$. Choose
$\gamma>0$ small enough for which $\gamma_{0}:=\gamma\lambda_{2}<1$. Assume
the following conditions are satisfied at initialization

$\displaystyle\alpha_{0}^{2}\geq$
$\displaystyle\frac{4}{1-\gamma_{0}}\lambda_{0}\|X\|_{F}\|\hat{y}(0)-y\|,$
(14) $\displaystyle\alpha_{0}^{3}\geq$
$\displaystyle\left[1+\frac{\gamma_{0}^{2}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}^{2}}{\lambda_{2}^{2}}\right]\frac{8\lambda_{3}}{(1-\gamma_{0})^{2}}\|X\|_{F}^{2}\|\hat{y}(0)-y\|,$
(15)

where
$\lambda_{0}:=\max\left\\{\frac{\lambda_{1}}{C_{3}},\frac{\gamma_{0}\lambda_{1}\lambda_{3}}{(1-\gamma_{0})C_{2}\lambda_{2}},\frac{\lambda_{3}}{C_{1}}\right\\}$.
Then for all $t\geq 0$ the followings hold:

  1. (i)

$\|W(t)\|\leq\lambda_{1}$, $\|A(t)\|\leq\lambda_{2}$,
$\|b(t)\|\leq\lambda_{3}$,

  2. (ii)

$\sigma_{\min}(Z(t))\geq\alpha_{0}/2$,

  3. (iii)

$L(t)\leq\exp\\{-\alpha_{0}^{2}t/2\\}L(0)$, where $L(t):=L(\theta(t))$.

Theorem 4.2 shows that the operator norm $A(t)$ is upper bounded by the
constant $\lambda_{2}$ throughout training, i.e., $\|A(t)\|\leq\lambda_{2}$.
By choosing $\gamma>0$ small enough, the forward propagation is kept well-
posed throughout training. Moreover, Theorem 4.2 also indicates that the
smallest singular value of $Z(t)$ is lower bounded by the constant
$\alpha_{0}/2$ during training. As a result, the training loss consistently
decreases to zero at a linear rate. It is worth noting that the network only
needs to have linear width $m$ in the sample size $N$ as long as the initial
conditions (14)-(15) are satisfied. Section 5 will provide concrete examples
to show that these initial conditions are indeed satisfied by using only
linear widths.

###  4.2 Discrete time analysis: convergence of gradient descent

By applying the Euler method to the gradient flow with stepsize $\eta>0$, we
obtain the gradient descent as follows
$\theta(k+1)=\theta(k)-\eta\frac{\partial L(k)}{\partial\theta}$. Unlike the
continuous analysis in the gradient flow, we don’t have the explicit formula
of the dynamics of the prediction $\hat{y}(k)$ in the discrete time analysis.
Instead, we need to first show the difference of equilibrium points $Z(k)$ in
two consecutive iterations. The following result serves this purpose by
providing a bound between two equilibrium points based on their parameters. We
defer the proof in Appendix A.5.

######  Lemma 4.3.

Given matrices $(W_{a},A_{a})$ and $(W_{b},A_{b})$, let
$\lambda_{1}:=\max\\{\|W_{a}\|,\|W_{b}\|\\}$, and
$\lambda_{2}:=\max\\{\|A_{a}\|,\|A_{b}\|\\}$. Choose $\gamma>0$ small enough
for which $\gamma_{0}:=\gamma\lambda_{2}<1$. Then the corresponding
equilibrium points exist and are denoted by $Z_{a}$ and $Z_{b}$, respectively.
Moreover, we have

$\displaystyle\|Z_{a}-Z_{b}\|\leq\frac{\|X\|_{F}}{1-\gamma_{0}}\left[\frac{\gamma_{0}}{1-\gamma_{0}}\frac{\lambda_{1}}{\lambda_{2}}\|A_{a}-A_{b}\|+\|W_{a}-W_{b}\|\right].$
(16)

Lemma 4.3 allows us to derive the relationship of the predictions $\hat{y}(k)$
between two consecutive iterations. The loss can be shown to consistently
decrease as long as the step size selected is small enough. To satisfy the
conditions of the step size, more initial conditions are needed to obtain the
convergence result for gradient descent. The convergence result of the
gradient descent is provided in the following theorem, and the proof is
included in Appendix A.6.

######  Theorem 4.4.

Let $C_{1},C_{2},C_{3}>0$ be given positive numbers. Denote
$\alpha_{0}:=\sigma_{\min}(Z(0))$, $\lambda_{1}:=\|W(0)\|+C_{1}$,
$\lambda_{2}:=\|A(0)\|+C_{2}$, and $\lambda_{3}:=\|b(0)\|+C_{3}$. Choose
$\gamma>0$ small enough for which $\gamma_{0}:=\gamma\lambda_{2}<1$. Assume
the following conditions are satisfied at initialization

$\displaystyle\alpha_{0}^{2}\geq$
$\displaystyle\frac{8}{1-\gamma_{0}}\lambda_{0}\|X\|_{F}\|\hat{y}(0)-y\|$ (17)
$\displaystyle\alpha_{0}^{3}\geq$
$\displaystyle\left[1+\frac{\gamma_{0}^{2}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}^{2}}{\lambda_{2}^{2}}\right]\frac{16\lambda_{3}}{(1-\gamma_{0})^{2}}\|X\|_{F}^{2}\|\hat{y}(0)-y\|$
(18) $\displaystyle\alpha_{0}^{2}\geq$
$\displaystyle\left[1+\frac{\gamma_{0}^{2}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}^{2}}{\lambda_{2}^{2}}\right]\frac{16\lambda_{3}^{2}}{(1-\gamma_{0})^{2}}\|X\|_{F}^{2}$
(19)

where $\lambda_{0}$ is defined in Theorem 4.2. Choose step size $\eta>0$ such
that

$\displaystyle\eta<\min\left\\{\frac{4}{\alpha_{0}^{2}},2(1-\gamma_{0})^{2}\overline{\lambda}\lambda_{1}^{-4}\lambda_{3}^{-2}\|X\|_{F}^{-2}\right\\},$
(20)

where
$\overline{\lambda}:=\left[1+\frac{\gamma_{0}^{2}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}^{2}}{\lambda_{2}^{2}}\right]\left(\left[1+\frac{\gamma_{0}^{2}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}^{2}}{\lambda_{2}^{2}}\right]\lambda_{1}^{-2}+\lambda_{3}^{-2}\right)^{-2}$.
Then the followings holds $\forall k\geq 0$

  1. (i)

$\|W(k)\|\leq\lambda_{1}$, $\|A(k)\|\leq\lambda_{2}$,
$\|b(k)\|\leq\lambda_{3}$,

  2. (ii)

$\sigma_{\min}(Z(k))\geq\alpha_{0}/2$,

  3. (iii)

$L(k)\leq(1-\eta\alpha_{0}^{2}/4)^{k}L(0)$, where $L(k):=L(\theta(k))$.

With some extra conditions on initialization and stepsize, Theorem 4.4 obtains
similar results for gradient descent as in Theorem 4.2 for gradient flow. In
particular, the well-posedness is guaranteed throughout training using a
singular constant $\gamma$. The smallest singular value $\sigma_{\min}(Z(t))$
remains greater than the constant $\alpha_{0}/2$ so that the loss $L(k)$
continuously decreases to zero at a linear rate. Similarly, these results only
need the implicit network with linear width in the sample size $N$. However,
these results are not directly comparable because it remains unclear how
likely the initial conditions (17)-(19) are satisfied. Thus, we reserve
Section 5 to further analyze on the satisfiability of the initial assumptions.
Specifically, we address this issue by providing concrete examples by which
the initial assumptions are indeed satisfied by only using linear widths.

##  5 Satisfiability of the initial assumptions

This section provides concrete examples for which the initial conditions
(17)-(19) of Theorem 4.4 are satisfied. Since the equilibrium point $Z(k)$ is
provided implicitly, the singular values of $Z(k)$ are difficult to determine,
even at initialization. What exacerbates the problem is the nonlinearity of
the ReLU activation $\sigma$. Fortunately, we can find a subset of
initialization that can easily induce the singular values of $Z$ by using
homogeneity and nonnegativity of $\sigma$. The following result provides an
explicit expression for the equilibrium point $Z$ in terms of the feature
matrix $\Phi$ as long as $A$ is appropriately initialized. The proof is given
in Appendix A.7.

######  Lemma 5.1.

Assume $A_{ij}\geq 0$ for all $i,j\in[m]$ and $\gamma>0$ is chosen for which
$\|A\|<\gamma^{-1}$. Then $Z=\Phi\left(I_{m}-\gamma A\right)^{-1}$.

Suppose that we are given $A(0)\in{\mbox{$\mathbb{R}$}}^{m\times m}$ with
$A_{ij}\geq 0$ for all $i,j\in[m]$ and $\|A(0)\|\leq\lambda_{2}$ for some
positive constant $\lambda_{2}$. By choosing $\gamma>0$ small enough for which
$\gamma_{0}=\gamma\lambda_{2}<1$, Lemma 5.1 implies that
$\alpha_{0}:=\sigma_{\min}[Z(0)]\geq\sigma_{\min}[\Phi(0)]/(1+\gamma_{0})$.
Therefore, the satisfiability of the assumption made in Theorem 4.4 is reduced
to find an appropriate $\Phi(0)$ or $W(0)$ for which
$\sigma_{\min}[\Phi(0)]/(1+\gamma_{0})$ satisfies the initial conditions
(17)-(19). Some mostly recent literature [26, 10, 21, 22, 20] have been
provided rigorously theoretical analyses and concrete examples to demonstrate
that the initial conditions are easily satisfied with different types of
initialization. The rest of this section provides concrete examples to show a
subset of initialization that satisfies these assumptions.

###  5.1 Linear width $m=\Omega(N)$ suffices the initial conditions

Specifically, we apply Theorem 4.4 for $C_{1}=C_{2}=C_{3}=1$. Let
$\theta(0):=\text{vec}\left(W(0),A(0),b(0)\right)$, where we choose $W(0)$
such that $\sigma_{\min}[\Phi(0)]>0$, $A(0)_{ij}\geq 0$ with
$\|A(0)\|=\|W(0)\|$, and $b(0)=0$. For $\sigma_{\min}[\Phi(0)]>0$, a concrete
example is that $W(0)$ is chosen according to LeCun’s initialization (e.g.,
see [20, Section 3.1],[21, Section 3.1]). It is easy to come up with examples
for which $A(0)_{ij}\geq 0$ and $\|A(0)\|=\mathcal{O}\left(\|W(0)\|\right)$.
For example, $A(0):=\|W(0)\|I_{m}$. Thus, for simplicity, we assume
$\|A(0)\|=\|W(0)\|$, then $\lambda_{2}=\lambda_{1}$. Since $b(0)=0$, we have
$\lambda_{3}=1$, and $\|\hat{y}(0)-y\|=\|y\|$. Next, we can choose $\gamma>0$
small enough for which $\gamma_{0}=1/2$. It follows from Lemma 5.1 that the
initial conditions (17)-(19) become

$\displaystyle\left(\frac{2}{3}\sigma_{\min}[\Phi(0)]\right)^{2}\geq$
$\displaystyle 16\lambda_{1}\|X\|_{F}\|y\|$ (21)
$\displaystyle\left(\frac{2}{3}\sigma_{\min}[\Phi(0)]\right)^{3}\geq$
$\displaystyle 128\|X\|_{F}^{2}\|y\|$ (22)
$\displaystyle\left(\frac{2}{3}\sigma_{\min}[\Phi(0)]\right)^{2}\geq$
$\displaystyle 128\|X\|_{F}^{2}.$ (23)

Let $\tilde{\theta}(0):=\beta\theta(0)$ for some $\beta>0$ and denote
$\alpha_{1}:=\sigma_{\min}[\Phi(0)]$. Condition (21) at $\tilde{\theta}(0)$
becomes

$\displaystyle\left(\frac{2}{3}\beta\alpha_{1}\right)^{2}\geq$ $\displaystyle
16(\beta\|W(0)\|+1)\|X\|_{F}\|y\|.$ (24)

The left-hand side (LHS) of the above inequality is a polynomial of degree $2$
in $\beta$, whereas the RHS is linear in $\beta$. Thus, as $\|W(0)\|$ is fixed
in the inequality above, the inequality (24) is satisfied as long as $\beta$
is chosen large enough. Similarly, the condition (22) and (23) are satisfied
for $\tilde{\theta}(0)$ if

$\displaystyle\left(\frac{2}{3}\beta\alpha_{1}\right)^{3}\geq$ $\displaystyle
128\|X\|_{F}^{2}\|y\|,$ (25)
$\displaystyle\left(\frac{2}{3}\beta\alpha_{1}\right)^{2}\geq$ $\displaystyle
128\|X\|_{F}^{2}.$ (26)

The RHS of the above inequalities are all constants, whereas the LHS are
polynomials of degree $3$ and $2$ in $\beta$, respectively. Thus, the
conditions (25) and (26) are also satisfied for large enough $\beta$. As a
result, all the initial conditions of Theorem 4.4 are satisfied at
$\tilde{\theta}(0)$ for large $\beta$.

###  5.2 Random initialization: width $m=\Omega(N^{2})$ suffices the initial
conditions

We show that conditions (17)-(19) are satisfied under similar random
initialization as LeCun’s Initialization [19], Xavier Initialization [14], and
He Initialization [15], provided a stronger condition $m=\Omega(N^{2})$. For
simplicity, we assume $\|x_{i}\|=1$ and
$\left|y_{i}\right|=\mathcal{O}\left(1\right)$ for all $i\in[N]$. Then
$\|X\|_{F}=\sqrt{N}$ and $\|y\|=\mathcal{O}\left(\sqrt{N}\right)$. Consider
initialization

$\displaystyle b(0)_{i}\overset{i.i.d.}{\sim}\mathcal{N}(0,1/m),\quad
A(0)_{ij}\overset{i.i.d.}{\sim}\left|\mathcal{N}\right|(0,1),\quad
W(0)_{ij}\overset{i.i.d.}{\sim}\mathcal{N}(0,1),$ (27)

where $\left|\mathcal{N}\right|$ stands for half-normal distribution. It
follows from Theorem 4.4.5 of [24] that (with high probability)
$\|b(0)\|=\mathcal{O}\left(1\right)$,
$\|W(0)\|=\mathcal{O}\left(\sqrt{m}\right)$, and
$\|A(0)\|=\mathcal{O}\left(\sqrt{m}\right)$. For $C_{1}=C_{2}=C_{3}=1$, we
have $\lambda_{1}=\mathcal{O}\left(\sqrt{m}\right)$,
$\lambda_{2}=\mathcal{O}\left(\sqrt{m}\right)$,
$\lambda_{3}=\mathcal{O}\left(1\right)$. We can choose $\gamma>0$ small enough
to ensure $\gamma_{0}=1/2$. By using standard concentration argument and Lemma
5.1, we obtain

$\displaystyle\|\hat{y}(0)\|=\|Z(0)b(0)\|=\mathcal{O}\left(\left\|\Phi(0)\left[I_{m}-\gamma
A(0)\right]^{-1}\right\|\right)=\mathcal{O}\left(\|X\|_{F}\right).$

Since $\|X\|_{F}=\sqrt{N}$ and $\|y\|=\mathcal{O}\left(\sqrt{N}\right)$, the
initial conditions (17)-(19) are reduced to

$\displaystyle\alpha_{1}^{2}=\Omega\left(N\sqrt{m}\right).$ (28)

Using Matrix-Chernoff inequality, one can easily show $\alpha_{1}^{2}\geq
m\lambda_{*}/4$ (see Lemma 5.2 of [22]) with a probability of at least
$1-\delta$, if $m=\tilde{\Omega}(N/\lambda_{*})$ holds, where
$\lambda_{*}=\lambda_{\min}\left(G_{*}\right)$ with
$G_{*}:=\mathbb{E}_{w\sim\mathcal{N}(0,I_{d})}\left[\sigma(Xw)\sigma(Xw)^{T}\right]$
and $\tilde{\Omega}$ omits logarithmic factors depending on $\delta$. It can
be shown $\lambda_{*}>0$ under some mild data assumption. For example,
$\lambda_{*}>0$ if no two data are parallel to each other by Lemma 3.2 of
[13]. Additionally, we have

$\displaystyle\lambda_{*}\leq\frac{\mathbf{tr}(G_{*})}{N}=\frac{\mathbb{E}\|\sigma(Xw)\|^{2}}{N}\leq\frac{\mathbb{E}\|Xw\|^{2}}{N}=\frac{\|X\|_{F}^{2}}{N}=1.$

Therefore, all the initial conditions are satisfied for
$m=\Omega(N^{2}\lambda_{*}^{-2})$. Then, if we additionally assume that the
data points follow some sub-Gaussian distribution, [21, Theorem 3.3] implies
that $\lambda_{*}=\Theta(1)$. Thus, all initial conditions are satisfied for
$m=\Omega(N^{2})$.

##  6 Related works

The convergence problem of gradient flow and gradient descent for standard
deep networks has been studied recently [10, 1, 27, 21, 2, 20]. To show the
convergence, the most common strategy is NTK method [17], where it can be
shown that dynamic of the prediction is governed by a gram matrix whose
smallest eigenvalues is kept strictly positive during training. The previous
work of [10, 1, 27, 2] requires all hidden layers have large widths since
their analyses rely on studying various quantities related to the changes in
the activation patterns during training. In [21], the authors claim one wide
hidden layer is enough if the network follows a pyramidal topology and the
derivative of activation is lower bounded. Moreover, they prove that linear
width is enough for overparameterization, while previous works need a width of
at least $\Omega(N^{8})$. The most recent study [20] further confirms linear
width suffices global convergence for ReLU activated networks by focusing the
analysis on the last hidden layer so that the pyramidal topology is not
necessary. However, none of the previous works can be applied to implicit
networks directly since implicit networks have infinitely many layers, and
equilibrium equations may not be well-posed during training.

The well-posedness of the forward propagation is the main challenge for any
infinite-depth or implicit-depth networks. A line of recent works [25, 6, 7,
4] have shown the instability of forward propagation in implicit networks, and
the probability of divergent forward pass gradually raises with training
epochs. Several strategies are suggested to deal with the well-posedness
problem, e.g., adding constraints [12] or regularization[6], reformulates the
forward propagation [25]. Unfortunately, none of them can theoretically
guarantee convergence. By using an extra softmax layer on the weight matrix,
[18] proves the convergence of the gradient flow, but the result only holds
for linear activation while nonlinearity is critical for learnability. To deal
with well-posedness, our strategy is similar to [13] by introducing a scalar
in front of the weight matrix. The well-posedness problem is resolved since we
can find a sufficiently small constant scalar. The novelty of this paper is
that implicit network considered in this paper has regular output while [13]
uses a skip connection to facilitate training process. Moreover, this is the
first work, to our best knowledge, that establishes the global convergence for
an infinite-depth (weight tying) network with linear width.

##  7 Experimental Results

Figure 1: We evaluate the impact of the width $m$ on the training loss, test
loss, and operator norm of the scaled matrix $\gamma A(k)$ on the modified
dataset of MNIST.

We evaluate our results use real-world datasets like MNST, FashionMNST,
CIFAR10, and SVHN. Experimental Setups. For each dataset, we use classes $0$
and $1$, and $500$ samples are randomly drawn from each class to generate the
training dataset of $N=1000$. All data samples are converted to gray scale and
resized to $28\times 28$. We also normalize each data to have unit norm. We
run $500$ epochs of gradient descent with a fixed step-size. It follows from
the analysis of Section 5 that we initialize
$W(0)_{ij}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,1/m)$. For simplicity, we
set $A(0)=I_{m}$ and $b(0)=0$, so that
$\lambda_{1}=\lambda_{2}=\lambda_{3}=\mathcal{O}\left(1\right)$. Let
$\theta(0):=\text{vec}\left(W(0),A(0),b(0)\right)$. Then we multiply
$\theta(0)$ by a large $\beta>0$ if $\theta(0)$ does not satisfies the initial
conditions (17)-(19). Moreover, it follows from Theorem 4.4 that we can set
$\eta=N^{-1}$, i.e., $\eta=10^{-3}$, to ensure consistent descent in loss
function.

###  7.1 Over-Parameterization Study

In this work, some of our finds are based on the setting of over-
parameterization. Thus, we first study how overparameterization affects the
convergence rates of the gradient descent method. In addition, we also study
the impact of overparameterized implicit networks on unseen test data.
Similarly, the test data is constructed by randomly selecting $500$ unseen
samples from each class. Third, we study how overparameterization affects the
changes of $\|A(k)\|$ from its initialization, since consistent $\|A(k)\|$ is
critical for the well-posedness of the forward propagation. Here, “operator
norm” in the plots denotes $\gamma\|A(k)\|$ with $\gamma=0.1$.

Due to the space limit, Figure 1 only includes the result of the dataset
MNIST, and the rest results for other real datasets are included in Appendix
B. The first sub-figure in Figure 1 shows that as $m$ becomes larger, better
convergence rates can be observed. Accordingly, the second sub-figure shows
that the neural networks can achieve lower test loss as $m$ becomes larger.
The third sub-figure shows that the operator norms are changed throughout the
training process, but overall the operator norms are approximately equal to
their initialization, i.e., $\gamma$, since we set $A(0)=I_{m}$.

###  7.2 Hyper-parameter Study of $\gamma$

Figure 2: We evaluate the impact of choice $\gamma$ on the training loss, test
loss, and operator norm of the scaled matrix $\gamma A(k)$ on the modified
dataset of MNIST.

In this work, we introduce the scalar $\gamma$ to solve the well-posedness of
implicit neural networks. Since $\gamma$ is critical for our theoretical
finds, it is worth to conduct a series of experiments to study the impact of
this new hyper-parameter. Recall that Theorem 4.4 shows $m=\Omega(N)$ is
enough to ensure the convergence of gradient descent. With $N=1000$, we set
$m=1000$ and keep the rest setups the same as before. We choose the values of
$\gamma$ from $\\{0.1,0.3,0.5,0.8\\}$ such that it covers a reasonable range
of values for $\gamma$. Here, we study how the different choices of hyper-
parameter $\gamma$ impact the training and test performances. We summarize the
results on MNIST in Figure 2. The rest results are included in Appendix B to
save the space. From Figure 2, we can see the lower training loss is obtained
by using relatively larger $\gamma$. Accordingly, the achieved test loss is
also lower when $\gamma$ becomes larger. These observations show that $\gamma$
has impact on training and testing performances.

$\gamma$ | 0.1 | 0.3 | 0.5 | 0.8
---|---|---|---|---
iteration $\\#$ | 6 | 9 | 15 | 47.5
Table 1: Impact of the choice $\gamma$ on forward iteration to converge on the
modified dataset of MNIST.

On the other hand, we also test how the choice of $\gamma$ impacts the number
of iteration forward propagation needed to find the equilibrium point. We stop
the forward propagation if either $\|Z^{\ell+1}-Z^{\ell}\|\leq 10^{-2}$ or
reaches the max iteration of $100$. Table 1 contains the averaged number of
forward iteration for different $\gamma$ values to find the equilibrium point
$Z$ on MNIST. Although Figure 2 shows the lower training and test losses can
be obtained by using relatively larger $\gamma$, Table 1 indicates larger
$\gamma$ results in more iterations needed in the forward propagation to find
the equilibrium point. Moreover, as the operator norm of $\gamma A(k)$ becomes
larger as growth of epoch, forward propagation takes even more iterations to
converge. This scenario is also observed in some previous works [7, 4, 6].
Thus, there is probably a trade-off between the convergence speed and well-
posedness of implicit neural networks.

###  7.3 Hyper-Parameter Study of Learning Rate

Figure 3: We evaluate the impact of the step size $\eta$ on the training loss,
test loss, and operator norm of the scaled matrix $\gamma A(k)$ on the
modified dataset of FashionMNIST.

In general, the step size or learning rate $\eta$ as a hyper-parameter is
significant on both training and test performance. However, the choice of
$\eta$ mostly likely is based on the practical experiments. Based on some
theoretical analyses, $\eta=\mathcal{O}\left(L^{-1}\right)$ is one of the
widely obtained result, where $L$ is the Lipschitz constant by assuming the
gradient of neural network is Lipschitz continuous. Unfortunately, the value
of $L$ is generally unknown or relatively large [21]. In Section 5, we show
$\eta=\mathcal{O}\left(N^{-1}\right)$ is enough to guarantee the convergence
of the gradient descent, which is much larger than many previous results
$\eta=\mathcal{O}\left(N^{-2}\right)$ [11, 10, 13], even the implicit neural
network could have infinitely many layers. Thus, we also construct a series of
numerical experiments to justify the choice of step size $\eta$. With almost
the same setup as before, we set the values for step size $\eta$ from
$\\{10^{-2},10^{-3},10^{-4}\\}$. The corresponding train losses and test
losses are illustrated in Figure 3. The first sub-figure in Figure 3 shows
that the objective value does not consistently decrease if we choose the step
size $\eta$ larger than the suggested value in Theorem 4.4. Accordingly, the
second sub-figure in Figure 3 shows instability of the test loss for large
step size. In addition, the operator norm of $\gamma A(k)$ raises fast to a
relatively large value, which probably causes divergence of the forward
propagation. Therefore, choosing step size $\eta$ larger than $N^{-1}$ is
probably causing the divergence of the forward propagation and the entire
training process. This experimentally indicates the result of the step size
$\eta$ in Theorem 4.4 is tight.

##  8 Conclusion

This paper studies the convergence problem of first-order methods such as
gradient descent for ReLU implicit networks with infinitely many layers.
Specifically, we provide sufficient conditions under which both gradient flow
and gradient descent converge to a global minimum at a linear rate. Moreover,
we show that these sufficient conditions can be indeed satisfied by some
initialization as long as the width $m$ is linear in the sample size $N$,
i.e., $m=\Omega(N)$, even when the implicit networks have infinitely many
layers and are activated by ReLU. Moreover, we also show that popular random
initializations satisfy sufficient conditions under a stronger condition where
the width is quadratic of the sample size, i.e., $m=\Omega(N^{2})$.

## References

  * [1] Zeyuan Allen-Zhu, Yuanzhi Li, and Zhao Song.  A convergence theory for deep learning via over-parameterization.  In International Conference on Machine Learning, pages 242–252. PMLR, 2019.
  * [2] Sanjeev Arora, Simon S Du, Wei Hu, Zhiyuan Li, Ruslan Salakhutdinov, and Ruosong Wang.  On exact computation with an infinitely wide neural net.  arXiv preprint arXiv:1904.11955, 2019.
  * [3] Shaojie Bai, J Zico Kolter, and Vladlen Koltun.  Trellis networks for sequence modeling.  arXiv preprint arXiv:1810.06682, 2018.
  * [4] Shaojie Bai, J Zico Kolter, and Vladlen Koltun.  Deep equilibrium models.  arXiv preprint arXiv:1909.01377, 2019.
  * [5] Shaojie Bai, Vladlen Koltun, and J Zico Kolter.  Multiscale deep equilibrium models.  arXiv preprint arXiv:2006.08656, 2020.
  * [6] Shaojie Bai, Vladlen Koltun, and J Zico Kolter.  Stabilizing equilibrium models by jacobian regularization.  arXiv preprint arXiv:2106.14342, 2021.
  * [7] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud.  Neural ordinary differential equations.  arXiv preprint arXiv:1806.07366, 2018.
  * [8] Raj Dabre and Atsushi Fujita.  Recurrent stacking of layers for compact neural machine translation models.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 6292–6299, 2019.
  * [9] Mostafa Dehghani, Stephan Gouws, Oriol Vinyals, Jakob Uszkoreit, and Łukasz Kaiser.  Universal transformers.  arXiv preprint arXiv:1807.03819, 2018.
  * [10] Simon Du, Jason Lee, Haochuan Li, Liwei Wang, and Xiyu Zhai.  Gradient descent finds global minima of deep neural networks.  In International Conference on Machine Learning, pages 1675–1685. PMLR, 2019.
  * [11] Simon S Du, Xiyu Zhai, Barnabas Poczos, and Aarti Singh.  Gradient descent provably optimizes over-parameterized neural networks.  arXiv preprint arXiv:1810.02054, 2018.
  * [12] Laurent El Ghaoui, Fangda Gu, Bertrand Travacca, Armin Askari, and Alicia Y Tsai.  Implicit deep learning.  arXiv preprint arXiv:1908.06315, 2, 2019.
  * [13] Tianxiang Gao, Hailiang Liu, Jia Liu, Hridesh Rajan, and Hongyang Gao.  A global convergence theory for deep relu implicit networks via over-parameterization.  arXiv preprint arXiv:2110.05645, 2021.
  * [14] Xavier Glorot and Yoshua Bengio.  Understanding the difficulty of training deep feedforward neural networks.  In Proceedings of the thirteenth international conference on artificial intelligence and statistics, pages 249–256. JMLR Workshop and Conference Proceedings, 2010.
  * [15] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.  In Proceedings of the IEEE international conference on computer vision, pages 1026–1034, 2015.
  * [16] Jiaoyang Huang and Horng-Tzer Yau.  Dynamics of deep neural networks and neural tangent hierarchy.  In International conference on machine learning, pages 4542–4551. PMLR, 2020.
  * [17] Arthur Jacot, Franck Gabriel, and Clément Hongler.  Neural tangent kernel: Convergence and generalization in neural networks.  arXiv preprint arXiv:1806.07572, 2018.
  * [18] Kenji Kawaguchi.  On the theory of implicit deep learning: Global convergence with implicit layers.  arXiv preprint arXiv:2102.07346, 2021.
  * [19] Yann A LeCun, Léon Bottou, Genevieve B Orr, and Klaus-Robert Müller.  Efficient backprop.  In Neural networks: Tricks of the trade, pages 9–48. Springer, 2012\.
  * [20] Quynh Nguyen.  On the proof of global convergence of gradient descent for deep relu networks with linear widths.  arXiv preprint arXiv:2101.09612, 2021.
  * [21] Quynh Nguyen and Marco Mondelli.  Global convergence of deep networks with one wide layer followed by pyramidal topology.  arXiv preprint arXiv:2002.07867, 2020.
  * [22] Quynh Nguyen, Marco Mondelli, and Guido F Montufar.  Tight bounds on the smallest eigenvalue of the neural tangent kernel for deep relu networks.  In International Conference on Machine Learning, pages 8119–8129. PMLR, 2021.
  * [23] Samet Oymak and Mahdi Soltanolkotabi.  Toward moderate overparameterization: Global convergence guarantees for training shallow neural networks.  IEEE Journal on Selected Areas in Information Theory, 1(1):84–105, 2020.
  * [24] Roman Vershynin.  High-dimensional probability: An introduction with applications in data science, volume 47.  Cambridge university press, 2018.
  * [25] Ezra Winston and J Zico Kolter.  Monotone operator equilibrium networks.  Advances in neural information processing systems, 33:10718–10728, 2020.
  * [26] Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals.  Understanding deep learning (still) requires rethinking generalization.  Communications of the ACM, 64(3):107–115, 2021.
  * [27] Difan Zou, Yuan Cao, Dongruo Zhou, and Quanquan Gu.  Gradient descent optimizes over-parameterized deep relu networks.  Machine Learning, 109(3):467–492, 2020.

##  Appendix A Appendix

###  A.1 Proof of Lemma 3.1

For each $\ell$, we have

$\displaystyle\left\|Z^{\ell+1}-Z^{\ell}\right\|_{F}=$
$\displaystyle\left\|\sigma(\gamma Z^{\ell}A+\Phi)-\sigma(\gamma
Z^{\ell-1}A+\Phi)\right\|_{F}$ $\displaystyle\leq$
$\displaystyle\gamma\|Z^{\ell}A-Z^{\ell-1}A\|_{F}\quad\text{by the Lipschitz
continuity of $\sigma$}$ $\displaystyle\leq$
$\displaystyle\gamma\|A\|\|Z^{\ell}-Z^{\ell-1}\|_{F}$ $\displaystyle<$
$\displaystyle\gamma_{0}\|Z^{\ell}-Z^{\ell-1}\|_{F}.$

Repeating the above argument $\ell$ times yields

$\displaystyle\left\|Z^{\ell+1}-Z^{\ell}\right\|_{F}\leq\gamma_{0}^{\ell}\left\|Z^{1}-Z^{0}\right\|_{F}=\gamma_{0}^{\ell}\left\|Z^{1}\right\|_{F}=\gamma_{0}^{\ell}\left\|\Phi\right\|_{F},$

where the second last equality is due to $Z^{0}=0$. Therefore, for all $p,q$
with $p\leq q$, we have

$\displaystyle\|Z^{p}-Z^{q}\|_{F}\leq$
$\displaystyle\|Z^{p}-Z^{p+1}\|_{F}+\cdots+\|Z^{q-1}-Z^{q}\|_{F}$
$\displaystyle\leq$
$\displaystyle\gamma_{0}^{p}\|\Phi\|_{F}+\cdots+\gamma_{0}^{q-1}\left\|\Phi\right\|_{F}$
$\displaystyle\leq$
$\displaystyle\gamma_{0}^{p}\|\Phi\|_{F}\left(1+\gamma_{0}+\gamma_{0}^{2}+\cdots\right)$
$\displaystyle=$ $\displaystyle\frac{\gamma_{0}^{p}}{1-\gamma_{0}}\|\Phi\|.$

Let $p,q\rightarrow\infty$, then $\|Z^{p}-Z^{q}\|_{F}\rightarrow 0$. Thus,
$\\{Z^{\ell}\\}_{\ell=0}^{\infty}$ is a Cauchy sequence in
${\mbox{$\mathbb{R}$}}^{N\times m}$. By the completeness of
${\mbox{$\mathbb{R}$}}^{N\times m}$, $\\{Z^{\ell}\\}_{\ell=0}^{\infty}$
converges to the unique limit $Z$, which is the equilibrium point. Now, let
$p=0$ and $q=\ell$, then we have

$\displaystyle\|Z^{\ell}\|_{F}\leq\frac{1}{1-\gamma_{0}}\left\|\Phi\right\|_{F}.$

Moreover, let $\ell\rightarrow\infty$, we obtain
$\|Z\|_{F}\leq\frac{1}{1-\gamma_{0}}\left\|\Phi\right\|_{F}$.

###  A.2 Proof of Lemma 3.2

Since the equilibrium point $Z$ is the root of the function $f$ defined by

$\displaystyle f(Z,A,W):=Z-\sigma(\gamma ZA+\Phi),$

then the differential of $f$ is given by

$\displaystyle df=$ $\displaystyle dZ-d\sigma(\gamma ZA+\Phi)$
$\displaystyle=$ $\displaystyle dZ-\sigma^{\prime}(U)\odot d(\gamma ZA+\Phi)$
$\displaystyle=$ $\displaystyle
dZ-\sigma^{\prime}(U)\odot\gamma(dZ)A-\sigma^{\prime}(U)\odot\gamma
ZdA-\sigma^{\prime}(U)\odot d\Phi,$

where $U:=\gamma ZA+\Phi$.

Taking vectorization on both sides yields

$\displaystyle\text{vec}\left(df\right)=$
$\displaystyle\text{vec}\left(dZ-\sigma^{\prime}(U)\odot\gamma(dZ)A-\sigma^{\prime}(U)\odot\gamma
ZdA-\sigma^{\prime}(U)\odot d\Phi\right)$ $\displaystyle=$
$\displaystyle\text{vec}\left(dZ\right)-\gamma
D\text{vec}\left(dZA\right)-\gamma
D\text{vec}\left(ZdA\right)-D\text{vec}\left(d\Phi\right)$ $\displaystyle=$
$\displaystyle\text{vec}\left(dZ\right)-\gamma D\left(A^{T}\otimes
I_{N}\right)\text{vec}\left(dZ\right)-\gamma D\left(I_{m}\otimes
Z\right)\text{vec}\left(dA\right)-D\text{vec}\left(d\Phi\right),$

where $D:=\mathop{\bf diag}[\text{vec}\left(\sigma^{\prime}(U)\right)]$. Thus,
we obtain

$\displaystyle\frac{\partial f}{\partial Z}=$
$\displaystyle\left[I_{Nm}-\gamma D(A^{T}\otimes I_{N})\right]^{T}$
$\displaystyle\frac{\partial f}{\partial A}=$
$\displaystyle-\gamma\left[D(I_{m}\otimes Z)\right]^{T}$
$\displaystyle\frac{\partial f}{\partial\Phi}=$ $\displaystyle-D^{T}.$

By reverse triangle inequality, we have

$\displaystyle\left\|I_{Nm}-\gamma D(A^{T}\otimes I_{N})\right\|\geq
1-\gamma\|D(A^{T}\otimes I_{N})\|\geq 1-\gamma\|A\|>1-\gamma_{0}>0,$ (29)

where we use $\|A\|\leq M$ and $\gamma_{0}=\gamma M<1$. Hence, the matrix
$Q:=\left[I_{Nm}-\gamma D(A^{T}\otimes I_{N})\right]$ is invertible.

Since the fixed point $Z$ is the root of $f$ and the matrix $Q$ is invertible,
by using the implicit function theorem, we have

$\displaystyle\frac{\partial Z}{\partial A}\frac{\partial f}{\partial
Z}+\frac{\partial f}{\partial A}=0,$

which implies

$\displaystyle\frac{\partial Z}{\partial A}=-\left(\frac{\partial f}{\partial
A}\right)\left(\frac{\partial f}{\partial
Z}\right)^{-1}=\gamma\left[D(I_{m}\otimes Z)\right]^{T}Q^{-T}.$

Similarly, we have

$\displaystyle\frac{\partial Z}{\partial\Phi}=-\left(\frac{\partial
f}{\partial\Phi}\right)\left(\frac{\partial f}{\partial
Z}\right)^{-1}=D^{T}Q^{-T}.$

Moreover, the definition of $\Phi$ implies that

$\displaystyle d\Phi=d\sigma(XW)=\sigma^{\prime}(XW)\odot XdW.$

Taking vectorization on both sides yields

$\displaystyle\text{vec}\left(d\Phi\right)=E\text{vec}\left(XdW\right)=E\left(I_{m}\otimes
X\right)\text{vec}\left(dW\right).$

where $E:=\mathop{\bf
diag}\left[\text{vec}\left(\sigma^{\prime}(XW)\right)\right]$. Therefore, we
have

$\displaystyle\frac{\partial\Phi}{\partial W}=\left[E\left(I_{m}\otimes
X\right)\right]^{T}.$

By using the chain rule, we have

$\displaystyle\frac{\partial Z}{\partial W}=\left(\frac{\partial\Phi}{\partial
W}\right)\left(\frac{\partial
Z}{\partial\Phi}\right)=\left[DE\left(I_{m}\otimes X\right)\right]^{T}Q^{-T}.$

Next, we have $\hat{y}=Zb$, so that

$\displaystyle d\hat{y}=dZb=(dZ)b+Zdb.$

Taking vectorization on both sides yields

$\displaystyle\text{vec}\left(d\hat{y}\right)=\left(b^{T}\otimes
I_{N}\right)\text{vec}\left(dZ\right)+Z\text{vec}\left(db\right),$

so that we have

$\displaystyle\frac{\partial\hat{y}}{\partial Z}=$
$\displaystyle\left(b^{T}\otimes I_{N}\right)^{T}$ (30)
$\displaystyle\frac{\partial\hat{y}}{\partial b}=$ $\displaystyle Z^{T}.$ (31)

By using the chain rule, we obtain

$\displaystyle\frac{\partial\hat{y}}{\partial A}=$
$\displaystyle\frac{\partial Z}{\partial A}\frac{\partial\hat{y}}{\partial
Z}=\gamma\left[D(I_{m}\otimes Z)\right]^{T}Q^{-T}\left(b^{T}\otimes
I_{N}\right)^{T}$ (32) $\displaystyle\frac{\partial\hat{y}}{\partial W}=$
$\displaystyle\frac{\partial Z}{\partial W}\frac{\partial\hat{y}}{\partial
Z}=\left[DE\left(I_{m}\otimes X\right)\right]^{T}Q^{-T}\left(b^{T}\otimes
I_{N}\right)^{T}.$ (33)

Note that

$\displaystyle dL=(\hat{y}-y)^{T}d\hat{y}$

so that

$\displaystyle\frac{\partial L}{\partial\hat{y}}=\hat{y}-y.$

Therefore, the partial derivative of $L$ with respect to $W$, $A$, and $b$ are
given by

$\displaystyle\frac{\partial L}{\partial W}=$
$\displaystyle\frac{\partial\hat{y}}{\partial W}\frac{\partial
L}{\partial\hat{y}}=\left[DE\left(I_{m}\otimes
X\right)\right]^{T}Q^{-T}\left(b^{T}\otimes I_{N}\right)^{T}(\hat{y}-y)$
$\displaystyle\frac{\partial L}{\partial A}=$
$\displaystyle\frac{\partial\hat{y}}{\partial A}\frac{\partial
L}{\partial\hat{y}}=\gamma\left[D(I_{m}\otimes
Z)\right]^{T}Q^{-T}\left(b^{T}\otimes I_{N}\right)^{T}(\hat{y}-y)$
$\displaystyle\frac{\partial L}{\partial b}=$
$\displaystyle\frac{\partial\hat{y}}{\partial b}\frac{\partial
L}{\partial\hat{y}}=Z^{T}(\hat{y}-y).$

###  A.3 Proof of Lemma 4.1

The gradient flows are given by

$\displaystyle\frac{d\text{vec}\left(W\right)}{dt}=-\frac{\partial L}{\partial
W},\quad\frac{d\text{vec}\left(A\right)}{dt}=-\frac{\partial L}{\partial
A},\quad\frac{db}{dt}=-\frac{\partial L}{\partial b}.$

By using the chain rule, we obtain the dynamics system of the equilibrium
point as follows

$\displaystyle\frac{d\text{vec}\left(Z\right)}{dt}=$
$\displaystyle\left(\frac{\partial Z}{\partial
A}\right)^{T}\left(\frac{d\text{vec}\left(A\right)}{dt}\right)+\left(\frac{\partial
Z}{\partial W}\right)^{T}\left(\frac{d\text{vec}\left(W\right)}{dt}\right)$
$\displaystyle=$ $\displaystyle-\left(\frac{\partial Z}{\partial
A}\right)^{T}\left(\frac{\partial L}{\partial A}\right)-\left(\frac{\partial
Z}{\partial W}\right)^{T}\left(\frac{\partial L}{\partial W}\right)$
$\displaystyle=$ $\displaystyle-\left[\left(\frac{\partial Z}{\partial
A}\right)^{T}\left(\frac{\partial Z}{\partial A}\right)+\left(\frac{\partial
Z}{\partial W}\right)^{T}\left(\frac{\partial Z}{\partial
W}\right)\right]\left(\frac{\partial\hat{y}}{\partial
Z}\right)\left(\frac{\partial L}{\partial\hat{y}}\right).$

Accordingly, the dynamics of the prediction $\hat{y}$ is given by

$\displaystyle\frac{d\hat{y}}{dt}=$
$\displaystyle\left(\frac{\partial\hat{y}}{\partial
b}\right)^{T}\left(\frac{db}{dt}\right)+\left(\frac{\partial\hat{y}}{\partial
Z}\right)^{T}\left(\frac{d\text{vec}\left(Z\right)}{dt}\right)$
$\displaystyle=$ $\displaystyle-\left\\{\left(\frac{\partial\hat{y}}{\partial
b}\right)^{T}\left(\frac{\partial\hat{y}}{\partial
b}\right)+\left(\frac{\partial\hat{y}}{\partial
A}\right)^{T}\left(\frac{\partial\hat{y}}{\partial
A}\right)+\left(\frac{\partial\hat{y}}{\partial
W}\right)^{T}\left(\frac{\partial\hat{y}}{\partial
W}\right)\right\\}\left(\frac{\partial L}{\partial\hat{y}}\right)$
$\displaystyle\triangleq$ $\displaystyle-H(t)(\hat{y}(t)-y),$

where

$\displaystyle\left(\frac{\partial\hat{y}}{\partial
b}\right)^{T}\left(\frac{\partial\hat{y}}{\partial b}\right)=Z(t)Z(t)^{T}$
$\displaystyle\left(\frac{\partial\hat{y}}{\partial
A}\right)^{T}\left(\frac{\partial\hat{y}}{\partial
A}\right)=\gamma^{2}\left\\{\left[D(I_{m}\otimes
Z)\right]^{T}Q^{-T}\left(b^{T}\otimes
I_{N}\right)^{T}\right\\}^{T}\left\\{\left[D(I_{m}\otimes
Z)\right]^{T}Q^{-T}\left(b^{T}\otimes I_{N}\right)^{T}\right\\}$
$\displaystyle\left(\frac{\partial\hat{y}}{\partial
W}\right)^{T}\left(\frac{\partial\hat{y}}{\partial
W}\right)=\left\\{\left[DE\left(I_{m}\otimes
X\right)\right]^{T}Q^{-T}\left(b^{T}\otimes
I_{N}\right)^{T}\right\\}^{T}\left\\{\left[DE\left(I_{m}\otimes
X\right)\right]^{T}Q^{-T}\left(b^{T}\otimes I_{N}\right)^{T}\right\\}.$

###  A.4 Proof of Theorem 4.2

###### Proof.

We make the inductive hypothesis as follows for all $0\leq s\leq t$

  1. (i)

$\|W(s)\|\leq\lambda_{1}$, $\|A(s)\|\leq\lambda_{2}$,
$\|b(s)\|\leq\lambda_{3}$,

  2. (ii)

$\sigma_{\min}(Z(s))\geq\alpha_{0}/2$,

  3. (iii)

$\|\hat{y}(s)-y\|^{2}\leq\exp\\{-(\alpha_{0}^{2}/2)s\\}\|\hat{y}(0)-y\|^{2}$.

For any $0\leq s\leq t$, we have

$\displaystyle\|\Phi(s)\|_{F}=\|\sigma(XW(s))\|_{F}\leq\|XW(s)\|_{F}\leq\|X\|_{F}\|W(s)\|\leq\lambda_{1}\|X\|_{F},$
(34)

where the last inequality follows from the inductive hypothesis (i). Then
Lemma 3.1 implies that

$\displaystyle\left\|Z(s)\right\|_{F}\leq\frac{1}{1-\gamma_{0}}\|\Phi(s)\|_{F}\leq\frac{\lambda_{1}}{1-\gamma_{0}}\left\|X\right\|_{F}.$
(35)

Then we can bound the partial derivative of $L$ with respect to $b$ as follows

$\displaystyle\|(\partial L/\partial b)(s)\|=$
$\displaystyle\|Z(s)^{T}(\hat{y}(s)-y)\|$ $\displaystyle\leq$
$\displaystyle\|Z(s)\|\|\hat{y}(s)-y\|$ $\displaystyle\leq$
$\displaystyle\frac{\lambda_{1}}{1-\gamma_{0}}\|X\|_{F}\|\hat{y}(s)-y\|,\quad\text{by
Eq.\eqref{eq: Z(s)}}$ (36) $\displaystyle\leq$
$\displaystyle\frac{\lambda_{1}}{1-\gamma_{0}}\|X\|_{F}\exp\left\\{-(\alpha_{0}^{2}/4)s\right\\}\|\hat{y}(0)-y\|,$
(37)

where the last inequality is due to the inductive hypothesis (iii). Then we
can bound the difference between $b(t)$ and $b(0)$ as follows

$\displaystyle\|b(t)-b(0)\|\leq\int_{0}^{t}\|(\partial L/\partial
b)(s)\|ds\leq\frac{4\lambda_{1}}{\alpha_{0}^{2}(1-\gamma_{0})}\|X\|_{F}\|\hat{y}(0)-y\|\leq
C_{3},$

where the last inequality follows from the initial assumptions (14)-(15). By
Weyl’s inequality, we obtain

$\displaystyle\|b(t)\|\leq\|b(t)-b(0)\|+\|b(0)\|\leq\lambda_{3}.$

With the similar argument, we can bound the partial derivative of $L$ with
respect to $W$

$\displaystyle\left\|(\partial L/\partial W)(s)\right\|=$
$\displaystyle\left\|\left[D(s)E(s)\left(I_{m}\otimes
X\right)\right]^{T}Q(s)^{-T}\left(b(s)^{T}\otimes
I_{N}\right)^{T}(\hat{y}(s)-y)\right\|$ $\displaystyle\leq$
$\displaystyle\|D(s)\|\|E(s)\|\|X\|\|Q(s)^{-1}\|\|b(s)\|\|\hat{y}(s)-y\|$
$\displaystyle\leq$
$\displaystyle\|X\|_{F}\cdot\frac{1}{1-\gamma_{0}}\cdot\lambda_{3}\cdot\|\hat{y}(s)-y\|$
$\displaystyle\leq$
$\displaystyle\frac{\lambda_{3}}{1-\gamma_{0}}\|X\|_{F}\exp\left\\{-(\alpha_{0}^{2}/4)s\right\\}\|\hat{y}(0)-y\|,$

and the difference between $W(t)$ and $W(0)$

$\displaystyle\|W(t)-W(0)\|\leq\|W(t)-W(0)\|_{F}\leq\int_{0}^{t}\left\|(\partial
L/\partial
W)(s)\right\|ds\leq\frac{4\lambda_{3}}{\alpha_{0}^{2}(1-\gamma_{0})}\|X\|_{F}\|\hat{y}(0)-y\|\leq
C_{1},$

which further implies

$\displaystyle\|W(t)\|\leq\|W(t)-W(0)\|+\|W(0)\|\leq\lambda_{1}.$

Similarly, we can bound the partial derivative of $L$ with respect to $A$ as
follows

$\displaystyle\left\|(\partial L/\partial A)(s)\right\|=$
$\displaystyle\left\|\gamma\left[D(s)(I_{m}\otimes
Z(s))\right]^{T}Q(s)^{-T}\left(b(s)^{T}\otimes
I_{N}\right)^{T}(\hat{y}(s)-y)\right\|$ $\displaystyle\leq$
$\displaystyle\gamma\|D(s)\|\|Z(s)\|\|Q(s)^{-1}\|\|b(s)\|\|\hat{y}(s)-y\|$
$\displaystyle\leq$
$\displaystyle\gamma\cdot\frac{\lambda_{1}}{1-\gamma_{0}}\|X\|_{F}\cdot\frac{1}{1-\gamma_{0}}\cdot\lambda_{3}\cdot\|\hat{y}(s)-y\|,\quad\text{By
$\eqref{eq: Z(s)}$ and inductive hypothesis (iii)}$ $\displaystyle=$
$\displaystyle\frac{\gamma_{0}}{(1-\gamma_{0})^{2}}\cdot\frac{\lambda_{1}\lambda_{3}}{\lambda_{2}}\cdot\|X\|_{F}\cdot\|\hat{y}(s)-y\|,\quad\text{By
}\gamma_{0}=\gamma\lambda_{2}$ $\displaystyle\leq$
$\displaystyle\frac{\gamma_{0}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}\lambda_{3}}{\lambda_{2}}\|X\|_{F}\exp\\{-(\alpha_{0}^{2}/4)s\\}\|\hat{y}(0)-y\|,$

and the difference between $A(t)$ and $A(0)$ as follows

$\displaystyle\|A(t)-A(0)\|\leq$ $\displaystyle\|A(t)-A(0)\|_{F}$
$\displaystyle\leq$ $\displaystyle\int_{0}^{t}\left\|(\partial L/\partial
A)(s)\right\|ds$ $\displaystyle\leq$
$\displaystyle\frac{4\gamma_{0}}{\alpha_{0}^{2}(1-\gamma_{0})^{2}}\frac{\lambda_{1}\lambda_{3}}{\lambda_{2}}\|X\|_{F}\|\hat{y}(0)-y\|$
$\displaystyle\leq$ $\displaystyle C_{2}.$

Therefore, we obtain

$\displaystyle\|A(t)\|\leq\|A(t)-A(0)\|+\|A(0)\|\leq\lambda_{2}.$

It follows from the inductive hypothesis and analysis above that
$\|A(s)\|\leq\lambda_{2}$ for all $0\leq s\leq t$. Since
$\gamma_{0}=\lambda_{2}\gamma<1$, it follows from Lemma 3.1 that the unique
equilibrium point $Z(s)$ always exists for all $0\leq s\leq t$. Then we can
derive the dynamics of the equilibrium point as follows

$\displaystyle\frac{d\text{vec}\left(Z\right)}{dt}=$
$\displaystyle\left(\frac{\partial Z}{\partial
A}\right)^{T}\left(\frac{d\text{vec}\left(A\right)}{dt}\right)+\left(\frac{\partial
Z}{\partial W}\right)^{T}\left(\frac{d\text{vec}\left(W\right)}{dt}\right)$
$\displaystyle=$ $\displaystyle\left(\frac{\partial Z}{\partial
A}\right)^{T}\left(-\frac{\partial L}{\partial A}\right)+\left(\frac{\partial
Z}{\partial W}\right)^{T}\left(-\frac{\partial L}{\partial W}\right)$
$\displaystyle=$ $\displaystyle-\left[\left(\frac{\partial Z}{\partial
A}\right)^{T}\left(\frac{\partial Z}{\partial A}\right)+\left(\frac{\partial
Z}{\partial W}\right)^{T}\left(\frac{\partial Z}{\partial
W}\right)\right]\left(\frac{\partial\hat{y}}{\partial
Z}\right)\left(\frac{\partial L}{\partial\hat{y}}\right),$

Note that for each $0\leq s\leq t$, we have

$\displaystyle\|(\partial Z/\partial A)(s)\|=$
$\displaystyle\left\|\gamma\left[D(s)(I_{m}\otimes
Z(s))\right]^{T}Q(s)^{-T}\right\|$ $\displaystyle\leq$
$\displaystyle\gamma\cdot\frac{\lambda_{1}}{(1-\gamma_{0})}\cdot\|X\|_{F}\cdot\frac{1}{(1-\gamma_{0})}$
$\displaystyle=$
$\displaystyle\frac{\gamma_{0}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}}{\lambda_{2}}\|X\|_{F},$

and

$\displaystyle\|(\partial Z/\partial W)(s)\|=$
$\displaystyle\left\|\left[D(s)E(s)\left(I_{m}\otimes
X\right)\right]^{T}Q(s)^{-T}\right\|$ $\displaystyle\leq$
$\displaystyle\|X\|_{F}\|Q(s)^{-1}\|$ $\displaystyle\leq$
$\displaystyle\frac{1}{(1-\gamma_{0})}\|X\|_{F}$

and

$\displaystyle\left\|\left(\frac{\partial\hat{y}}{\partial
Z}\right)\left(\frac{\partial
L}{\partial\hat{y}}\right)\right\|\leq\left\|b(s)\right\|\|\hat{y}(s)-y\|\leq\lambda_{3}\exp\left\\{-(\alpha_{0}^{2}/4)s\right\\}\|\hat{y}(0)-y\|.$

Therefore, we obtain

$\displaystyle\left\|\frac{d\text{vec}\left(Z\right)}{ds}\right\|\leq$
$\displaystyle\Big{(}\|(\partial Z/\partial A)(s)\|^{2}+\|(\partial Z/\partial
W)(s)\|^{2}\Big{)}\left\|\left(\frac{\partial\hat{y}}{\partial
Z}\right)\left(\frac{\partial L}{\partial\hat{y}}\right)\right\|$
$\displaystyle\leq$
$\displaystyle\left[\frac{\gamma_{0}^{2}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}^{2}}{\lambda_{2}^{2}}+1\right]\frac{\lambda_{3}}{(1-\gamma_{0})^{2}}\|X\|_{F}^{2}\exp\left\\{-(\alpha_{0}^{2}/4)s\right\\}\|\hat{y}(0)-y\|,$

and

$\displaystyle\left\|Z(t)-Z(0)\right\|\leq\left\|Z(t)-Z(0)\right\|_{F}\leq$
$\displaystyle\int_{0}^{t}\left\|\frac{d\text{vec}\left(Z\right)}{ds}\right\|ds$
$\displaystyle\leq$
$\displaystyle\left[\frac{\gamma_{0}^{2}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}^{2}}{\lambda_{2}^{2}}+1\right]\frac{4\lambda_{3}}{\alpha_{0}^{2}(1-\gamma_{0})^{2}}\|X\|^{2}\|\hat{y}(0)-y\|$
$\displaystyle\leq$ $\displaystyle\alpha_{0}/2,$

where the last inequality follows from the initial assumption.

Since $\sigma_{\min}{Z(0)}=\alpha_{0}$, Weyl’s inequality implies that

$\displaystyle\sigma_{\min}\left\\{Z(t)\right\\}\geq\sigma_{\min}\left\\{Z(0)\right\\}-\|Z(0)-Z(t)\|\geq\alpha_{0}/2.$

As a result, we have

$\displaystyle\lambda_{\min}\left\\{H(t)\right\\}\geq\lambda_{\min}\left\\{Z(t)Z(t)^{T}\right\\}=\sigma_{\min}\\{Z(t)\\}^{2}\geq\alpha_{0}^{2}/4,$

and so that

$\displaystyle\frac{d}{dt}L(\theta(t))=$
$\displaystyle\frac{d}{dt}\left(\frac{1}{2}\|\hat{y}(t)-y\|^{2}\right)$
$\displaystyle=$ $\displaystyle-(\hat{y}(t)-y)^{T}\frac{d\hat{y}(t)}{dt}$
$\displaystyle=$ $\displaystyle-(\hat{y}(t)-y)^{T}H(t)(\hat{y}(t)-y)$
$\displaystyle\leq$ $\displaystyle-(\alpha_{0}^{2}/4)\|\hat{y}(t)-y\|^{2}$
$\displaystyle=$ $\displaystyle-(\alpha_{0}^{2}/2)L(\theta(t)).$

Solving the above ordinary differential equation yields

$\displaystyle
L(\theta(t))\leq\exp\left\\{-(\alpha_{0}^{2}/2)t\right\\}L(\theta(0)).$

This completes the proof. ∎

###  A.5 Proof of Lemma 4.3

Since $\|A_{a}\|<\gamma^{-1}$ and $\|A_{b}\|<\gamma^{-1}$, Lemma 3.1 implies
the corresponding equilibrium points $Z_{a}$ and $Z_{b}$ are uniquely
determined. Note that $Z_{i}=\lim_{\ell\rightarrow\infty}Z_{i}^{\ell}$ for
$i\in\\{a,b\\}$. Then for any $\ell\geq 1$, we have

$\displaystyle\|Z_{a}^{\ell+1}-Z_{b}^{\ell+1}\|_{F}=$
$\displaystyle\|\sigma(\gamma Z_{a}^{\ell}A_{a}+\Phi_{a})-\sigma(\gamma
Z_{b}^{\ell}A_{b}+\Phi_{b})\|_{F}$ $\displaystyle\leq$ $\displaystyle\|\gamma
Z_{a}^{\ell}A_{a}+\Phi_{a}-\gamma Z_{b}^{\ell}A_{b}-\Phi_{b}\|_{F}$
$\displaystyle\leq$
$\displaystyle\gamma\|Z_{a}^{\ell}A_{a}-Z_{a}^{\ell}A_{b}\|_{F}+\gamma\|Z_{a}^{\ell}A_{b}-Z_{b}^{\ell}A_{b}\|_{F}+\|\Phi_{a}-\Phi_{b}\|_{F}$
$\displaystyle\leq$
$\displaystyle\gamma\|Z_{a}^{\ell}\|\|A_{a}-A_{b}\|_{F}+\gamma\|Z_{a}^{\ell}-Z_{b}^{\ell}\|_{F}\|A_{b}\|+\|\Phi_{a}-\Phi_{b}\|_{F}.$

In the rest of the proof, we will bound each term above. Lemma 3.1 implies
that

$\displaystyle\gamma\|Z_{a}^{\ell}\|\|A_{a}-A_{b}\|_{F}\leq\gamma\cdot\frac{1}{1-\gamma_{0}}\|\Phi_{a}\|_{F}\cdot\|A_{a}-A_{b}\|\leq\frac{\gamma_{0}}{1-\gamma_{0}}\frac{\lambda_{1}}{\lambda_{2}}\|X\|_{F}\|A_{a}-A_{b}\|.$

Since $\|A_{b}\|\leq\lambda_{2}$, the second term has the following inequality

$\displaystyle\gamma\|Z_{a}^{\ell}-Z_{b}^{\ell}\|\|A_{b}\|_{F}\leq\gamma_{0}\|Z_{a}^{\ell}-Z_{b}^{\ell}\|.$

By using the Lipschitz continuity of $\sigma(\cdot)$, we have

$\displaystyle\|\Phi_{a}-\Phi_{b}\|_{F}=\|\sigma(XW_{a})-\sigma(XW_{b})\|_{F}\leq\|X\|_{F}\|W_{a}-W_{b}\|.$

Combining all results together yields

$\displaystyle\|Z_{a}^{\ell+1}-Z_{b}^{\ell+1}\|_{F}\leq\gamma_{0}\|Z_{a}^{\ell}-Z_{b}^{\ell}\|_{F}+\left(\frac{\gamma_{0}}{1-\gamma_{0}}\frac{\lambda_{1}}{\lambda_{2}}\|X\|_{F}\|A_{a}-A_{b}\|+\|X\|_{F}\|W_{a}-W_{b}\|\right)$

Let
$C:=\frac{\gamma_{0}}{1-\gamma_{0}}\frac{\lambda_{1}}{\lambda_{2}}\|X\|_{F}\|A_{a}-A_{b}\|+\|X\|_{F}\|W_{a}-W_{b}\|$
and apply the same argument $\ell$ times yields

$\displaystyle\|Z_{a}^{\ell+1}-Z_{b}^{\ell+1}\|_{F}\leq\gamma_{0}^{\ell}\|Z_{a}^{1}-Z_{b}^{1}\|_{F}+\left(1+\gamma_{0}+\cdots+\gamma_{0}^{\ell-1}\right)C.$

Let $\ell\rightarrow\infty$. The continuity of the operator norm implies that

$\displaystyle\|Z_{a}-Z_{b}\|_{F}\leq\frac{1}{1-\gamma_{0}}C.$

###  A.6 Proof of Theorem 4.4

###### Proof.

We make the inductive hypothesis as follows for all $0\leq s\leq k$

  1. (i)

$\|W(s)\|\leq\lambda_{1}$, $\|A(s)\|\leq\lambda_{2}$,
$\|b(s)\|\leq\lambda_{3}$,

  2. (ii)

$\sigma_{\min}(Z(s))\geq\alpha_{0}/2$,

  3. (iii)

$L(\theta(s))\leq(1-\eta\alpha_{0}^{2}/4)^{s}L(\theta(0))$.

For any $0\leq s\leq k$, the inductive hypothesis implies that

$\displaystyle\|\Phi(s)\|=\|\sigma(XW(s))\|\leq\lambda_{1}\|X\|.$ (38)

By using Eq. (35), the partial derivative of $L$ with respect to $b$ can be
bounded as follows

$\displaystyle\|(\partial L/\partial b)(s)\|=$
$\displaystyle\|Z(s)^{T}(\hat{y}(s)-y)\|$ $\displaystyle\leq$
$\displaystyle\|Z(s)\|\|\hat{y}(s)-y\|$ $\displaystyle\leq$
$\displaystyle\frac{\lambda_{1}}{1-\gamma_{0}}\|X\|_{F}\|\hat{y}(s)-y\|,$ (39)
$\displaystyle\leq$
$\displaystyle\frac{\lambda_{1}}{1-\gamma_{0}}\|X\|_{F}\|\hat{y}(0)-y\|\beta^{s},$
(40)

where $\beta:=\sqrt{1-\eta\alpha_{0}^{2}/4}$, and the second inequality is due
to Eq.(38), and the last inequality is by inductive hypothesis (iii). Then the
difference between $b(k+1)$ and $b(0)$ is given by

$\displaystyle\left\|b(k+1)-b(0)\right\|\leq$
$\displaystyle\eta\sum_{s=0}^{k}\left\|(\partial L/\partial b)(s)\right\|$
$\displaystyle\leq$
$\displaystyle\eta\frac{\lambda_{1}}{1-\gamma_{0}}\|X\|_{F}\|\hat{y}(0)-y\|\cdot\sum_{s=0}^{k}\beta^{s}$
$\displaystyle\leq$
$\displaystyle\frac{\eta}{1-\beta}\frac{\lambda_{1}}{1-\gamma_{0}}\|X\|_{F}\|\hat{y}(0)-y\|$
$\displaystyle\leq$
$\displaystyle\frac{8}{\alpha_{0}^{2}}\frac{\lambda_{1}}{1-\gamma_{0}}\|X\|_{F}\|\hat{y}(0)-y\|$
$\displaystyle\leq$ $\displaystyle C_{3},$

where the last inequality is due to the initial assumption. Then we obtain

$\displaystyle\|b(k+1)\|\leq\|b(k+1)-b(0)\|+\|b(0)\|\leq\lambda_{3}.$

With the similar argument, we can bound the partial derivative of $L$ with
respect to $W$

$\displaystyle\left\|(\partial L/\partial W)(s)\right\|=$
$\displaystyle\left\|\left[D(s)E(s)\left(I_{m}\otimes
X\right)\right]^{T}Q(s)^{-T}\left(b(s)^{T}\otimes
I_{N}\right)^{T}(\hat{y}(s)-y)\right\|$ $\displaystyle\leq$
$\displaystyle\|X\|_{F}\cdot\|Q(s)^{-1}\|\cdot\|b(s)\|\cdot\|\hat{y}(s)-y\|$
$\displaystyle\leq$
$\displaystyle\frac{\lambda_{3}}{1-\gamma_{0}}\|X\|_{F}\|\hat{y}(s)-y\|$
$\displaystyle\leq$
$\displaystyle\frac{\lambda_{3}}{1-\gamma_{0}}\|X\|_{F}\|\hat{y}(0)-y\|\cdot\beta^{s},$
(41)

and the difference between $W(k+1)$ and $W(0)$ is given by

$\displaystyle\|W(k+1)-W(0)\|\leq$ $\displaystyle\eta\sum_{s=0}^{k}\|(\partial
L/\partial W)(s)\|$ $\displaystyle\leq$
$\displaystyle\eta\frac{\lambda_{3}}{1-\gamma_{0}}\|X\|_{F}\|\hat{y}(0)-y\|\cdot\sum_{s=0}^{k}\beta^{s}$
$\displaystyle\leq$
$\displaystyle\frac{\eta}{1-\beta}\frac{\lambda_{3}}{1-\gamma_{0}}\|X\|_{F}\|\hat{y}(0)-y\|$
$\displaystyle\leq$
$\displaystyle\frac{8}{\alpha_{0}^{2}}\frac{\lambda_{3}}{1-\gamma_{0}}\|X\|_{F}\|\hat{y}(0)-y\|$
$\displaystyle\leq$ $\displaystyle C_{1},$

where the last inequality is due to the initial assumption (17). This implies

$\displaystyle\|W(k+1)\|\leq\|W(k)-W(0)\|+\|W(0)\|\leq\lambda_{1}.$

Similarly, we can bound the partial derivative of $L$ with respect to $A$ as
follows

$\displaystyle\left\|(\partial L/\partial A)(s)\right\|=$
$\displaystyle\left\|\gamma\left[D(s)(I_{m}\otimes
Z(s))\right]^{T}Q(s)^{-T}\left(b(s)^{T}\otimes
I_{N}\right)^{T}(\hat{y}(s)-y)\right\|$ $\displaystyle\leq$
$\displaystyle\gamma\|Z(s)\|\cdot\|Q(s)^{-1}\|\cdot\|b(s)\|\cdot\|\hat{y}(s)-y\|$
$\displaystyle\leq$
$\displaystyle\frac{\gamma}{(1-\gamma_{0})}\lambda_{1}\|X\|_{F}\cdot\frac{1}{1-\gamma_{0}}\cdot\lambda_{3}\cdot\|\hat{y}(s)-y\|$
$\displaystyle\leq$
$\displaystyle\frac{\gamma_{0}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}\lambda_{3}}{\lambda_{2}}\|X\|_{F}\cdot\|\hat{y}(s)-y\|$
$\displaystyle\leq$
$\displaystyle\frac{\gamma_{0}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}\lambda_{3}}{\lambda_{2}}\|X\|_{F}\cdot\|\hat{y}(0)-y\|\cdot\beta^{s},$
(42)

and the difference between $A(k+1)$ and $A(0)$ is given by

$\displaystyle\|A(k+1)-A(0)\|\leq$
$\displaystyle\eta\sum_{s=0}^{k}\left\|(\partial L/\partial A)(s)\right\|$
$\displaystyle\leq$
$\displaystyle\eta\frac{\gamma_{0}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}\lambda_{3}}{\lambda_{2}}\|X\|\cdot\|\hat{y}(0)-y\|\sum_{s=0}^{k}\beta^{s}$
$\displaystyle\leq$
$\displaystyle\frac{\eta}{1-\beta}\frac{\gamma_{0}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}\lambda_{3}}{\lambda_{2}}\|X\|\cdot\|\hat{y}(0)-y\|$
$\displaystyle\leq$
$\displaystyle\frac{8}{\alpha_{0}^{2}}\frac{\gamma_{0}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}\lambda_{3}}{\lambda_{2}}\|X\|\cdot\|\hat{y}(0)-y\|$
$\displaystyle\leq$ $\displaystyle C_{2},$

where the last inequality is due to the initial assumption (17). Therefore, we
obtain

$\displaystyle\|A(k+1)\|\leq\|A(k+1)-A(0)\|+\|A(0)\|\leq\lambda_{2}.$

Since $\|A(k+1)\|\leq\lambda_{2}$ and $\gamma_{0}=\gamma\lambda_{2}<1$, Lemma
3.1 implies the unique equilibrium point $Z(k+1)$ exists.

It follows from Lemma 4.3, Eq.(41), and Eq.(42) that for any $s$ and $\ell$,
we have

$\displaystyle\|Z(s+1)-Z(s)\|_{F}\leq\left[1+\frac{\gamma_{0}^{2}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}^{2}}{\lambda_{2}^{2}}\right]\eta\frac{\lambda_{3}}{(1-\gamma_{0})^{2}}\|X\|_{F}^{2}\|\hat{y}(s)-y\|.$
(43)

Therefore, we have

$\displaystyle\|Z(k)-Z(0)\|_{F}\leq$
$\displaystyle\sum_{s=0}^{k}\|Z(s+1)-Z(s)\|_{F}$ $\displaystyle\leq$
$\displaystyle\sum_{s=0}^{k}\left[1+\frac{\gamma_{0}^{2}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}^{2}}{\lambda_{2}^{2}}\right]\eta\frac{\lambda_{3}}{(1-\gamma_{0})^{2}}\|X\|^{2}\|\hat{y}(0)-y\|\cdot\beta^{s}$
$\displaystyle\leq$
$\displaystyle\left[1+\frac{\gamma_{0}^{2}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}^{2}}{\lambda_{2}^{2}}\right]\eta\frac{\lambda_{3}}{(1-\gamma_{0})^{2}}\|X\|_{F}^{2}\|\hat{y}(0)-y\|\cdot\frac{1}{1-\beta}$
$\displaystyle\leq$
$\displaystyle\left[1+\frac{\gamma_{0}^{2}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}^{2}}{\lambda_{2}^{2}}\right]\frac{\lambda_{3}}{(1-\gamma_{0})^{2}}\|X\|_{F}^{2}\|\hat{y}(0)-y\|\cdot\frac{8}{\alpha_{0}^{2}}$
$\displaystyle\leq$ $\displaystyle\alpha_{0}/2,$

where the last inequality is due to the initial assumption Eq.(17)-Eq.(19). As
a result, we obtain

$\displaystyle\sigma_{\min}\left[Z(k)\right]\geq\sigma_{\min}\left[Z(0)\right]-\|Z(k)-Z(0)\|\geq\alpha_{0}/2.$
(44)

Now, we are ready to derive the linear convergence result. Note that

$\displaystyle L(\theta(k+1))=$
$\displaystyle\frac{1}{2}\left\|\hat{y}(k+1)-y\right\|^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\|\hat{y}(k+1)-\hat{y}(k)+\hat{y}(k)-y\|^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\|\hat{y}(k+1)-\hat{y}(k)\|^{2}+\left\langle\hat{y}(k+1)-\hat{y}(k),\hat{y}(k)-y\right\rangle+L(\theta(k)).$

In the rest of the proof, we will bound each term above. Note that

$\displaystyle\|\hat{y}(k+1)-\hat{y}(k)\|=$
$\displaystyle\|Z(k+1)b(k+1)-Z(k)b(k)\|$ $\displaystyle\leq$
$\displaystyle\|Z(k+1)b(k+1)-Z(k)b(k+1)\|+\|Z(k)b(k+1)-Z(k)b(k)\|$
$\displaystyle\leq$
$\displaystyle\|Z(k+1)-Z(k)\|\|b(k+1)\|+\|Z(k)\|\|b(k+1)-b(k)\|$
$\displaystyle\leq$
$\displaystyle\left[1+\frac{\gamma_{0}^{2}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}^{2}}{\lambda_{2}^{2}}\right]\eta\frac{\lambda_{3}}{(1-\gamma_{0})^{2}}\|X\|_{F}^{2}\|\hat{y}(k)-y\|\cdot\lambda_{3}+\frac{1}{1-\gamma_{0}}\left\|\Phi(k)\right\|_{F}\cdot\eta\|(\partial
L/\partial b)(k)\|$ $\displaystyle\leq$
$\displaystyle\left[1+\frac{\gamma_{0}^{2}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}^{2}}{\lambda_{2}^{2}}\right]\eta\frac{\lambda_{3}^{2}}{(1-\gamma_{0})^{2}}\|X\|_{F}^{2}\|\hat{y}(k)-y\|+\eta\frac{\lambda_{1}^{2}}{(1-\gamma_{0})^{2}}\|X\|_{F}^{2}\|\hat{y}(k)-y\|$
$\displaystyle\leq$
$\displaystyle\eta\underbrace{\left(\left[1+\frac{\gamma_{0}^{2}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}^{2}}{\lambda_{2}^{2}}\right]\lambda_{1}^{-2}+\lambda_{3}^{-2}\right)\frac{1}{(1-\gamma_{0})^{2}}\lambda_{1}^{2}\lambda_{3}^{2}\|X\|_{F}^{2}}_{:=Q_{1}}\cdot\|\hat{y}(k)-y\|.$

Let $g:=Z(k)b(k+1)$, then we have

$\displaystyle\left\langle\hat{y}(k+1)-\hat{y}(k),\hat{y}(k)-y\right\rangle=\left\langle\hat{y}(k+1)-g,\hat{y}(k)-y\right\rangle+\left\langle
g-\hat{y}(k),\hat{y}(k)-y\right\rangle,$

where

$\displaystyle\left\langle\hat{y}(k+1)-g,\hat{y}(k)-y\right\rangle=$
$\displaystyle\left\langle Z(k+1)b(k+1)-Z(k)b(k+1),\hat{y}(k)-y\right\rangle$
$\displaystyle\leq$ $\displaystyle\|Z(k+1)-Z(k)\|\|b(k+1)\|\|\hat{y}(k)-y\|$
$\displaystyle\leq$
$\displaystyle\left[1+\frac{\gamma_{0}^{2}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}^{2}}{\lambda_{2}^{2}}\right]\eta\frac{\lambda_{3}}{(1-\gamma_{0})^{2}}\|X\|_{F}^{2}\|\hat{y}(k)-y\|\cdot\lambda_{3}\cdot\|\hat{y}(k)-y\|$
$\displaystyle=$
$\displaystyle\eta\underbrace{\left[1+\frac{\gamma_{0}^{2}}{(1-\gamma_{0})^{2}}\frac{\lambda_{1}^{2}}{\lambda_{2}^{2}}\right]\frac{\lambda_{3}^{2}}{(1-\gamma_{0})^{2}}\|X\|_{F}^{2}}_{:=Q_{2}}\cdot\|\hat{y}(k)-y\|^{2}$

and

$\displaystyle\left\langle g-\hat{y}(k),\hat{y}(k)-y\right\rangle=$
$\displaystyle\left\langle Z(k)b(k+1)-Z(k)b(k),\hat{y}(k)-y\right\rangle$
$\displaystyle=$
$\displaystyle-\eta(\hat{y}(k)-y)^{T}\left[Z(k)Z(k)^{T}\right](\hat{y}(k)-y)$
$\displaystyle\leq$
$\displaystyle-\eta\cdot(\alpha_{0}^{2}/4)\|\hat{y}(k)-y\|^{2}.$

Combining all the results together, we have

$\displaystyle L(\theta(k+1))\leq$
$\displaystyle\left[1-\eta\frac{\alpha_{0}^{2}}{2}+\eta^{2}Q_{1}^{2}+2\eta
Q_{2}\right]L(\theta(k))$ $\displaystyle\leq$
$\displaystyle\left[1-\eta\left(\frac{\alpha_{0}^{2}}{2}-4Q_{2}\right)\right]L(\theta(k)),\quad\text{By
Eq.\eqref{eq:step size}}$ $\displaystyle\leq$
$\displaystyle\left[1-\eta\frac{\alpha_{0}^{2}}{4}\right]L(\theta(k)),$

where the last inequality follows from Eq.(19). ∎

###  A.7 Proof of Lemma 5.1

Note that $\Phi_{ij}\geq 0$ by the definition of Eq.(3). Since $Z^{0}=0$, we
have

$\displaystyle Z^{1}=\sigma(\gamma Z^{0}A+\Phi)=\sigma(\Phi)=\Phi.$

Since $A_{ij}\geq 0$, we have

$\displaystyle Z^{2}=\sigma(\gamma Z^{1}A+\Phi)=\sigma(\gamma\Phi
A+\Phi)=\gamma\Phi A+\Phi=\Phi\left(\gamma A+I_{m}\right)$

Repeating the same argument $\ell$ times, we have

$\displaystyle
Z^{\ell+1}=\Phi\left(\hat{A}^{\ell}+\hat{A}^{\ell-1}+\cdots+\hat{A}+I_{m}\right),$

where $\hat{A}:=\gamma A$. Let $\ell\rightarrow\infty$, then Neumann series
implies that

$\displaystyle Z=\Phi(I_{m}-\hat{A})^{-1}.$

##  Appendix B More experimental results

In this section, we provide more experimental results to justify our
theoretical findings.

Figure 4: We evaluate the impact of the width $m$ on the training loss, test
loss, and operator norm of the scaled matrix $\gamma A(k)$ on the modified
dataset of MNIST, FashionMNIST, CIFAR10, and SVHN

Figure 5: We evaluate the impact of the learning rate $\eta$ on the training
loss, test loss, and operator norm of the scaled matrix $\gamma A(k)$ on the
modified dataset of MNIST, FashionMNIST, CIFAR10, and SVHN

Figure 6: We evaluate the impact of choice $\gamma$ on the training loss, test
loss, and operator norm of the scaled matrix $\gamma A(k)$ on the modified
dataset of MNIST, FashionMNIST, CIFAR10, and SVHN