A GLOBAL CONVERGENCE THEORY FOR DEEP RELU IMPLICIT NETWORKS VIA OVER-PARAMETERIZATION
Year: 2022

Abstract:
Implicit deep learning has received increasing attention recently, since it generalizes the recursive prediction rules of many commonly used neural network architectures. Its prediction rule is provided implicitly based on the solution of an equilibrium equation. Although many recent studies have experimentally demonstrates its superior performances, the theoretical understanding of implicit neural networks is limited. In general, the equilibrium equation may not be well-posed during the training. As a result, there is no guarantee that a vanilla (stochastic) gradient descent (SGD) training nonlinear implicit neural networks can converge. This paper fills the gap by analyzing the gradient flow of Rectified Linear Unit (ReLU) activated implicit neural networks. For an m-width implicit neural network with ReLU activation and n training samples, we show that a randomly initialized gradient descent converges to a global minimum at a linear rate for the square loss function if the implicit neural network is over-parameterized. It is worth noting that, unlike existing works on the convergence of (S)GD on finitelayer over-parameterized neural networks, our convergence results hold for implicit neural networks, where the number of layers is infinite.

heading: 1 INTRODUCTION
text: 1) Background and Motivation: In the last decade, implicit deep learning (El Ghaoui et al., 2019) have attracted more and more attention. Its popularity is mainly because it generalizes the recursive rules of many widely used neural network architectures. A line of recent works (Bai et al., 2019; El Ghaoui et al., 2019; Bai et al., 2020) have shown that the implicit neural network architecture is a wider class that includes most current neural network architectures as special cases, such as feed-forward neural networks, convolution neural networks, residual networks, and recurrent neural networks. Moreover, implicit deep learning is also well known for its competitive performance compared to other regular deep neural networks but using significantly fewer computational resources (Dabre & Fujita, 2019; Dehghani et al., 2018; Bai et al., 2018).
Although a line of literature has been shown the superior performance of implicit neural networks experimentally, the theoretical understanding is still limited. To date, it is still unknown if a simple first-order optimization method such as (stochastic) gradient descent can converge on an implicit neural network activated by a nonlinear function. Unlike a regular deep neural network, an implicit neural network could have infinitely many layers, resulting in the possibility of divergence of the
forward propagation (El Ghaoui et al., 2019; Kawaguchi, 2021). The main challenge in establishing the convergence of implicit neural network training lies in the fact that, in general, the equilibrium equation of implicit neural networks cannot be solved in closed-form. What exacerbates the problem is the well-posedness of the forward propagation. In other words, the equilibrium equation may have zero or multiple solutions. A line of recent studies have suggested a number of strategies to handle this well-posedness challenge, but they all involved reformulation or solving subproblems in each iteration. For example, El Ghaoui et al. (2019) suggested to reformulate the training as the Fenchel divergence formulation and solve the reformulated optimization problem by projected gradient descent method. However, this requires solving a projection subproblem in each iteration and convergence was only demonstrated numerically. By using an extra softmax layer, Kawaguchi (2021) established global convergence result of gradient descent for a linear implicit neural network. Unfortunately, their result cannot be extended to nonlinear activations, which are critical to the learnability of deep neural networks.
This paper proposes a global convergence theory of gradient descent for implicit neural networks activated by nonlinear Rectified Linear Unit (ReLU) activation function by using overparameterization. Specifically, we show that random initialized gradient descent with fixed stepsize converges to a global minimum of a ReLU implicit neural network at a linear rate as long as the implicit neural network is overparameterized. Recently, over-parameterization has been shown to be effective in optimizing finite-depth neural networks (Zou et al., 2020; Nguyen & Mondelli, 2020; Arora et al., 2019; Oymak & Soltanolkotabi, 2020). Although the objective function in the training is nonsmooth and non-convex, it can be shown that GD or SGD converge to a global minimum linearly if the width m of each layer is a polynomial of the number of training sample n and the number of layers h, i.e., m = poly(n, h). However, these results cannot be directly applied to implicit neural networks, since implicit neural networks have infinitely many hidden layers, i.e., h → ∞, and the well-posedness problem surfaces during the training process. In fact, Chen et al. (2018); Bai et al. (2019; 2020) have all observed that the time and number of iterations spent on forward propagation are gradually increased with the the training epochs. Thus, we have to ensure the unique equilibrium point always exists throughout the training given that the width m is only polynomial of n.
2) Preliminaries of Implicit Deep Learning: In this work, we consider an implicit neural network with the transition at the `-th layer in the following form (El Ghaoui et al., 2019; Bai et al., 2019):
z` = σ ( γ√ m Az`−1 + φ(x) ) , (1)
where φ : Rd → Rm is a feature mapping function that transforms an input vector x ∈ Rd to a desired feature vector φ , φ(x), z` ∈ Rm is the output of the `-th layer, A ∈ Rm×m is a trainable weight matrix, σ(u) = max{0, u} is the ReLU activation function, and γ ∈ (0, 1) is a fixed scalar to scaleA. As will be shown later in Section 2.1, γ plays the role of ensuring the existence of the limit z∗ = lim`→∞ z
`. In general, the feature mapping function φ is a nonlinear function, which extracts features from the low-dimensional input vector x. In this paper, we consider a simple nonlinear feature mapping function φ given by
φ(x) , 1√ m σ(Wx), (2)
where W ∈ Rm×d is a trainable parameter matrix. As ` → ∞, an implicit neural network can be considered an infinitely deep neural network. Consequently, z∗ is not only the limit of the sequence {z`}∞`=0 with z0 = 0, but it is also the equilibrium point (or fixed point) of the equilibrium equation:
z∗ = σ(γ̃Az∗ + φ), (3)
where γ̃ , γ/ √ m. In implicit neural networks, the prediction ŷ for the input vector x is the combination of the fixed point z∗ and the feature vector φ, i.e.,
ŷ = uTz∗ + vTφ, (4)
where u,v ∈ Rm are trainable weight vectors. For simplicity, we use θ , vec (A,W ,u,v) to group all training parameters. Given a training data set {(xi, yi)}ni=1, we want to minimize
L(θ) = n∑ i=1 1 2 (ŷi − yi)2 = 1 2 ‖ŷ − y‖2, (5)
where ŷ and y are the vectors formed by stacking all the prediction and labels.
3) Main Results: Our results are based on the following observations. We first analyze the forward propagation and find that the unique equilibrium point always exists if the scaled matrix γ̃A in Eq. (3) has an operator norm less than one. Thus, the well-posedness problem is reduced to finding a sequence of scalars {γk}∞k=1 such that γ̃kA(k) is appropriately scaled. To achieve this goal, we show that the operator norm A(k) is uniformly upper bounded by a constant over all iterations. Consequently, a fixed scalar γ is enough to ensure the well-posedness of Eq. (3). Our second observation is from the analysis of the gradient descent method with infinitesimal step-size (gradient flow). By applying the chain rule with the gradient flow, we derive the dynamics of prediction ŷ(t) which is governed by the spectral property of a Gram matrix. In particular, if the smallest eigenvalue of the Gram matrix is lower bounded throughout the training, the gradient descent method enjoys a linear convergence rate. Along with some basic functional analysis results, it can be shown that the smallest eigenvalue of the Gram matrix at initialization is lower bounded if no two data samples are parallel. Although the Gram matrix varies in each iteration, the spectral property is preserved if the Gram matrix is close to its initialization. Thus, the convergence problem is reduced to showing the Gram matrix in latter iterations is close to its initialization. Our third observation is that we find random initialization, over-parameterization, and linear convergence jointly enforce the (operator) norms of parameters upper bounded by some constants and close to their initialization. Accordingly, we can use this property to show that the operator norm ofA is upper bounded and the spectral property of the Gram matrix is preserved throughout the training. Combining all these insights together, we can conclude that the random initialized gradient descent method with a constant step-size converges to a global minimum of the implicit neural network with ReLU activation.
The main contributions of this paper are summarized as follows:
(i) By scaling the weight matrix A with a fixed scalar γ, we show that the unique equilibrium point z∗ for each x always exists during the training if the parameters are randomly initialized, even for the nonlinear ReLU activation function.
(ii) We analyze the gradient flow of implicit neural networks. Despite the non-smooth and nonconvexity of the objective function, the convergence to a global minimum at a linear rate is guaranteed if the implicit neural network is over-parameterized and the data is non-degenerate.
(iii) Since gradient descent is discretized version of gradient flow, we can show gradient descent with fixed stepsize converges to a global minimum of implicit neural networks at a linear rate under the same assumptions made by the gradient flow analysis, as long as the stepsize is chosen small enough.
Notation: For a vector x, ‖x‖ is the Euclidean norm of x. For a matrix A, ‖A‖ is the operator norm of A. If A is a square matrix, then λmin(A) and λmax(A) denote the smallest and largest eigenvalue ofA, respectively, and λmax(A) ≤ ‖A‖. We denote [n] , {1, 2, · · · , n}.

heading: 2 WELL-POSEDNESS AND GRADIENT COMPUTATION
text: In this section, we provide a simple condition for the equilibrium equation (3) to be well-posed in the sense that the unique equilibrium point exists. Instead of backpropagating through all the intermediate iterations of a forward pass, we derive the gradients of trainable parameters by using the implicit function theorem. In this work, we make the following assumption on parameter initialization. Assumption 1 (Random Initialization). The entries Aij and Wij are randomly initialized by the standard Gaussian distribution N (0, 1), and ui and vi are randomly initialized by the symmetric Bernoulli or Rademacher distribution. Remark 2.1. This initialization is similar to the approaches widely used in practice (Glorot & Bengio, 2010; He et al., 2015). The result obtained in this work can be easily extended to the case where the distributions forAij ,Wij , ui, and vi are replaced by sub-Gaussian random variables.

heading: 2.1 FORWARD PROPAGATION AND WELL-POSEDNESS
text: In a general implicit neural network, Eq. (3) is not necessarily well-posed, since it may admit zero or multiple solutions. In this work, we show that scaling the matrix A with γ̃ = γ/ √ m guarantees
the existence and uniqueness of the equilibrium point z∗ with random initialization. This follows from a foundational result in random matrix theory as restated in the following lemma. Lemma 2.1 (Vershynin (2018), Theorem 4.4.5). Let A be an m × n random matrix whose entries Aij are independent, zero-mean, and sub-Gaussian random variables. Then, for any t > 0, we have ‖A‖ ≤ CK( √ m + √ n + t) with probability at least 1 − 2e−t2 . Here C > 0 is a fixed constant, and K = maxi,j ‖Aij‖ψ2 . Under Assumption 1, Lemma 2.1 implies that, with exponentially high probability, ‖A‖ ≤ c √ m for some constant c > 0. By scalingA by a positive scalar γ̃, we show that the transition Eq. (1) is a contraction mapping. Thus, the unique equilibrium point exists with detailed proof in Appendix A.1. Lemma 2.2. If ‖A‖ ≤ c √ m for some c > 0, then for any γ0 ∈ (0, 1), the scalar γ , min{γ0, γ0/c} uniquely determines the existence of the equilibrium z∗ for every x, and ‖z`‖ ≤ 11−γ0 ‖φ‖ for all `.
Lemma 2.2 indicates that equilibria always exist if we can maintain the operator norm of the scaled matrix (γ/ √ m)A less than 1 during the training. However, the operator norms of matrix A are changed by the update of the gradient descent. It is hard to use a fixed scalar γ to scale the matrix A over all iterations, unless the operator norm of A is bounded. In Section 3, we will show ‖A‖ ≤ 2c √ m always holds throughout the training, provided that ‖A(0)‖ ≤ c √ m at initialization and the width m is sufficiently large. Thus, by using the scalar γ = min{γ0, γ0/(2c)} for any γ0 ∈ (0, 1), equilibria always exist and the equilibrium equation Eq. (3) is well-posed.

heading: 2.2 BACKWARD GRADIENT COMPUTING
text: For a regular network with finite layers, one needs to store all intermediate parameters and apply backpropagation to compute the gradients. In contrast, for an implicit neural network, we can derive the formula of the gradients by using the implicit function theorem. Here, the equilibrium point z∗ is a root of the function f given by f(z,A,W ) , z − σ (γ̃Az + φ). The essential challenge is to ensure the partial derivative ∂f/∂z at z∗ is always invertible throughout the training. The following lemma shows that the partial derivative ∂f/∂z at z∗ always exists with the scalar γ, and the gradient derived in the lemma always exists. The gradient derivations for each trainable parameters are provided in the following lemma and the proof is provided in Appendix A.2. Lemma 2.3 (Gradients of an Implicit Neural Network). Define the scalar γ , min{γ0, γ0/c}. If ‖A‖ ≤ c √ m for some constant c > 0, then for any γ0 ∈ (0, 1), the following results hold:
(i) The partial derivatives of f with respect to z,A, andW are ∂f
∂z = [Im − γ̃ diag(σ′(γ̃Az + φ))A] T , (6)
∂f ∂A = −γ̃
[ zT ⊗ diag(σ′(γ̃Az + φ)) ]T , (7)
∂f ∂W = − 1√ m
[ xT ⊗ diag(σ′(Wx)) ]T diag(σ′(γ̃Az + φ))T , (8)
where γ̃ , γ/ √ m.
(ii) For any vector v, the following inequality holds λmin {Im − γ̃ diag(σ′(v))A} > 1− γ0 > 0, (9)
which further implies ∂f/∂z is invertible at the equilibrium point z∗.
(iii) The gradient of the objective function L with respect toA andW are given by ∇uL = n∑ i=1 (ŷi − yi)zi, ∇vL = n∑ i=1 (ŷi − yi)φi, (10)
∇AL = n∑ i=1 γ̃(ŷi − yi)DTi [Im − γ̃DiA] −T uzTi , (11)
∇WL = n∑ i=1 1√ m (ŷi − yi)ETi { DTi [Im − γ̃DiA] −T u+ v } xTi , (12)
where zi is the equilibrium point for the training data (xi, yi), Di , diag(σ′(γ̃Azi + φi)), and Ei , diag(σ′(Wxi)).

heading: 3 GLOBAL CONVERGENCE OF THE GRADIENT DESCENT METHOD
text: In this section, we establish the global convergence results of the gradient descent method in implicit neural networks. In Section 3.1, we first study the dynamics of the prediction induced by the gradient flow, that is, the gradient descent with infinitesimal step-size. We show that the dynamics of the prediction is controlled by a Gram matrix whose smallest eigenvalue is strictly positive over iterations with high probability. Based on the findings in gradient flow, we will show that the random initialized gradient descent method with a constant step-size converges to a global minimum at a linear rate in Section 3.2.

heading: 3.1 CONTINUOUS TIME ANALYSIS
text: The gradient flow is equivalent to the gradient descent method with an infinitesimal step-size. Thus, the analysis of gradient flow can serve as a stepping stone towards understanding discrete-time gradient-based algorithms. Following previous works on gradient flow of different machine learning models (Saxe et al., 2013; Du et al., 2019; Arora et al., 2018; Kawaguchi, 2021) the gradient flow of implicit neural networks is given by the following ordinary differential equations:
dvec (A) dt = −∂L(t) ∂A , dvec (W ) dt = −∂L(t) ∂W , du dt = −∂L(t) ∂u , dv dt = −∂L(t) ∂v , (13)
where L(t) represents the value of the objective function L(θ) at time t.
Our results relies on the analysis of the dynamics of prediction ŷ(t). In particular, Lemma 3.1 shows that the dynamics of prediction ŷ(t) is governed by the spectral property of a Gram matrixH(t). Lemma 3.1 (Dynamics of Prediction ŷ(t)). Suppose ‖A(s)‖ ≤ c √ m for all 0 ≤ s ≤ t. Let X ∈ Rn×d, Φ(s) ∈ Rn×m, andZ(s) ∈ Rn×m be the matrices whose rows are the training data xi, feature vectors φi, and equilibrium points zi at time s, respectively. With scalar γ , min{γ0, γ0/c} for any γ0 ∈ (0, 1), the dynamics of prediction ŷ(t) is given by
dŷ dt =−
[ (γ2M(t) + In) ◦Z(t)Z(t)T +Q(t) ◦XXT + Φ(t)Φ(t)T ] (ŷ − y),
,−H(t)(ŷ − y), (14) where ◦ is the Hadamard product (i.e., element-wise product), and matrices M(t) ∈ Rn×n and Q(t) ∈ Rn×n are defined as follows:
M(t)ij , 1
m uT (Im − γ̃DiA)−1DiDTj (Im − γ̃DjA)−Tu, (15)
Q(t)ij , 1
m
( DTi (Im − γ̃DiA)−1u+ v )T EiE T j ( DTj (Im − γ̃DjA)−Tu+ v ) . (16)
Note that the matrixH(t) is clearly positive semidefinite since it is the sum of three positive semidefinite matrices. If there exists a constant λ > 0 such that λmin{H(t)} ≥ λ > 0 for all t, i.e., H(t) is positive definite, then the dynamics of the loss function L(t) satisfies the following inequality
L(t) ≤ exp{−λt}L(0), which immediately indicates that the objective value L(t) is consistently decreasing to zero at a geometric rate. With random initialization, we will show that H(t) is positive definite as long as the number of parameters m is sufficiently large and no two data points are parallel to each other. In particular, by using the nonlinearity of the feature map φ, we will show that the smallest eigenvalue of the Gram matrix G(t) , Φ(t)Φ(t)T is strictly positive over all time t with high probability. As a result, the smallest eigenvalue ofH(t) is always strictly positive.
Clearly, G(t) is a time-varying matrix. We first analyze its spectral property at its initialization. When t = 0, it follows from the definition of the feature vector φ in Eq. (2) that
G(0) = Φ(0)Φ(0)T = 1
m σ(XW (0)T )σ(XW (0)T )T =
1
m m∑ r=1 σ(Xwr(0))σ(Xwr(0)) T .
By Assumption 1, each vector wr(0) follows the standard multivariate normal distribution, i.e., wr(0) ∼ N(0, Id). By letting m→∞, we obtain the covariance matrixG∞ ∈ Rn×n as follows:
G∞ , Ew∼N(0,Id)[σ(Xw)σ(Xw) T ]. (17)
Here G∞ is a Gram matrix induced by the ReLU activation function and the random initialization. The following lemma shows that the smallest eigenvalue of G∞ is strictly positive as long as no two data points are parallel. Moreover, later in Lemmas 3.3 and 3.4, we conclude that the spectral property ofG∞ is preserved in the Gram matricesG(0) andG(t) during the training, as long as the number of parameter m is sufficiently large.
Lemma 3.2. Assume ‖xi‖ = 1 for all i ∈ [n]. If xi 6‖ xj for all i 6= j, then λ0 , λmin{G∞} > 0.
Proof. The proof follows from the Hermite Expansions of the matrixG∞ and the complete proof is provided in Appendix A.4.
Assumption 2 (Training Data). Without loss of generality, we can assume that eachxi is normalized to have a unit norm, i.e., ‖xi‖ = 1, for all i ∈ [n]. Moreover, we assume xi 6‖ xj for all i 6= j.
In most real-world datasets, it is extremely rare that two data samples are parallel. If this happens, by adding some random noise perturbation to the data samples, Assumption 2 can still be satisfied. Next, we show that at the initialization, the spectral property of G∞ is preserved in G(0) if m is sufficiently large. Specifically, the following lemma shows that if m = Ω̃(n2), then G(0) has a strictly positive smallest eigenvalue with high probability. The proof follows from the standard concentration bound for Gaussian random variables, and we relegate the proof to Appendix A.5.
Lemma 3.3. Let λ0 = λmin(G∞) > 0. If m = Ω ( n2 λ20 log ( n δ )) , then with probability of at least 1− δ, it holds that ‖G(0)−G∞‖2 ≤ λ04 , and hence λmin(G(0)) ≥ 3 4λ0.
During training, G(t) is time-varying, but it can be shown to be close to G(0) and preserve the spectral property of G∞, if the matrices W (t) has a bounded operator norm and it is not far away fromW (0). This result is formally stated below and its proof is provided in Appendix A.6. Lemma 3.4. Suppose ‖W (0)‖ ≤ c √ m, and λmin{G(0)} ≥ 34λ0. For any matrix W ∈ R m×d that satisfies ‖W ‖ ≤ 2c √ m and ‖W −W (0)‖ ≤ √ mλ0
16c‖X‖2 , R, the matrix defined by G , 1 mσ(XW T )σ(XW T )T satisfies ‖G−G(0)‖ ≤ λ04 and λmin(G) ≥ λ0 2 .
The next lemma shows three facts: (1) The smallest eigenvalue of G(t) is strictly positive for all t ≥ 0; (2) The objective value L(t) converges to zero at a linear convergence rate; (3) The (operator) norms of all trainable parameters are upper bounded by some constants, which further implies that unique equilibrium points in matrices Z(t) always exist. We prove this lemma by induction, which is provided in Appendix A.7. Lemma 3.5. Suppose that ‖u(0)‖ = √ m, ‖v(0)‖ = √ m, ‖W (0)‖ ≤ c √ m, ‖A(0)‖ ≤ c √ m, and
λmin{G(0)} ≥ 34λ0 > 0. If m = Ω ( c2n‖X‖2 λ30 ‖ŷ(0)− y‖2 ) and 0 < γ ≤ min{ 12 , 1 4c}, then for any t ≥ 0, the following results hold:
(i) λmin(G(t)) ≥ λ02 ,
(ii) ‖u(t)‖ ≤ 16c √ n
λ0 ‖ŷ(0)− y‖,
(iii) ‖v(t)‖ ≤ 8c √ n
λ0 ‖ŷ(0)− y‖,
(iv) ‖W (t)‖ ≤ 2c √ m,
(v) ‖A(t)‖ ≤ 2c √ m,
(vi) ‖ŷ(t)− y‖2 ≤ exp{−λ0t}‖ŷ(0)− y‖2.
By using simple union bounds in Lemma 2.1 and 3.5, we immediately obtain the global convergence result for the gradent flow as follows.
Theorem 3.1 (Convergence Rate of Gradient Flow). Suppose that Assumptions 1 and 2 hold. If we set the number of parameter m = Ω ( n2 λ20 log ( n δ )) and choose 0 < γ ≤ min{ 12 , 1 4c}, then with probability at least 1− δ over the initialization, we have ‖ŷ(t)− y‖2 ≤ exp{−λ0t}‖ŷ(0)− y‖2, ∀t ≥ 0.
This theorem establishes the global convergence of the gradient flow. Despite the nonconvexity of the objective function L(θ), Theorem 3.1 shows that if m is sufficiently large, then the objective value is consistently decreasing to zero at a geometric rate. In particular, Theorem 3.1 requires m = Ω̃(n2), which is similar or even better than recent results for the neural network with finite layers (Nguyen & Mondelli, 2020; Oymak & Soltanolkotabi, 2020; Allen-Zhu et al., 2019; Zou & Gu, 2019; Du et al., 2019). In particular, previous results showed that m is a polynomial of the number of training sample n and the number of layers h, i.e., m = Ω(nαhβ) with α ≥ 2 and β ≥ 12 (Nguyen & Mondelli, 2020, Table 1). These results do not apply in our case since we have infinitely many layers. By taking advantage of the nonlinear feature mapping function, we establish the global convergence for the gradient flow with m independent of depth h.

heading: 3.2 DISCRETE TIME ANALYSIS
text: In this section, we show that the randomly initialized gradient descent method with a fixed stepsize converges to a global minimum at a linear rate. With similar argument used in the analysis of gradient flow, we can show the (operator) norms of the training parameters are upper bounded by some constants. It is worth noting that, unlike the analysis of gradient flow, we do not have an explicit formula for the dynamics of prediction ŷ(t) in the discrete time analysis. Instead, we have to show the difference between the equilibrium pointsZ(k) in two consecutive iterations are bounded. Based on this, we can further bound the changes in the predictions ŷ(k) between two consecutive iterations. Another challenge is to show the objective value consistently decreases over iterations. A general strategy is to show the objective function is (semi-)smooth with respect to parameter θ, and apply the descent lemma to the objective function (Nguyen & Mondelli, 2020; Allen-Zhu et al., 2019; Zou & Gu, 2019). In this section, we take advantage of the nonlinear feature mapping function in Eq. (4). Consequently, we are able to obtain a Polyak-Łojasiewicz-like condition (Karimi et al., 2016; Nguyen, 2021), which allows us to provide a much simpler proof.
The following lemma establishes the global convergence for the gradient descent method with a fixed step-size when the operator norms ofA(0) andW (0) are bounded and λmin(G(0)) > 0. The proof is proved in Appendix A.8. Lemma 3.6 (Gradient Descent Convergence Rate). Suppose ‖u(0)‖ = √ m, ‖v(0)‖ = √ m, ‖W (0)‖ ≤ c √ m, ‖A(0)‖ ≤ c √ m, and λmin{G(0)} ≥ 34λ0 > 0. If 0 < γ ≤ min{ 1 2 , 1 4c},
m = Ω ( c2n‖X‖2
λ30 ‖ŷ(0)− y‖2
) , and stepsize α = O ( λ0/n 2 ) , then for any k ≥ 0, we have
(i) λmin(G(k)) ≥ λ02 ,
(ii) ‖u(k)‖ ≤ 32c √ n
λ0 ‖ŷ(0)− y‖,
(iii) ‖v(k)‖ ≤ 16c √ n
λ0 ‖ŷ(0)− y‖,
(iv) ‖W (k)‖ ≤ 2c √ m,
(v) ‖A(k)‖ ≤ 2c √ m,
(vi) ‖ŷ(k)− y‖2 ≤ (1− αλ0/2)k‖ŷ(0)− y‖2.
By using simple union bounds to combine Lemma 2.1 and 3.6, we obtain the global convergence result for the gradient descent. Theorem 3.2 (Convergence Rate of Gradient Descent). Suppose that Assumption 1 and 2 hold. If we set m = Ω ( n2 λ0 log ( n δ )) , 0 < γ ≤ min{ 12 , 1 4c}, and choose step-size α = O ( λ0/n 2 ) , then with probability at least 1− δ over the initialization, we have ‖ŷ(k)− y‖2 ≤ (1− αλ0/2)k‖ŷ(0)− y‖2, ∀k ≥ 0.

heading: 4 EXPERIMENTAL RESULTS
text: In this section, we use real-world datasets MNST, FashionMNST, CIFAR10, and SVHN to evaluate our theoretical findings. We initialize the entries of parameters A, W , u, and v by standard Gaussian or symmetric Bernoulli distribution independently as suggested in Assumption 1. For each dataset, we only use classes 0 and 1, and 500 samples are randomly drawn from each class to generate the training dataset with n = 1000 samples. All data samples are converted to gray scale and resized into a 28 × 28 pixel image. We also normalize each data to have unit norm. If two parallel samples are observed, we add a random Gaussian noise perturbation to one of them. Thus, Assumption 2 is also satisfied. We run 1000 epochs of gradient descent with a fixed step-size. We test three metrics with different width m. We first test how the extent of over-parameterization
affects the convergence rates. Then, we test the relation between the extent of over-parameterization and the operator norms between matrixA(k) and its initialization. Note that the “operator norm” in the plots denotes the operator norm of the scaled matrix (γ/ √ m)‖A(k)‖. Third, we test the extent of over-parameterization and the performance of the trained neural network on the unseen test data. Similar to the training dataset, We randomly select 500 samples from each class as the test dataset.
From Figure 1, the figures in the first column show that as m becomes larger, we have better convergence rates. The figures in the second column show that as m becomes larger, the neural networks achieve lower test loss. The figures in the third column show that the operator norms are slightly larger for larger m but overall the operator norms are approximately equal to its initialization. The bell curve from the classical bias-variance trade-off does not appear. This opens the door to a new research direction in implicit neural networks in the analyses of generalization error and bias-variance trade-off.

heading: 5 RELATED WORKS
text: Implicit models has been explored explored by the deep learning community for decades. For example, Pineda (1987) and ALMEIDA (1987) studied implicit differentiation techniques for training recurrent dynamics, also known as recurrent back-propagation (Liao et al., 2018). Recently, there has been renewed interested in the implicit models in the deep learning community (El Ghaoui et al., 2019; Gould et al., 2019). For example, Bai et al. (2019) introduces an implicit neural network called deep equilibrium model for the for the task of sequence modeling. By using implicit ODE solvers, Chen et al. (2018) proposed neural ordinary differential equation as an implicit residual network with continuous-depth. Other instantiations of implicit modeling include optimization layers (Djolonga & Krause, 2017; Amos & Kolter, 2017), differentiable physics engines (de Avila Belbute-Peres et al., 2018; Qiao et al., 2020), logical structure learning (Wang et al., 2019), and continuous generative models (Grathwohl et al., 2019).
The theoretical study of training finite-layer neural networks via over-parameterization has been an active research area. Jacot et al. (2018) showed the trajectory of the gradient descent method can be characterized by a kernel called neural tangent kernel for smooth activation and infinitely width neural networks. For a finite-width neural network with smooth or ReLU activation, Arora et al. (2019); Du et al. (2019); Li & Liang (2018) showed that the dynamics of the neural network is governed by a Gram matrix. Consequently, Zou et al. (2020); Du et al. (2019); Allen-Zhu et al. (2019); Nguyen & Mondelli (2020); Arora et al. (2019); Zou & Gu (2019); Oymak & Soltanolkotabi (2020) showed that (stochastic) gradient descent can attain global convergence for training a finite-layer neural network when the width m is a polynomial of the sample size n and the depth h. However, their results cannot be applied directly to implicit neural networks since implicit neural networks have infinite layers and the equilibrium equation may not be well-posed. Our work establishes the well-posedness of the equilibrium equation even if the width m is only square of the sample size n.

heading: 6 CONCLUSION AND FUTURE WORK
text: In this paper, we provided a convergence theory for implicit neural networks with ReLU activation in the over-parameterization regime. We showed that the random initialized gradient descent method with fixed step-size converges to a global minimum of the loss function at a linear rate if the width m = Ω̃(n2). In particular, by using a fixed scalar γ ∈ (0, 1) to scale the random initialized weight matrixA, we proved that the equilibrium equation is always well-posed throughout the training. By analyzing the gradient flow, we observe that the dynamics of the prediction vector is controlled by a Gram matrix whose smallest eigenvalue is lower bounded by a strictly positive constant as long as m = Ω̃(n2). We envision several potential future directions based on the observations made in the experiments. First, we believe that our analysis can be generalized to implicit neural networks with other scaling techniques and initialization. Here, we use a scalar γ/ √ m to ensure the existence of the equilibrium point z∗ with random initialization. With an appropriate normalization, the global convergence for identity initialization can be obtained. Second, we believe that the width m not only improves the convergence rates but also the generalization performance. In particular, our experimental results showed that with a larger m value, the test loss is reduced while the classical bell curve of bias-variance trade-off is not observed.

heading: ACKNOWLEDGMENTS
text: This work has been supported in part by National Science Foundation grants III-2104797, DMS1812666, CAREER CNS-2110259, CNS-2112471, CNS-2102233, CCF-2110252, CNS-2120448, CCF-19-34884 and a Google Faculty Research Award.

heading: A APPENDIX

heading: A.1 PROOF OF LEMMA 2.2
text: Proof. Let γ0 ∈ (0, 1). Set γ , min{γ0, γ0/c}. Denote γ̃ = γ/ √ m. Then∥∥z`+1 − z`∥∥ =∥∥σ (γ̃Az` + φ)− σ (γ̃Az`−1 + φ)∥∥
≤γ̃ ∥∥Az` −Az`−1∥∥ , σ is 1-Lipschitz continuous
=γ̃ ∥∥A(z` − z`−1)∥∥
≤γ̃‖A‖‖z` − z`−1‖, ≤γ̃c √ m‖z` − z`−1‖, ‖A‖ ≤ c √ m
=γ0‖z` − z`−1‖.
Applying the above argument ` times, we obtain
‖z`+1 − z`‖ ≤ γ`0‖z1 − z0‖ = γ`0‖z1‖ = γ`0‖σ(φ)‖ ≤ γ`0‖φ‖,
where we use the fact z0 = 0. For any positive integers p, q with p ≤ q, we have
‖zp − zq‖ ≤‖zp − zp+1‖+ · · ·+ ‖zq−1 − zq‖ ≤γp0‖φ‖+ · · ·+ γ q 0‖φ‖
≤γp0‖φ‖ ( 1 + γ0 + γ 2 0 + · · · ) =
γp0 1− γ0 ‖φ‖.
Since γ0 ∈ (0, 1), we have ‖zp − zq‖ → 0 as p → ∞. Hence, {z`}∞`=1 is a Cauchy sequence. Since Rm is complete, the equilibrium point z∗ is the limit of the sequence {z`}∞`=1, so that z exists and is unique. Moreover, let q →∞, then we obtain ‖zp − z∗‖ ≤ γ p
1−γ ‖φ‖, so that the fixed-point iteration converges to z linearly.
Let p = 0 and q = `, then we obtain ‖z`‖ ≤ 11−γ0 ‖φ‖.

heading: A.2 PROOF OF LEMMA 2.3
text: Proof. (i) To simplify the notations, we denote D , diag(σ′(γ̃Az + φ)), and E , diag(σ′(Wx)). The differential of f is given by
df =d(z − γ̃σ(γ̃Az + φ)) =dz −Dd(γ̃Az + φ) = [Im − γ̃DA] dz − γ̃D(dA)z −Ddφ.
Taking vectorization on both sides yields
vec (df) = [Im − γ̃DA] vec (dz)− vec (γ̃DdAz)−Dvec (dφ) = [Im − γ̃DA] vec (dz)− γ̃[zT ⊗D]vec (dA)−Dvec (dφ) .
Therefore, the partial derivative of f with respect to z,A, and φ are given by
∂f ∂z = [Im − γ̃DA]T ,
∂f ∂A = −γ̃
[ zT ⊗D ]T ,
∂f ∂φ = −DT .
It follows from the definition of the feature vector φ in Eq. (2) that
dφ = 1√ m dσ(Wx) = 1√ m E(dW )x = 1√ m
[ xT ⊗E ] vec (dW ) .
Thus, the partial derivative of φ with respect toW is given by ∂φ
∂W = 1√ m
[ xT ⊗E ]T . (18)
By using the chain rule, we obtain the partial derivative of f with respect toW as follows ∂f
∂W =
∂φ
∂W
∂f ∂φ = − 1√ m
[ xT ⊗E ]T DT .
(ii) Let v be an arbitrary vector, and u be an arbitrary unit vector. The reverse triangle inequality implies that
‖(Im − γ̃ diag(σ′(v))A)u‖ ≥‖u‖ − ‖γ̃ diag(σ′(v))Au‖ ≥‖u‖ − γ̃‖diag(σ′(v))‖‖A‖‖u‖ (a)
≥ (1− γ0)‖u‖ =1− γ0 > 0,
where (a) is due to |σ′(v)| ≤ 1 and ‖A‖op ≤ c √ m. Therefore, taking infimum on the left-hand side over all unit vector u yields the desired result.
(iii) Since f(z∗,A,W ) = 0, taking implicit differentiation of f with respect toA at z∗ gives us( ∂z
∂A ∣∣∣∣ z=z∗ )( ∂f ∂z ∣∣∣∣ z=z∗ ) + ( ∂f ∂A ∣∣∣∣ z=z∗ ) = 0.
The results in part (i)-(ii) imply the smallest eigenvalue of ∂f∂z ∣∣ z∗
is strictly positive, so that it is invertible. Therefore, we have
∂z∗ ∂A = −
( ∂f
∂A ∣∣∣∣ z=z∗ )( ∂f ∂z ∣∣∣∣ z=z∗ )−1 = γ̃ [ zT ⊗D ]T [Im − γ̃DA]−T . (19)
Similarly, we obtain the partial derivative of z∗ with respect toW as follows
∂z∗ ∂W = −
( ∂f
∂W ∣∣∣∣ z=z∗ )( ∂f ∂z ∣∣∣∣ z=z∗ )−1 = 1√ m [ xT ⊗E ]T DT [Im − γ̃DA]−T . (20)
To further simplify the notation, we denote z to be the equilibrium point z∗ by omitting the superscribe, i.e., z = z∗. Let ŷ = uTz + vTφ be the prediction for the training data (x,y). The differential of ŷ is given by
dŷ = d ( uTz + vTφ ) = uT dz + zdu+ vT dφ+ φT dv.
The partial derivative of ŷ with respect to u, v, z, and φ are given by ∂ŷ
∂z = u,
∂ŷ ∂u = z, ∂ŷ ∂v = φ, ∂ŷ ∂φ = v. (21)
Let ` = 12 (ŷ − y) 2. Then ∂`/∂ŷ = (ŷ − y). By chain rule, we have
∂` ∂u = ∂ŷ ∂u ∂` ∂ŷ = z(ŷ − y), (22)
∂` ∂φ = ∂ŷ ∂v ∂` ∂ŷ = φ(ŷ − y). (23)
By using (19)-(20) and chain rule, we obtain ∂`
∂A = ∂z ∂A ∂` ∂z
= ∂z
∂A
∂ŷ
∂z
∂` ∂ŷ = γ̃(ŷ − y)
[ zT ⊗D ]T [Im − γ̃DA]−T u, (24)
and
∂` ∂W = ∂z ∂W ∂ŷ ∂z ∂` ∂ŷ + ∂φ ∂W ∂ŷ ∂φ ∂` ∂ŷ
= 1√ m
(ŷ − y)[xT ⊗E]T [ DT (Im − γ̃DA)−Tu+ v ] . (25)
Since L = ∑n i=1 `i with `i = `(ŷi, yi), we have dL = ∑n i=1 d`i and ∂L/∂`i = 1. Therefore, we obtain
∂L ∂A = n∑ i=1 ∂`i ∂A = n∑ i=1 γ̃(ŷi − yi) [ zTi ⊗Di ]T [Im − γ̃DiA]−T u, (26)
∂L
∂W = n∑ i=1 ∂`i ∂W = n∑ i=1 1√ m (ŷi − yi)[xTi ⊗Ei]T [ DTi (Im − γ̃DiA)−Tu+ v ] , (27)
∂L ∂u = n∑ i=1 ∂`i ∂u = n∑ i=1 (ŷi − yi)zi, (28)
∂L ∂v = n∑ i=1 ∂`i ∂v = n∑ i=1 (ŷi − yi)φi. (29)

heading: A.3 PROOF OF LEMMA 3.1
text: Proof. Let zi denote the i-th equilibrium point for the i-the data sample xi. By using (19), (20), (26) and (27), we obtain the dynamics of the equilibrium point zi as follows
dzi dt = ( ∂zi ∂A )T dvec (A) dt + ( ∂zi ∂W )T dvec (W ) dt
= ( ∂zi ∂A )T ( − ∂L ∂A ) + ( ∂zi ∂W )T ( − ∂L ∂W ) =− γ̃2
n∑ j=1 (ŷj − yj) [Im − γ̃DiA]−1 [zTi ⊗Di] [ zTj ⊗Dj ]T [Im − γ̃DjA]−T u
− 1 m n∑ j=1 (ŷj − yj) [Im − γ̃DiA]−1Di [ xTi ⊗Ei ] [xTj ⊗Ej ]T [ DTj (Im − γ̃DjA)−Tu+ v ] =− γ̃2
n∑ j=1 (ŷj − yj) [Im − γ̃DiA]−1DiDTj [Im − γ̃DjA] −T uzTi zj
− 1 m n∑ j=1 (ŷj − yj) [Im − γ̃DiA]−1DiEiETj [ DTj (Im − γ̃DjA)−Tu+ v ] xTi xj .
By using (18) and 27, we obtain the dynamics of the feature vector φi
dφi dt = ( ∂φi ∂W )T dvec (W ) dt
= ( ∂φi ∂W )T ( − ∂L ∂W ) =− 1
m n∑ j=1 (ŷi − yi)EiETj [DTj (Im − γ̃DjA)−Tu+ v]xTi xj .
By chain rule, the dynamics of the prediction ŷi is given by
dŷi dt = ( ∂ŷi ∂zi )T dzi dt + ( ∂ŷi ∂φi )T dφi dt + ( ∂ŷi ∂u )T du dt + ( ∂ŷi ∂v )T dv dt
=− γ̃2 n∑ j=1 (ŷj − yj) [ uT (Im − γ̃DiA)−1DiDTj (Im − γ̃DjA)−Tu ] (zTi zj)
− 1 m n∑ j=1 (ŷj − yj) [( DTi (Im − γ̃DiA)−1u+ v )T EiE T j ( DTj (Im − γ̃DjA)−Tu+ v )] (xTi xj)
− n∑ j=1 (ŷj − yj)(zTi zj)
− n∑ j=1 (ŷj − yj)(φTi φj).
Define the matricesM(t) ∈ Rn×n andQ(t) ∈ Rn×n as follows
M(t)ij , 1
m uT (Im − γ̃DiA)−1DiDTj (Im − γ̃DjA)−Tu,
Q(t)ij , 1
m
( DTi (Im − γ̃DiA)−1u+ v )T EiE T j ( DTj (Im − γ̃DjA)−Tu+ v ) .
Let X ∈ Rn×d, Φ(t) ∈ Rn×m, and Z(t) ∈ Rn×m be the matrices whose rows are the training data xi, feature vectors φi, and equilibrium points zi at time t, respectively. The dynamics of the prediction vector ŷ is given by
dŷ dt = −
[( γ2M(t) + In ) ◦Z(t)Z(t)T +Q(t) ◦XXT + Φ(t)Φ(t)T ] (ŷ(t)− y).

heading: A.4 PROOF OF LEMMA 3.2

heading: A.4.1 REVIEW OF HERMITE EXPANSIONS
text: To make the paper self-contained, we review the necessary background about the Hermite polynomials in this section. One can find each result in this section from any standard textbooks about functional analysis such as MacCluer (2008); Kreyszig (1978), or most recent literature (Nguyen & Mondelli, 2020, Appendix D) and (Oymak & Soltanolkotabi, 2020, Appendix H).
We consider an L2-space defined by L2(R, dP ), where dP is the Gaussian measure, that is,
dP = p(x)dx, where p(x) = 1√ 2π e− x2 2 .
Thus, L2(R, dP ) is a collection of functions f for which∫ ∞ −∞ |f(x)|2 dP (x) = ∫ ∞ −∞ |f(x)|2 p(x)dx = Ex∼N(0,1) |f(x)| 2 <∞.
Lemma A.1. The ReLU activation σ ∈ L2(R, dP ).
Proof. Note that∫ ∞ −∞ |σ(x)|2 p(x)dx ≤ ∫ ∞ −∞ |x|2 p(x)dx = Ex∼N(0,1) |x| 2 = Var(x) = 1.
For any functions f, g ∈ L2(R, dP ), we define an inner product
〈f, g〉 := ∫ ∞ −∞ f(x)g(x)dP (x) = ∫ ∞ −∞ f(x)g(x)p(x)dx = Ex∼N(0,1)[f(x)g(x)].
Furthermore, the induced norm ‖ · ‖ is given by
‖f‖2 = 〈f, f〉 = ∫ ∞ −∞ |f(x)|2 dP (x) = Ex∼N(0,1) |f(x)| 2 .
This L2 space has an orthonormal basis with respect to the inner product defined above, called normalized probabilist’s Hermite polynomials {hn(x)}∞n=0 that are given by
hn(x) = 1√ n! (−1)nex 2/2Dn(e−x 2/2), where Dn(e−x 2/2) =
dn
dxn e−x 2/2.
Lemma A.2. The normalized probabilist’s Hermite polynomials is an orthonormal basis of L2(R, dP ): 〈hm, hn〉 = δmn.
Proof. Note that Dn(e−x 2/2) = e−x 2/2Pn(x) for a polynomial with degree of n and leading term is (−1)nxn. Thus, we can consider hn(x) = 1√n! (−1) nPn(x).
Assume m < n
〈hn, hm〉 =Ex∼N(0,1)[hn(x)hm(x)]
= ∫ ∞ −∞ hn(x)hm(x) 1√ 2π e−x 2/2dx,
= 1√
2π √ n! (−1)n ∫ ∞ −∞ Dn(e−x 2/2)hm(x)dx, rewrite hn(x) by its definition
= 1√
2π √ n! √ m! (−1)n+m ∫ ∞ −∞ Dn(e−x 2/2)Pm(x)dx, rewrite hm by the polynomial form
= 1√
2π √ n! √ m! (−1)2n+m ∫ ∞ −∞ e−x 2/2Dn[Pm(x)]dx, integration by parts n times
There is no boundary terms because the super exponential decay of e−x 2/2 at infinity. Since m < n, then Dn(Pm) = 0 so that 〈hm, hn〉=0. If m = n, then Dn(Pm) = (−1)nn!. Thus, 〈hn, hn〉 = 1.
Remark: Since {hn} is an orthonormal basis, for every f ∈ L2(R, dP ), we have
f(x) = ∞∑ n=0 〈f, hn〉hn(x)
in the sense that
lim N→∞ ∥∥∥∥∥f(x)− N∑ n=0 〈f, hn〉hn(x) ∥∥∥∥∥ 2 = lim N→∞ Ex∼N(0,1) ∣∣∣∣∣f(x)− N∑ n=0 〈f, hn〉hn(x) ∣∣∣∣∣ 2 = 0
Lemma A.3. f ∈ L2(R, dP ) if and only if ∑∞ n=0 |〈f, hn〉| 2 <∞.
Proof. Note that 〈f, f〉 = ∫ ∞ −∞ |f(x)|2 dP (x)
= ∫ ∞ −∞ ( ∞∑ i=0 〈f, hi〉hi(x) ) ∞∑ j=0 〈f, hj〉hj(x)  dP (x) =
∞∑ i,j=0 〈f, hi〉 〈f, hj〉 ∫ ∞ −∞ hi(x)hj(x)dP (x)
= ∞∑ i=1 |〈f, hi〉|2 .
Lemma A.4. Consider a Hilbert spaceH with inner product 〈·, ·〉. If ‖fn−f‖ → 0 and ‖gn−g‖ → 0, then 〈f, g〉 = limn→∞ 〈fn, gn〉.
Proof. Observe that
|〈f, g〉 − 〈fn, gn〉| ≤ |〈f, g〉 − 〈fn, g〉|+ |〈fn, g〉 − 〈fn, gn〉| ≤‖f‖‖g − gn‖+ ‖fn‖‖g − gn‖.
Let n→∞, then the continuity of ‖ · ‖ implies the desired result.
Lemma A.5. Let {hn(x)} be the normalized probabilist’s Hermite polynomials. For any fixed number t, we have
ext−t 2/2 = ∞∑ n=0 tn√ n! hn(x). (30)
Proof. First, we show f(x) = ext−t 2/2 ∈ H , L2(R, dP ).
〈f, f〉 =Ex∼N(0,1) |f(x)| 2
= ∫ ∞ −∞ e2xt−t 2 1√ 2π e−x 2/2dx
=et 2 ∫ ∞ −∞ 1√ 2π exp { − (x− 2t) 2 2 } dx, x ∼ N(2t, 1)
=et 2 <∞.
Thus f(x) ∈ H . Then f(x) = ∑∞ n=0 〈f, hn〉hn(x). Note that
〈f, hn〉 =Ex∼N(0,1)[f(x)hn(x)]
= ∫ ∞ −∞ ext−t 2/2 · 1√ n! (−1)nex 2/2Dn(e −x2/2) · 1√ 2π e−x 2/2dx
= 1√ n! (−1)n 1√ 2π ∫ ∞ −∞ ext−t 2/2 ·Dn(e−x 2/2)dx, integration by parts n times
= 1√ n! (−1)2n 1√ 2π ∫ ∞ −∞ ext−t 2/2tn · e−x 2/2dx = tn√ n! ∫ ∞ −∞ 1√ 2π e−(x−t) 2/2dx, x ∼ N(t, 1) = tn√ n! .
Lemma A.6. Let a, b ∈ Rd with ‖a‖ = ‖b‖ = 1, then Ew∼N(0,Id)[hn(〈a,w〉)hm(〈b,w〉)] = 〈a, b〉 n δmn.
Proof. Given fixed numbers s and t, we define two functions f(w) = e〈a,w〉t−t 2/2 and g(w) = e〈b,w〉s−s 2/2. Let x = 〈a,w〉 and y = 〈b, w〉. Then we have
f(w) =e〈a,w〉t−t 2/2 = ext−t 2/2 = ∞∑ n=0 tn√ n! hn(x) = ∞∑ n=0 tn√ n! hn(〈a,w〉),
g(w) =e〈b,w〉s−s 2/2 = eys−s 2/2 = ∞∑ n=0 sn√ n! hn(y) = ∞∑ n=0 sn√ n! hn(〈b,w〉).
Define a Hilbert space Hd = L2(Rd, dP ), where dP is the multivariate Gaussian measure, equipped with inner product 〈f, g〉 , Ew∼N(0,Id)[f(w)g(w)]. Clearly, f, g ∈ Hd. Define sequences {fN} and {gN} as follows
fN (w) = N∑ n=0 tn√ n! hn(〈a,w〉) and gN (w) = N∑ n=0 sn√ n! hn(〈b,w〉).
Since ‖f − fN‖ → 0 and ‖g − gN‖ → 0, we have Ew∼N(0,Id)[f(w)g(w)] = 〈f, g〉
= lim N→∞
〈fN , gN 〉
= lim N→∞ Ew∼N(0,Id)[fN (w)gN (w)]
= lim N→∞ N∑ n,m=0 tnsm√ n! √ m! Ew∼N(0,Id)[hn(〈a,w〉)gn(〈b,w〉)]
Note that the LHS is also given by
Ew∼N(0,Id)[f(w)g(w)] =e −t2/2−s2/2Ew∼N(0,Id)[e 〈a,w〉t+〈b,w〉s]
=e−t 2/2−s2/2Ew∼N(0,Id)[e
∑d i=1 wi(ait+bis)]
=e−t 2/2−s2/2 d∏ i=1 Ewi∼N(0,1)[e wi(ait+bis)]
=e−t 2/2−s2/2 d∏ i=1 Mwi(ait+ bis) =e〈a,b〉st
= ∞∑ n=0 〈a, b〉n (st)n n! .
Since s and t are arbitrary numbers, matching the coefficients yields
Ew∼N(0,Id)[hn(〈a,w〉)hm(〈b,w〉)] = 〈a, b〉 n δmn.

heading: A.4.2 LOWER BOUND THE SMALLEST EIGENVALUES OF G∞
text: The result in this subsection is similar to the results in (Nguyen & Mondelli, 2020, Appendix D) and (Oymak & Soltanolkotabi, 2020, Appendix H). The key difference is the assumptions made on the training data. In particular, Oymak & Soltanolkotabi (2020) assumes the training data is δ-separable, i.e., min{‖xi − xj‖, ‖xi + xj‖} ≥ δ > 0 for all i 6= j, and Nguyen & Mondelli (2020) assumes the data xi follows some sub-Gaussian random variable, while we assume no two data are parallel to each other, i.e., xi 6‖ xj for all i 6= j.
Lemma A.7. Given an activation function σ, if σ ∈ L2(R, dP ) and ‖xi‖ = 1 for all i ∈ [n], then
G∞ = ∞∑ k=0 |〈σ, hk〉|2 ( XXT ◦ · · · ◦XXT )︸ ︷︷ ︸ k times , (31)
where ◦ is elementwise product.
Proof. Observe
G∞ij =Ew∼N(0,Id) [σ(〈w,xi〉)σ(〈w,xj〉)]
= ∞∑ k,`=0 〈σ, hk〉 〈σ, h`〉Ew∼N(0,Id) [hk(〈w,xi〉)h`(〈w,xj〉)]
= ∞∑ k,`=0 〈σ, hk〉 〈σ, h`〉 · 〈xi,xj〉k δk`
= ∞∑ k=0 〈σ, hk〉2 〈xi,xj〉k
Note that the tensor product of xi and xi is xi ⊗ xi ∈ Rd 2×1, so that
〈xi,xj〉k = 〈 xi ⊗ · · · ⊗ xi︸ ︷︷ ︸
k times
,xj ⊗ · · · ⊗ xj︸ ︷︷ ︸ k times
〉
Here we introduce the (row-wise) Khatri–Rao product of two matrices A ∈ Rk×m, B ∈ Rk×n. Then
A ∗B = A1∗ ⊗B1∗... Ak∗ ⊗Bk∗  ∈ Rk×mn where Ai∗ indicates the i-th row of matrix A. Therefore, the i-th row of X ∗ · · · ∗X , X∗n is xi ⊗ · · · ⊗ xi. As a result, we obtain a more compact form of (31) as follows
G∞ = ∞∑ k=0 |〈σ, hk〉|2 (X∗k)(X∗k)T . (32)
Lemma A.8. If σ(x) is a nonlinear function and |σ(x)| ≤ |x| and , then
sup{n : 〈σ, hn〉 > 0} =∞.
Proof. It is equivalent to show σ(x) is not a finite linear combination of polynomials. We prove by contradiction. Suppose σ(x) = a0 + a1x + · · · + anxn. Since σ(0) = 0 = a0, then σ(x) = a1x+ · · ·+ anxn. Observe that
lim x→∞ |σ(x)| |x| = lim x→∞ |a1x+ · · ·+ anxn| |x|
= lim x→∞ ∣∣a1 + · · ·+ anxn−1∣∣ , =∞
which contradicts |σ(x)||x| ≤ 1 for all x 6= 0.
Lemma A.9. If xi 6‖ xj for all i 6= j, then there exists k0 > 0 such that λmin [ (X∗k)(X∗k)T ] > 0 for all k ≥ k0. Therefore, λmin(G∞) > 0.
Proof. To simplify the notation, denoteK = (X∗k)T ∈ Rkd×n. Since xi 6‖ xj and ‖xi‖ = 1, then let δ , max{|〈xi,xj〉|} = max{|cos θij |} and δ ∈ (0, 1), where θij is the angle between xi and xj . For any unit vector v ∈ Rn, we have
vT (X∗k)(X∗k)Tv =‖Kv‖2 = ∥∥∥∥∥ n∑ i=1 viK∗i ∥∥∥∥∥ 2
= n∑ i=1 n∑ j=1 vivj 〈K∗i,K∗j〉
= n∑ i=1 n∑ j=1 vivj 〈xi,xj〉k
= n∑ i=1 v2i ‖xi‖2k + ∑ i6=j vivj 〈xi,xj〉k
=1 + ∑ i 6=j vivj 〈xi,xj〉k ,
where the last equality is because ‖xi‖ = 1 and ‖v‖ = 1. Note that∣∣∣∣∣∣ ∑ i 6=j vivj 〈xi,xj〉k ∣∣∣∣∣∣ ≤ ∑ i 6=j |vi| |vj | |〈xi,xj〉|k
≤δk ∑ i6=j |vi| |vj | , by |〈xi,xj〉| ≤ δ
≤δk (
n∑ i=1 |vi| )2 ≤nδk, by Cauchy-Schwart’s inequlity.
By inverse triangle inequality, we have
‖Kv‖2 ≥ 1− nδk.
Choose k0 ≥ log n/ log(1/δ), then λmin{(X∗k)(X∗k)T } > 0 for all k ≥ k0.

heading: A.5 PROOF OF LEMMA 3.3
text: Proof. Since ‖xi‖ = 1 and wr(0) ∼ N (0, Id), we have xTi wr(0) ∼ N(0, 1) for all i ∈ [n]. Let Xir , σ [ xTi wr(0) ] and Z ∼ N(0, 1), then for any |λ| ≤ 1/ √ 2, we have
E exp{X2irλ2} = E exp{σ [ xTi wr(0) ]2 λ2} ≤ E exp{Z2λ2} = 1/ √ 1− 2λ2 ≤ e2t 2 ,
where the first inequality is due tot |σ(x)| ≤ |x|, and the last inequality is by using the numerical inequality 1/(1 − x) ≤ e2x. Choose λ ≤ ( log √ 2 )1/2
, we obtain E{X2irλ2} ≤ 2. By using Markov’s inequality, we have for any t ≥ 0
P {|Xir| ≥ t} = P { |Xir|2 / log √ 2 ≥ t2/ log √ 2 } ≤ 2 exp { −t2 log √ 2 } ≤ 2 exp { −t2/4 } .
Therefore Xir is a sub-Gaussian random variable with sub-Gaussian norm ‖Xi‖ψ2 ≤ 2 (Vershynin, 2018, Proposition 2.5.2). Then XirXjr is a sub-exponential random variable with sub-exponential norm ‖XirXjr‖ψ1 ≤ 4 (Vershynin, 2018, Lemma 2.7.7). Observe that
Gij(0) = φi(0) Tφj(0) =
1
m m∑ r=1 σ [ xTi wr(0) ] σ [ xTj wr(0) ] = 1 m m∑ r=1 XirXjr.
Since G∞ij = E [Gij(0)], (Vershynin, 2018, Exercise 2.7.10) implies that Gij(0) − G∞ij is also a zero-mean sub-exponential random variable. It follows from the Bernstein’s inequality that
P { ‖G(0)−G∞‖2 ≥
λ0 4
} ≤P { ‖G(0)−G∞‖F ≥
λ0 4 } =P { ‖G(0)−G∞‖2F ≥ ( λ0 4 )2}
=P  n∑
i,j=1
∣∣Gij(0)−G∞ij ∣∣2 ≥ (λ04 )2
≤ n∑
i,j=1
P {∣∣Gij(0)−G∞ij ∣∣2 ≥ (λ04n )2}
= n∑ i,j=1 P {∣∣Gij(0)−G∞ij ∣∣ ≥ λ04n } ≤n2 · 2 exp { −cλ20m/n2
} ≤δ,
where c > 0 is some constant, and we use the facts ‖X‖2 ≤ ‖X‖F , and P{ ∑n i=1 xi ≥ ε} ≤∑n
i=1 P{xi ≥ ε/n}.

heading: A.6 PROOF OF LEMMA 3.4
text: Proof. By using the 1-Lipschitz continuity of σ(x), we have
‖G−G(0)‖ = 1 m ‖σ(XW T )σ(XW T )T − σ(XW (0)T )σ(XW (0)T )T ‖
≤ 1 m ‖σ(XW T )σ(XW T )T − σ(XW T )σ(XW (0)T )T ‖
+ 1
m ‖σ(XW T )σ(XW (0)T )T − σ(XW (0)T )σ(XW (0)T )T ‖
= 1
m ‖σ(XW T )‖‖σ(XW T )− σ(XW (0)T )‖ + 1
m ‖σ(XW T )− σ(XW (0)T )‖‖σ(XW (0)T )‖
≤ 1 m ‖X‖‖W ‖‖X‖‖W −W (0)‖+ 1 m ‖X‖‖W −W (0)‖‖X‖‖W (0)‖ ≤ 4c√ m ‖X‖2‖W −W (0)‖
≤λ0 4 .

heading: A.7 PROOF OF LEMMA 3.5
text: Proof. It suffices to show the result holds for γ = min{γ0, γ0/2c}, where γ0 = 1/2. Note that Lemma 2.1 still holds if one chooses a larger c. Thus, we choose c % √ λ0/‖X‖. We prove by the induction. Suppose that for 0 ≤ s ≤ t, the followings hold
(i) λmin(G(s)) ≥ λ02 ,
(ii) ‖u(s)‖ ≤ 16c √ n
λ0 ‖ŷ(0)− y‖,
(iii) ‖v(s)‖ ≤ 8c √ n
λ0 ‖ŷ(0)− y‖,
(iv) ‖W (s)‖ ≤ 2c √ m,
(v) ‖A(s)‖ ≤ 2c √ m,
(vi) ‖ŷ(s)− y‖2 ≤ exp{−λ0s}‖ŷ(0)− y‖2,
Since λmin(G(s)) ≥ λ02 , we have
d dt ‖ŷ(t)− y‖2 =− 2(ŷ(t)− y)TH(t)(ŷ(t)− y)
≤− λ0‖ŷ(t)− y‖2
Solving the ordinary differential equation yields
‖ŷ(t)− y‖2 ≤ exp{−λ0t}‖ŷ(0)− y‖2.
By using the inductive hypothesis ‖W (s)‖ ≤ 2c √ m, we have
‖φi(s)‖ = ∥∥∥∥ 1√mσ(W (s)xi) ∥∥∥∥ ≤ 1√m‖W (s)‖‖xi‖ ≤ 2c. It follows from Lemma 2.2 with γ0 = 1/2 that
‖z∗i (s)‖ ≤ 2‖φi(s)‖ ≤ 4c.
Note that
‖∇vL(s)‖ ≤ n∑ i=1 |ŷi(s)− yi| ‖φi(s)‖
≤2c n∑ i=1 |ŷi(s)− yi| ≤2c √ n‖ŷ(s)− y‖ ≤2c √ n exp{−λ0s/2}‖y(0)− y‖
and so
‖v(t)− v(0)‖ ≤ ∫ t 0 ‖∇vL(s)‖ds
≤2c √ n‖y(0)− y‖ ∫ t 0 exp{−λ0s/2}ds ≤4c √ n
λ0 ‖ŷ(0)− y‖,
Since vi(0) follows symmetric Bernoulli distribution, then ‖v(0)‖ = √ m and we obtain
‖v(t)‖ ≤ ‖v(t)− v(0)‖+ ‖v(0)‖ ≤ 8c √ n
λ0 ‖ŷ(0)− y‖,
where the last inequality is due to m = Ω ( c2n‖X‖2
λ30 ‖ŷ(0)− y‖2
) and c % √ λ0/‖X‖.
Similarly, we have
‖∇uL(s)‖ ≤ n∑ i=1 |ŷi(s)− yi| ‖z∗i ‖
≤4c √ n‖ŷ(s)− y‖ ≤4c √ n exp{−λ0s/2}‖ŷ(0)− y‖,
so that
‖u(t)− u(0)‖ ≤ ∫ t 0 ‖∇uL(s)‖ds ≤ 8c √ n λ0 ‖ŷ(0)− y‖.
Since ui(0) follows symmetric Bernoulli distribution, then ‖u(0)‖ = √ m and we obtain
‖u(t)‖ ≤ ‖u(t)− u(0)‖+ ‖u(0)‖ ≤ 16c √ n
λ0 ‖ŷ(0)− y‖.
Note that ‖∇WL(s)‖ ≤ n∑ i=1 1√ m |ŷi(s)− yi| ‖Ei(s)‖ ( ‖Ui(s)−1u(s)‖+ ‖v(s)‖ ) ‖xi‖
≤64c √ n
λ0 ‖ŷ(0)− y‖. n∑ i=1 |ŷi(s)− yi|
≤64cn λ0 ‖ŷ(0)− y‖ · ‖ŷ(s)− y‖ ≤64cn λ0 ‖ŷ(0)− y‖2 · exp{−λ0s/2},
so that
‖W (t)−W (0)‖ ≤ ∫ t 0 ‖∇WL(s)‖ds
≤ 128cn λ20 √ m ‖ŷ(0)− y‖2 ≤ λ0 √ m
16c‖X‖2
≤R. Therefore, we obtain
‖W (t)‖ ≤ ‖W (t)−W (0)‖+ ‖W (0)‖ ≤ 2c √ m,
Moreover, it follows from Lemma 3.4 that λmin{G(t)} ≥ λ02 . Note that
‖∇AL(s)‖ ≤ n∑ i=1 γ√ m |ŷi(s)− yi| ‖Di‖‖Ui(s)−1‖‖u(s)‖‖z∗i ‖
≤32c √ n
λ0 √ m ‖ŷ(0)− y‖ · n∑ i=1 |ŷi(s)− yi|
≤ 32cn λ0 √ m ‖ŷ(0)− y‖ · ‖ŷ(s)− y‖
≤ 32cn λ0 √ m ‖ŷ(0)− y‖2 · exp{−λ0s/2},
so that
‖A(t)−A(0)‖ ≤ ∫ t 0 ‖∇AL(s)‖ds
≤ 64cn λ20 √ m ‖ŷ(0)− y‖2.
Then
‖A(t)‖ ≤ ‖A(t)−A(0)‖+ ‖A(0)‖ ≤ 2c √ m.

heading: A.8 PROOF OF LEMMA 3.6
text: In this section, we prove the result for discrete time analysis or result for gradient descent. Assume ‖A(0)‖ ≤ c √ m and ‖W (0)‖ ≤ c √ m. Further, we assume λmin(G(0)) ≥ 34λ0 and we assume m = Ω ( c2n‖X‖2
λ30 ‖ŷ(0)− y‖2
) and choose 0 < γ ≤ min{1/2, 1/4c}. Moreover, we assume the
stepsize α = O ( λ0/n 2 ) . We make the inductive hypothesis as follows for all 0 ≤ s ≤ k
(i) λmin(G(s)) ≥ λ02 ,
(ii) ‖u(s)‖ ≤ 32c √ n
λ0 ‖ŷ(0)− y‖,
(iii) ‖v(s)‖ ≤ 16c √ n
λ0 ‖ŷ(0)− y‖,
(iv) ‖W (s)‖ ≤ 2c √ m,
(v) ‖A(s)‖ ≤ 2c √ m,
(vi) ‖ŷ(s)− y‖2 ≤ (1− αλ0/2)s‖ŷ(0)− y‖2.
Proof. Note that Lemma 2.1 still holds if one chooses a larger c. Thus, we choose c % √ λ0/‖X‖. By using the inductive hypothesis, we have for any 0 ≤ s ≤ k
‖φi(s)‖ = ‖ 1√ m σ(W (s)xi)‖ ≤ 1√ m ‖W (s)‖ ≤ 2c
and
‖Φ(s)‖ ≤ ‖Φ(s)‖F = ( n∑ i=1 ‖φi(s)‖2 )1/2 ≤ 2c √ n. (33)
By using Lemma 2.2, we obtain the upper bound for the equilibrium point zi(s) for any 0 ≤ s ≤ k as follows
‖zi(s)‖ ≤ 1
1− γ0 ‖φi(s)‖ = 2‖φi(s)‖ ≤ 4c,
where the last inequality is because we choose γ0 = 1/2, and
‖Z(s)‖ ≤ ‖Z(s)‖F = ( n∑ i=1 ‖zi(s)‖2 )1/2 = 4c √ n. (34)
By using the upper bound of φi(s), we obtain for any 0 ≤ s ≤ k
‖∇vL(s)‖ ≤ n∑ i=1 |ŷi(s)− yi| ‖φi(s)‖
≤2c n∑ i=1 |ŷi(s)− yi| ≤2c √ n‖ŷ(s)− y‖ ≤2c √ n(1− αλ0/2)s/2‖ŷ(0)− y‖.
Let β , √ 1− αλ0/2. Then the upper bound of ‖∇vL(s)‖ can be written as
‖∇vL(s)‖ ≤ 2c √ nβs‖ŷ(0)− y‖, (35)
and
‖v(k + 1)− v(0)‖ ≤ k∑ s=0 ‖v(s+ 1)− v(s)‖ = α k∑ s=0 ‖∇vL(s)‖
≤α · 2c √ n‖ŷ(0)− y‖ · k∑ s=0 βs
= 2(1− β2) λ0 · 2c √ n‖ŷ(0)− y‖1− β k+1 1− β ≤8c √ n
λ0 ‖ŷ(0)− y‖,
where the last inequality we use the facts β < 1. By triangle inequality, we obtain
‖v(k + 1)‖ ≤ ‖v(k + 1)− v(0)‖+ ‖v(0)‖ ≤ 16c √ n
λ0 ‖ŷ(0)− y‖,
which proves the result (iii). Similarly, we can upper bound the gradient of u
‖∇uL(s)‖ ≤ n∑ i=1 |ŷi(s)− yi| ‖zi‖ ≤ 4c √ n‖ŷ(s)− y‖ ≤ 4c √ nβs‖ŷ(0)− y‖ (36)
so that
‖u(k + 1)− u(0)‖ ≤ 16c √ n
λ0 ‖ŷ − y‖,
and
‖u(k)‖ ≤ ‖u(k)− u(0)‖+ ‖u(0)‖ ≤ 32c √ n
λ0 ‖ŷ(0)− y‖.
The result (ii) is also obtained.
By using the inductive hypothesis, we can upper bound the gradient ofW as follows
‖∇WL(s)‖ ≤ n∑ i=1 1√ m |ŷi(s)− yi| ‖Ei(s)‖ ( ‖Ui(s)−1u(s)‖+ ‖v(s)‖ ) ‖xi‖
≤128c √ n
λ0 √ m ‖ŷ(0)− y‖ n∑ i=1 |ŷi(s)− yi|
≤ 128cn λ0 √ m ‖ŷ(0)− y‖ · ‖ŷ(s)− y‖
≤ 128cn λ0 √ m ‖ŷ(0)− y‖2 · βs, (37)
so that
‖W (k + 1)−W (0)‖ ≤α k∑ s=0 ‖∇WL(s)‖
≤α · 128cn λ20 √ m ‖ŷ(0)− y‖2 · k∑ s=0 βs
≤ 512cn λ20 √ m ‖ŷ(0)− y‖2 ≤ √ mλ0
16c‖X‖2
≤R,
where the third inequality holds is because m is large, i.e., m = Θ ( c2n‖X‖2
λ30 ‖ŷ(0)− y‖2
) . To
simplify the notation, we assume
m = Cc2n‖X‖2
λ30 ‖ŷ(0)− y‖2 (38)
for some large number C > 0. Moreover, we obtain
‖W (k + 1)‖ ≤ ‖W (k + 1)−W (0)‖+ ‖W (0)‖ ≤ 2c √ m,
Therefore, it follows from Lemma 3.4 that λmin{G(k + 1)} ≥ λ02 . Thus, the results (i) and (iv) are established.
By using similar argument, we can upper bound the gradient ofA as follows Note that
‖∇AL(s)‖ ≤ n∑ i=1 γ√ m |ŷi(s)− yi| ‖Di‖‖Ui(s)−1‖‖u(s)‖‖z∗i ‖
≤64c √ n
λ0 √ m ‖ŷ(0)− y‖ · n∑ i=1 |ŷi(s)− yi|
≤ 64cn λ0 √ m ‖ŷ(0)− y‖ · ‖ŷ(s)− y‖
≤ 64cn λ0 √ m ‖ŷ(0)− y‖2 · βs,
so that
‖A(k + 1)−A(0)‖ ≤α k∑ s=0 ‖∇AL(s)‖
≤α · 64cn λ0 √ m ‖ŷ(0)− y‖2 · k∑ s=0 βs
≤ 256cn λ20 √ m ‖ŷ(0)− y‖2.
Since m = Cc 2n‖X‖2 λ30 ‖ŷ(0)− y‖2 and c, C > 0 are large enough, we have
‖A(k + 1)‖ ≤ ‖A(k + 1)−A(0)‖+ ‖A(0)‖ ≤ 2c √ m.
Therefore, the result (v) is obtained and the equilibrium points zi(k + 1) exists for all i ∈ [n]. To establish the result (vi), we need to derive the bounds between equilibrium points Z(k) and feature vectors Φ(k). We firstly bound the difference between equilibrium points zi(k + 1) and zi(k). For any ` ≥ 1, we have
‖z`+1i (k + 1)− z `+1 i (k)‖ =‖σ [ γ̃A(k + 1)z`i (k + 1) + φi(k + 1) ] − σ [ γ̃A(k)z`i (k) + φi(k) ] ‖
≤‖γ̃A(k + 1)z`i (k + 1) + φi(k + 1)− γ̃A(k)z`i (k)− φi(k)‖ ≤γ̃‖A(k + 1)z`i (k + 1)−A(k)z`i (k)‖+ ‖φi(k + 1)− φi(k)‖,
where the first term can be bounded as follows
γ̃‖A(k + 1)z`i (k + 1)−A(k)z`i (k)‖ ≤γ̃‖A(k + 1)−A(k)‖‖z`i (k + 1)‖+ γ̃‖A(k)‖‖z`i (k + 1)− z`i (k)‖ ≤γ̃α‖∇AL(k)‖(4c) + γ̃‖A(k)‖‖z`i (k + 1)− z`i (k)‖
≤64αcn λ0m ‖ŷ(0)− y‖2βk + (1/2)‖z`i (k + 1)− z`i (k)‖,
and the second term is bounded as follows
‖φi(k + 1)− φi(k)‖
= 1√ m ‖σ[W (k + 1)xi]− σ[W (k)xi]‖
≤ 1√ m ‖W (k + 1)−W (k)‖‖xi‖
≤ α√ m ‖∇WL(k)‖
≤128αcn λ0m ‖ŷ(0)− y‖2 · βk.
Thus, we obtain
‖z`+1i (k + 1)− z `+1 i (k)‖ ≤(1/2)‖z ` i (k + 1)− z`i (k)‖+
256αcn
λ0m ‖ŷ(0)− y‖2 · βk
≤(1/2)`‖z1i (k + 1)− z1i (k)‖+ 256αcn
λ0m ‖ŷ(0)− y‖2 · βk · ∞∑ j=0 2−j
≤(1/2)`‖z1i (k + 1)− z1i (k)‖+ 512αcn
λ0m ‖ŷ(0)− y‖2 · βk.
By letting `→∞, we obtain
‖zi(k + 1)− zi(k)‖ ≤ 512αcn
λ0m ‖ŷ(0)− y‖2 · βk.
By using the Cauchy-Schwartz’s inequality, we have
‖Z(k + 1)−Z(k)‖ ≤ ‖Z(k + 1)−Z(k)‖F ≤ 512αcn3/2
λ0m ‖ŷ(0)− y‖2 · βk. (39)
In addition, we will also bound the difference in φi(k + 1) and φi(k). Note that
‖φi(k + 1)− φi(k)‖ = 1√ m ‖σ[W (k + 1)xi]− σ[W (k)xi]‖ ≤ 128αcn λ0m ‖ŷ(0)− y‖2 · βk,
so that
‖Φ(k + 1)−Φ(k)‖ ≤ ‖Φ(k + 1)−Φ(k)‖F ≤ 128αcn3/2
λ0m ‖ŷ(0)− y‖2 · βk (40)
Now, we are ready to establish the result (vi). Note that
‖ŷ(k + 1)− y‖2 =‖ŷ(k + 1)− ŷ(k) + ŷ(k)− y‖2
=‖ŷ(k + 1)− ŷ(k)‖2 + 2 〈ŷ(k + 1)− ŷ(k), ŷ(k)− y〉+ ‖ŷ(k)− y‖2.
In the rest of this proof, we will bound each term in the above inequality. By the prediction rule of ŷ, we can bound the difference between ŷ(k + 1) and ŷ(k) as follows
‖ŷ(k + 1)− ŷ(k)‖ =‖Z(k + 1)u(k + 1) + Φ(k + 1)v(k + 1)−Z(k)u(k)−Φ(k)v(k)‖ ≤‖Z(k + 1)u(k + 1)−Z(k)u(k)‖+ ‖Φ(k + 1)v(k + 1)−Φ(k)v(k)‖,
where the first term can be bounded as follows by using (34), 36, 38, 39, hypothesis (ii), and a large constant C0 > 0
‖Z(k + 1)‖‖u(k + 1)− u(k)‖+ ‖Z(k + 1)−Z(k)‖‖u(k)‖ =α‖Z(k + 1)‖‖∇uL(k)‖+ ‖Z(k + 1)−Z(k)‖‖u(k)‖ ≤αC0c2n‖ŷ(0)− y‖ · βk,
and the second term is bounded as follows by using (33), 35, 40, 38, hypothesis (iii), and a large constant C0 > 0
‖Φ(k + 1)‖‖v(k + 1)− v(k)‖+ ‖Φ(k + 1)−Φ(k)‖‖v(k)‖ =α‖Φ(k + 1)‖‖∇vL(k)‖+ ‖Φ(k + 1)−Φ(k)‖‖v(k)‖ ≤αC0c2n‖ŷ(0)− y‖ · βk.
Therefore, we have
‖ŷ(k + 1)− ŷ(k)‖ ≤ αC0c2n‖ŷ(0)− y‖ · βk, (41)
where the scalar 2 is absorbed in C0 and the constant C0 is difference from C.
Let g , Z(k)u(k + 1) + Φ(k)v(k + 1). Then we have
〈ŷ(k + 1)− ŷ(k), ŷ(k)− y〉 = 〈ŷ(k + 1)− g, ŷ(k)− y〉+ 〈g − ŷ(k), ŷ(k)− y〉 .
Let us bound each term individually. By using the Cauchy-Schwartz inequality, we have
〈ŷ(k + 1)− g, ŷ(k)− y〉 = 〈(Z(k + 1)−Z(k))u(k + 1), ŷ(k)− y〉+ 〈(Φ(k + 1)−Φ(k))v(k + 1), ŷ(k)− y〉 ≤ (‖Z(k + 1)−Z(k)‖‖u(k + 1)‖+ ‖Φ(k + 1)−Φ(k)‖‖v(k + 1)‖) ‖ŷ(k)− y‖ ≤αC0c2n‖ŷ(0)− y‖ · βk‖ŷ(k)− y‖, by (34), 36, 38, 39 ≤αC0c2n · β2k‖ŷ(0)− y‖2. (42)
By using ∇uL(k) = Z(k)T (ŷ(k) − y), ∇vL(k) = Φ(k)T (ŷ(k) − y) and λmin(G(k)) ≥ λ0/2, we get
〈g − ŷ(k), ŷ(k)− y〉 = −α(ŷ(k)− y)T [ Z(k)Z(k)T + Φ(k)Φ(k)T ] (ŷ(k)− y)
≤ −αλ0 2 ‖ŷ(k)− y‖2. (43)
By combining the inequalities (41), 42, 43, we obtain ‖ŷ(k + 1)− y‖2 ≤ ( 1− α [ λ0 − C0c2n− αC20c4n2 ]) β2k‖ŷ(0)− y‖2
≤ (
1− αλ0 2
) · β2k‖ŷ(0)− y‖2
= ( 1− αλ0
2
)k+1 ‖ŷ(0)− y‖2,
where the second inequality is by α = O ( λ0 n2 ) . This proves the result (vi) and complete the proof.

Review summaries:
 In this paper, the authors theoretically analyze the convergence of gradient descent for an implicit neural network with infinite layers with ReLU activation.  The authors show the unique fixed point of the infinite-layered mapping when the weight matrix $\boldsymbol{A}$ has a properly bounded spectral norm.  Using implicit differentiation, the authors show the partial gradient at the fixed point.  Furthermore, the authors show the linear convergence rate by proving the strictly positive-definite of the Gram matrix $\boldsymbol{G}(t)$  (and $\boldsymbol{H}(t)$). 

The paper presents a proof of exponential convergence to global optimality in the over-parametrization settings for an implicit model with scaled weights parameters. Although existing work has established similar proofs for feedforward explicit neural networks, such methods don't work with non-linearly activated implicit models where the well-posedness issue poses challenges to the training process. The authors shows that by scaling the weights, well-posedness can be ensured. The convergence result is obtained first on continuous settings and is then extended to discrete settings. Numerical experiments on real datasets confirms the finding.