On the Theory of Implicit Deep Learning: Global Convergence with Implicit Layers
Authors: Kenji Kawaguchi
Venue: ICLR 2021 Spotlight (2021)
Keywords: Implicit Deep Learning
Deep Equilibrium Models
Gradient Descent
Learning Theory
Non-Convex Optimization

Abstract: A deep equilibrium model uses implicit layers, which are implicitly defined through an equilibrium point of an infinite sequence of computation. It avoids any explicit computation of the infinite sequence by finding an equilibrium point directly via root-finding and by computing gradients via implicit differentiation. In this paper, we analyze the gradient dynamics of deep equilibrium models with nonlinearity only on weight matrices and non-convex objective functions of weights for regression and classification. Despite non-convexity, convergence to global optimum at a linear rate is guaranteed without any assumption on the width of the models, allowing the width to be smaller than the output dimension and the number of data points. Moreover, we prove a relation between the gradient dynamics of the deep implicit layer and the dynamics of trust region Newton method of a shallow explicit layer. This mathematically proven relation along with our numerical observation suggests the importance of understanding implicit bias of implicit layers and an open problem on the topic. Our proofs deal with implicit layers, weight tying and nonlinearity on weights, and differ from those in the related literature.



Published as a conference paper at ICLR 2021
ON THE THEORY OF IMPLICIT DEEP LEARNING: GLOBAL CONVERGENCE WITH IMPLICIT LAYERS
Kenji Kawaguchi Harvard University Cambridge, MA 02138, USA kkawaguchi @ fas. ha ewwrd. edu
ABSTRACT
A deep equilibrium model uses implicit layers, which are implicitly defined through equilibrium point of infinite sequence of computation. It avoids any explicit computation of the infinite sequence by finding an equilibrium point directly via root-finding and by computing gradients via implicit differentiation. In this paper, we analyze the gradient dynamics of deep equilibrium models with nonlinearity only on weight matrices and non-convex objective functions of weights for regression and classification. Despite non-convexity, convergence to global optimum a linear rate is guaranteed without any assumption on the width of the models, allowing the width to be smaller than the output dimension and the number of data points. Moreover, we prove relation between the gradient dy- namics of the deep implicit layer and the dynamics of trust region Newton method shallow explicit layer. This mathematically proven relation along with our numerical observation suggests the importance of understanding implicit bias of implicit layers and an open problem on the topic. Our proof's deal with implicit layers, weight tying and nonlinearity weights, and differ from those in the re- lated literature.
INTRODUCTION
A feedforward deep neural network consists of a stack of H layers, where H is the depth of the network. The value for the depth H is typically hyperparameter and is chosen by network designers (e.g., ResNet-101 in He et al. 2016). Each layer computes some transformation of the output of the previous layer. Surprisingly, several recent studies achieved results competitive with the state-of- the-art performances by using the same transformation for each layer with weight tying (Dabre & Fujita, 2019; Bai et al., 2019b; Dehghani et al., 2019). In general terms, the output of the I-th layer with weight tying can be written by
( = 1122111);; xx, 0) for 1 = 1,2, 11-,,.,,,HH-- 1,
H
(1)
where x is the input to the neural network, 20)) is the output of the 1-th layer (with 20)) x), represents the trainable parameters that are shared among different layers (i.e., weight tying), and 2(1-1) + h22(1-1); () is some continuous function that transforms (1-1) given and 0. With weight tying, the memory requirement does not increase as the depth increases in the forward pass. However, the efficient backward pass compute gradients for training the network usually requires to store the values of the intermediate layers. Accordingly, the overall computational re- quirement typically increases as the finite depth H increases even with weight tying.
Instead of using finite depth H, Bai et al. (2019a) recently introduced the deep equilibriam model that is equivalent to running an infinitely deep feedforward network with weight tying. Instead of running the layer-by-laye computation in equation (1), the deep equilibrium model uses root- finding to directly compute a fixed point =* = lim =!0) where the limit can be ensured to exist by choice of h. We can train the deep equilibrium model with gradient-based optimization by analytically backpropagating through the fixed point using implicit differentiation (e.g.. Griewank & Walther, 2008; Bell & Burke, 2008; Christianson, 1994). With numerical experiments, Bai et al. (2019a) showed that the deep equilibrium model can improve performance over previous state-of- the-art models while significantly reducing memory consumption.

---

Published as a conference paper at ICLR 2021
Despite the remarkable performances of deep equilibrium models, our theoretical understanding of its properties is yet limited. Indeed, immense efforts are still underway to mathematically understand deep linear networks, which have finite values for the depth H without weight tying (Saxe et al., 2014; Kawaguchi, 2016; Hardt & Ma, 2017: Laurent& & Brecht, 2018; Arora et al., 2018; Bartlett et al., 2019; Du & Hu, 2019; Arora et al., 2019a; Zou al., 2020b). In deep linear networks, the funetion h each layer is linear in A and linear in x; i... the map (Ir, 0) + h22(t-1); xz, 8) is bilinear. Despite this linearity, several key properties of deep learning are still present in deep linear networks. For example, the gradient dynamies is nonlinear and the objective function is non- convex. Accordingly, understanding gradient dynamics of deep linear networks is considered to be valuable step towards the mathematical understanding of deep neural networks (Saxe al., 2014; Arora et al., 2018; 2019a).
In this paper, inspired by the previous studies of deep linear networks, we initiate a theoretical study of gradient dynamics of deep equilibrium linear models as a step towards theoretically understand- ing general deep equilibrium models. As we shall see in Section 2, the function h at each layer is nonlinear in 0 for deep equilibrium linear models, whereas it is linear for deep linear networks. This additional nonlinearity is essential to enforce the existence of the fixed point z". The additional non- linearity, the infinite depth, and weight tying are the three key proprieties of deep equilibrium linear models that are absent in deep linear networks. Because of these three differences, we cannot rely on the previous proofs and results in the literature of deep lincar networks. Furthermore, we analyze gradient dynamics, whereas Kawaguchi (2016); Hardt & Ma (2017); Laurent & Brecht (2018) stud- ied the loss landscape of deep linear networks. We also consider general class of loss functions for both regression and classification, whereas Saxe et al. (2014); Arora et al. (2018); Bartlett et al. (2019); Arora et al. (2019a); Zou et al. (2020b) analyzed gradient dynamics of deep linear networks in the setting of the square loss.
Accordingly, we employ different approaches in our analysis and derive qualitatively and quantita- tively different results when compared with previous studies. In Section 2, we provide theoretical and numerical observations that further motivate us to study deep equilibrium linear models. In Scc- tion 3, we mathematically prove convergence of gradient dynamics to global minima and the exact relationship between the gradient dynamics of deep equilibrium linear models and that of the adap- tive trust region method. Section 5 gives a review of related literature, which strengthens the main motivation of this paper along with the above discussion (in Section 1). Finally, Section 6 presents concluding remarks on our results, the limitation of this study, and future research directions.
2 PRELIMINARIES
We begin by defining the notation. We are given a training dataset ((z,, y)){ of n samples where I, E X' C Rm and y; E y S Rm are the i-th input and the i-th target output, respectively. We would like to learn hypothesis (or predictor) from parametric family H = {fo R. - Rm 0 E θ} by minimizing the objective function L (called the empirical loss) over 0 E θ: C())
= {(fe(x:), vs), where 0 is the parameter vector and & R x y + R20 is the loss function that measures the difference between the prediction (axx) and the target y; for each sample. For example, when the parametric family of interest is the class of linear models H. {o + W (o(x)) W E m, Xm], },, the objective function L can be rewritten as:
Lo(W) = Σ ((W
(2)
where the feature map ó is an arbitrary fixed function that is allowed to be nonlinear and is chosen by model designers to transforms an input & E Th= into the desired features Φ(x) E IR.. We use vec(W) E Rm to represent the standard vectorization of a matrix W E Rm,, xm.
Instead of linear models, our interest in this paper lies deep equilibrium models. The output =* of the last hidden layer of a deep equilibrium model is defined by
lim 200 lim h(=((-"); 0) = 1x;";,,,0),
1+30
where the last equality follows from the continuity of 2 + h(2; x. 0) (i.e., the limit commutes with the continuous function). Thus, =* can be computed by solving the equation 2 h(2*; 22, 8) with- out running the infinitely deep layer-by-layer computation. The gradients with respect to parameters are computed analytically via backpropagation through x* using implicit differentiation.
2

---

Published as a conference paper at ICLR 2021
2.1 DEEP EQUILIBRIUM LINEAR MODELS
A deep equilibrium linear model is an instance of the family of deep equilibrium models and is defined by setting the function h at each layer as follows:
h(2 1); x. 0) = (o))
d(x),
(4)
where 0 = (A, B) with two trainable parameter matrices A E Rm and B E R" Xm Along with a positive real number 7 (0, 1), the nonlinear function is used ensure the existence of the fixed point and is defined by (A4)ss = ET_, (xp(A23))
exp(A,,) The class of deep equilibrium linear
models is given by H. = {o + B (limy-+oo (z, ()) A E RmX" B E R=',** where (x, A) = (A)= + Φ(x). Therefore, the objective function for deep equilibrium linear models can be written as
I(A = Σe 1=1 [1000
e B
BBling lim (12E .b6
(5)
The outputs of deep equilibrium linear models fo(x) B (limy 20 (x, A)) are nonlinear and non-multilinear in the optimization variable A. This is in contrast to linear models and deep linear networks. From the optimization viewpoint, linear models W WΦ(x;) are called linear because they linear in the optimization variables W. Deep linear networks W(H) .... W()), x are multi- linear in the optimization variables (W()), W/ ..... W(")) (this holds also when replace by 6(x)). This difference creates a challenge in the analysis of deep equilibrium linear models.
Following previous works gradient dynamics of different machine learning models (Saxe et al., 2014; Ji & Telgarsky, 2020), we consider the process learning deep equilibrium linear models via gradient flow: OL (A., a
AL aB
di 4. - DA B.),
- B. di.
Ar, B,), A1 2 0,
(6)
where (A,, B,) represents the model parameters at time with an arbitrary initialization (Ag, Bo). Throughout this paper, feature map é and real number y (0, 1) are given and arbitrary (except in experimental observations) and we omit their universal quantifiers for the purpose of brevity.
2.2 PRELIMINARY OBSERVATION FOR ADDITIONAL MOTIVATION
Our analysis chiefly motivated as a step towards mathematically understanding general deep equi- librium models (as discussed in Sections I and 5). In addition to the main motivation, this section provides supplementary motivations through theoretical and numerical preliminary observations.
In general deep equilibrium models, the limit, lim; 209), is not ensured exist (see Appendix C). In this view, the class of deep equilibrium linear models is one instance where the limit is guaranteed to exist for any values of model parameters as stated in Proposition 1:
Proposition 1. Given any (z, A). the sequence (21) (x,A)); in Euclidean space RT converges. Proof. We use the nonlinearity a to ensure the convergence in our proof in Appendix A.5.
Proposition 1 shows that we can indeed define the deep equilibrium linear model with lim 20) =z. A). Therefore, understanding this model is a sensible starting point for theory of general deep equilibrium models.
As our analysis has been mainly motivated for theory, it would be of additional value to discuss whether the model would also make sense in practice, at least potentially the future. Consider an (unknown) underling data distribution P(x, y) = (()r)))(((r) Intuitively, if the mean of the P(y(=) approximately given by a (true unknown) deep equilibrium linear model, then would make sense to use the parametric family of deep equilibrium linear models to have the inductive bias in practice. To confirm this intuition, we conducted numerical simulations. To generate datasets, we first drew uniformly random 200 input images for input data points 2; from standard image dataset CIFAR-10, CIFAR-100 or Kuzushiji-MNIST (Krizhevsky & Hinton, 2009; Clanuwat et al., 2019).
We then generated targets as y: = B" (limy 202) (x; 4²)) + 5, where 6, 1id N/ 1). Each entry of the true (unknown) matrices A and B" was independently drawn from the standard normal distribution. For each dataset generated in this way, we used stochastic gradient descent (SGD) to train linear models, fully-connected feedforward deep neural networks with ReLU nonlinearity
3

---

Published as a conference paper at ICLR 2021
Linear (best)
550
Linear (werst)
DELM (best)
DELM
Enst
DELM
a6o 500
3000
vie
ares
5900
200
4500
5000
3900
200
wid
500
(a) Modified Kuzushij--MNIS) Linear v.s. DELM
(b) Modified Kuzushiji-MNiSY DNN v.s. DELM
Linear (best)
UNN (He
(warst)
DNN (11-3)
DELM (bwel)
DELM (worst)
DEL (heet)
DNN (H-4]
DELM (wers))
torth
1991
(c) Modified CIFAR-10: Linear v.s. DELM
(d) Modified CIFAR-10: DNN v.s. DELM
Linear (best)
Linwar (worst)
CoR tos
DELM (best)
DELM (worst)
DEL (beet)
DELM twars)
1000
3200
0003
3000
1006
2003 2003
0003
1200
1000
2000
2003
4200
3000
3000
2006
2006
40:3
1000
(e) Modified CIFAR-100: Linear v.s. DELM
(1) Modified CIFAR-100: DNN v.s. DELM
Figure I: Preliminary observations for additional motivation to theoretically understand deep equi- librium linear models. The figure shows test and train losses versus the number of epochs for linear models, deep equilibrium linear models (DELMs), and deep neural networks with ReLU (DNNs).
(DNNs), and deep equilibrium linear models. For all models, we fixed Φ(x) = x. See Appendix D for more details of the experimental settings.
The results of this numerical test are presented in Figure 1. In the figure, the plotted lines indicate the mean values over five random trials whereas the shaded regions represent error bars with one standard deviation. The plots for linear models and deep equilibrium linear models are shown with the best and worst learning rates (separately for each model in terms of the final test errors at epoch = 5000) from the set of learning rates SiR= = {0.01, 0.005, 0.001, 0.0005, 0.0001, 0.00005). The plots for DNNs are shown with the best learning rates (separately for each depth H) from the set SLR- As can be seen, all models preformed approximately the same at initial points, but deep equilibrium linear models outperformed both linear models and DNNs in test eтrors after training, confirming our intuition above. Moreover, we confirmed qualitatively same behaviors with four more datasets as well as for DNNs with and without bias terms in Appendix D. These observations additionally motivated us to study deep equilibrium linear models to obtain our main results in the next section. The purpose of these experiments to provide a secondary motivation for our theoretical analyses.
MAIN RESULTS
In this section, we establish mathematical properties of gradient dynamics for deep equilibrium linear models by directly analyzing its trajectories. We prove linear convergence to global minimum in Section 3.1 and further analyze the dynamics from the viewpoint of trust region in Section 3.2.
3.1 CONVERGENCE ANALYSIS
We begin in Section 3.1.1 with a presentation of the concept of the Lojyak--ojsiewiizzzz (PL) inequal- ity and additional notation. The PL inequality is used to regularize the choice of the loss functions f in our main convergence theorem for a general class of losses in Section 3.1.2. We conclude in Section 3.1.3 by providing concrete examples of the convergence theorem with the square loss and the logistic loss, where the PL inequality is no longer required as the PL incquality is proven to be satisfied by these loss functions.
3.1.1 THE
INEQUALITY AND ADDITIONAL NOTATION
In our context, the notion of the PI. inequality is formally defined as follows:
4

---

Published as a conference paper at ICLR 2021
Definition 1. The function Ly is said to satisfy the Polyak-Lojasies (PL) inequality with radius
(1))) R E (0, oo] and parameter k > 0 if 2 (vec(()) 2 x (vec(W)) - L for all [[W]] < R, where L.0 (vec(-)) == (o(+) and L == inf [ [ (W)).
With any radius R > 0 sufficiently large (such that it covers the domain of Lg), Definition 1 becomes equivalent to the definition of the PL inequality in the optimization literature (e.g., Polyak. 1963; Karimi et al., 2016). See Appendix C for additional explanations on the equivalence. In general, the non-convex objective function L of deep equilibrium linear models does not satisfy the PL inequality. Therefore, we cannot assume the inequality on L. However, in order obtain linear convergence for a general class of the loss functions L, we need some assumption on f: otherwise, we can choose a loss to violate the convergence. Accordingly, we will regularize the choice of the loss & through the PL. inequality on the function Lg W + Ei= ((Wo(z;) U:).
The PL inequality with radius R. E (0, oo] (Definition 1) leads the notion of the global minimum value in the domain corresponding to the radius in our analysis: L inl A AERmXm BBBB L(A, B), where Bn { E Rm | 11 < (1 R). With R oc, this recovers the global minimum value L in the unconstrained domain as Li L' == inf AgRm _gR L(A, B). Furthermore, if a global minimum AA,,B) E R= x Rm exists, there exists R < oo such that for any R E: [R. oo), we have 8 E BR and thus Li = L In other words, if a global minimum exists, using a (sufficiently large) finite radius < suffices to obtain L L,".
We close this subsection by introducing additional notation. For a real symmetric matrix M, we use (M) to represent its smallest eigenvalue. For an arbitrary matrix M E Rdxdd', we let rank(Mf) be its rank, ||p be its matrix norm induced by the vector p-norm, miin (M) be its smallest singular value (i.e., the min(d, (d, d' >-th largest singular value), M.j be its j-th column vector in Rd, and M,. be its i-th row vector in IRd", For d E Ns-o, we denote by la the identify matrix in Rdxd We define the Jacobian matrix JR,t E Rm xm of the vector-valued function A + oA)) sk
by (Jk,))s DAje A for all 2 0 and k = I,..., m. Finally, we define the feature matrix E RmxX by tx = ((T)) for ls = 1, 1,...,m, m and i 1, I....,n.
(A)ak)
3.1.2 MAIN CONVERGENCE THEOREM
Using the PL inequality only on the loss function & through Lo (Definition 1), we present our main theorem a guarantee on linear convergence to global minimum for the gradient dynamics of the non-convex objective L for deep equilibrium linear models:
Theorem 1. Let 8 my x y + R2 be arbitrary such that the function q - Elg, We) is differ- entiable for any E [l....,n] n) (with an arbitrary my E Ny50 and an arbitrary y). Then, for any
T > 0, R (0, oo] and > 0 such that [Bllll < (1-7)R for all t E [0,T] and Lo satisfies the
PL inequality with the radius R and the parameter K, the following holds:
L(AT, BT) S Liz + (L(Ao, Bo) Lā,R) -2nd7T
where dy infre Amin (D.) and D, a positive definite matrix defined by
D. : Σ (o, )ax @ (Im, +
JU,JE, TB.)
(8)
with U,: Im (Ai). Furthermore, XT Miitt for any T 2 (limy= Av m1+))
Proof. The additional nonlinearity a creates complex interaction among hidden neurons. This interaction is difficult to be factorized out for the gradients of L with respect to A. This different from but analogous to the challenge deal with nonlinear activations in the loss landscape of (non- overparameterizzed) deep nonlincar networks, for which previous works have made assumptions of sparse connections to factorize the interaction (Kawaguchi et al., 2019). In contrast, we do not rely on sparse connections. Instead, we observe that although it is difficult to factorize this complex interaction (due to the nonlinearity σ) in the space of loss landscape, we can factorize it in the space of gradient dynamics. See Appendix A. for the proof overview and the complete proof.
Theorem shows that in the worst case for Xr, the optimality gap decreases exponentially towards 2x zero I(AA, Br) T S Coe mitt57? where Co I(Ao,, Bo) L Therefore, for any
5

---

Published as a conference paper at ICLR 2021
desired accuracy > 0, setting Coe
mi++5)²
T
S c and solving for T yield that L(Ao, Bo) Lo.
m(1 + 7)2 log
(9) Theorem also states that the rate of convergence improves further depending on the quality of the matrix Dj (defined in equation (8) in terms of its smallest eigenvalue
L(Ar,Br) Se
for any T
2x
the partieular trajectory (Ar, B,) up to the specific time S T; i.c.. Ar infre(o.2) Amin (Dr). This opens up the direction of future work for further improvement of the convergence rate through the design of initialization (Ag, BB) to maximize Xr for trajectories generated from a specific initialization scheme.
3.1.3 EXAMPLES: SQUARE LOss AND LOGISTIC Loss
The main convergence theorem in the previous subsection is stated for any radius R E oo] and parameter & > 0 that satisfy the conditions on [Bella and the PL inequality (see Theorem 1). The values of these variables are not completely specified there they depend on the choice of the loss functions &. In this subsection, we show that these values can be specified further and the condition on PL inequality can be discarded by considering specific choice of loss functions f.
In particular, by using the square loss for l, we prove that we can set R = do and K = 20 ( (D) Corollary 1. Let E(g,y) my E Nyo). Assume that rank ④) = min(n, m.. Then for any T > 0,
||g yı 11 where yi E R for
1.2,...,M n (with an arbitrary
L(AT. Br) L + (L(Ao: Bo) - L 4a (0)" XT
where o muin (+) Lg = inf w Km LWW,, and Xy inf Amiu (Dt) 2 m(1+1)²
Proof. This statement follows from Theorem 1. The conditions on || Be and the PL inequality (in Theorem 1) are now discarded by using the property of the square loss f. See Appendix A.3 for the complete proof.
In Corollary I, the global linear convergence is established for the square loss without the notion of the radius R as we set R = co. Even with the square loss, the objective function L is non-convex. Despite the con-ooveex,y, Corollary shows that for any desired accuracy >
e for any T> m(1 + 7) L(Au, Bo) Lg
L(AT, BT) L Se
log
(10)
Corollary allows both cases of m S n and m > n. In the case of over peee-aarametrizzaioon m > n, the covariance matrix ΦΦ E mm (or XX with ((x) = x) is always rank deficient because rank( = rank (b) < n < m. This implies that the Hessian of Lo is always rank deficient, because the Hessian of Lo is V (vee"(ve)(W))
2[ 2 Imy] E mymmmmmm (see Appendix A.3 for its derivation) and because rank([04 @ Im,,)) rank T) rank(Tmy) S my, < my,. Since the strong convexity on a twice differentiable function requires its Hessian to be of full rank, this means that the objective Ly for linear models is not strongly convex in the case of over- parameterization m > n. Nevertheless, we establish the linear convergence to global minimum for deep equilibrium linear models in Corollary for both cases of m > n and m S n by using Theorem 1.
For the logistic loss for f, the following corollary proves the global convergence at a linear rate: Corollary 2. Let E(g, Rk) = y/ log( Ite-q) - (1 y:) log(1 ( + Tlgll1 with an arbitrary
Then for any T > 0 and R E (0, oo) such that |B < (1 - ( Ror all E T), the following holds:
0 where yi {0, 1} for i 1.2..... n. Assume that rank(
L(AT, Br) S 1 + (I.(Ao: Bo) - T.Ō.R) e
where muin
( (Φ) > Ar == inf inf Amin m(1 and
( 1
)
p(R) = inf WR, SE ())
1 1 e Wo((a)
1
1+ Wo(z))
20 0.
Proof. This statement follows from Theorem 1 by proving that the condition on PL inequality is satisfied with the parameter x = (2T + p(R))ominn ( See Appendix A.4 for the complete proof.

---

Published as a conference paper at ICLR 2021
In Corollary 2, we can also set R remove the notion of the radius R from the statement of the global convergence for the logistic loss. By setting R = co, Corollary 2 states that for any T > 0,
L(AT, Br) S L + (L(Ao, Bo) Lg) (rammnn()))
for the logistic loss. For any T > 0, this implies that for any desired accuracy c <> 0,
for any T> mm1 + 7)2 log L(Ao, Bo) Lö
L(Ar. BT) I
(11)
4ro ($))
In practice, may want to set >> 0 to regularize the parameters (for generalization) and to ensure the existence of global minima (for optimization and identifiability). That is, if we set r = 0 instead, the global minima may not exist in any bounded space, due to the property of the logistic loss. This is consistent with Corollary in that if T = 0, equation (11) does not hold and we must consider the convergence the global minimum value ir defined in bounded domain with a radius R < co. In the case of r 0 and R < co, Corollary 2 implies that for desired accuracy c > O,
for any T ml + 7)2
L(Au, Ba) Lo.R
L(Ar, Br) L. S
2p( (R)) ($)2 log
(12)
where we have R(R) > O because R < co. Therefore, Corollary 2 establish the linear convergence global minimum with both cases of 7 > 0 and 7 0 for the logistic loss.
3.2 UNDERSTANDING DYNAMICS THROUGH TRUST REGION NEWTON METHOD
In this subsection, we analyze the dynamics of deep equilibrium linear models the space of the hypothesis, fa.
B, (oo 20) (I, A,)). For any functions g and 9 with domain X' S R' we write g = g if g(z) g(z) for all x E X.
The following theorem shows that the dynamics of deep equilibrium linear models fe, can be written as fe.o. = + V where J. is scalar and Vn follows the dynamics of a trust region Newton method of shallow models with the (non-standard) adaptive trust region V. This suggests potential benefits of deep equilibrium linear models in two aspects: when compared to shallow models, it can sometimes accelerate optimization via the effect of the implicit trust region method (but not necessarily as the trust region method does not necessarily accelerate optimization) and induces novel implicit bias for generalization via the non-standard implicit trust region V4-
Theorem 2.
Let RT x y + Rzo be arbitrary such that the function q + E(q, y:) differen- tiable for any i {1, [l,...,m}] n) with my = 1 and (an arbitrary Y). Then for any time t 2 0, there exist real number 8, > 0 such that for any d, E (0, dc).
d fa.
1 ViO,
Vid,
J.
voc(VV) E argmii ()),, "
(13)
dt
where V1 = for Rm, |lv|aa, s 6ell vee(B,U, ).., : == U, (S, 5.Fi) U.T > 0, and
T.((w)= == I. (veo (BUU- + (vec (vec(B.U, 1) u
2" T V2 I ( (B,U, 1)v.
Here, F. Em Volalfn (Ix)) (lime 20 (x;,,, Ar)) (lim 200 (zi, A+)) with k((g) Egg y1) and S. : Im + 2² diag(1) of) with of E R and (vf)x- == 11 (B.U,-) Vk.
Proof. This is proven with the Karush-Kuhn-Tun mization problem: minimizovev, (v) See Appendix A.2.
(KKT) conditions for the constrained opti-
When many global minima exist, difference in the gradient dynamics can lead to a significant discrepancy in the learned models: i.c., two different gradient dynamics can find significantly differ- ent global minima with different behaviors for generalization and test accuracies (Kawaguchi et al., 2017). In machine learning, this is an important phenomenon called implicit bias - inductive bias induced implicitly through gradient dynamics and is the subject of emerging active research area (Gunasekar et al.. 2017: Soudry et al., 2018; Gunasekar et al., 2018; Woodworth et al., 2020; Moroshko et al., 2020).
7

---

Published as a conference paper at ICLR 2021
As can be seen in Theorem 2, the gradient dynamics of deep equilibrium linear models fa, differs from that of linear models Wi with any adaptive learning rates, fixed preconditioners, and existing variants of Newton methods. This is consistent with our experiments in Section 2.2 and Appendix D where the dynamics of deep equilibrium linear models resulted in the learned predictors with higher test accuracies, when compared to linear models with any learning rates. In this regard, Theorem 2 provides a partial explanation (and a starting point of the theory) for the observed generaliza- tion behaviors, whereas Theorem (with Corollaries and 2) provides the theory for the global convergence observed in the experiments.
Theorem 2, along with our experimental results. suggests the importance of theoretically under- standing implicit bias of the dynamics with the time-dependent trust region. In Appendix B, we show that Theorem suggests a new type of implicit bias towards a simple function as a result of infinite depth, whereas understanding this implicit bias in more details is left as an open problem for future work.
4 EXPERIMENTS
In this section, we conduct experiments to further verify and demonstrate our theory. To compare with the previous findings, we use the same synthetic data that in the previous work (Zou et al., 2020b); i.e., we randomly generate 2; E R 10 from the standard normal distribution and set yi = -2; + 0.14, for all i E (1.2...., n with standard normal distribution. We set d(x) = x and use the square loss E(g. y.) = |1q yi 1i1l2- As in the previous work, 2020b) and report the results in Figure 2 (a). As can be seen in the figure, deep equilibrium linear models converges to the global minimum value with all initialization and random trials, whereas linear ResNet converges to suboptimal value with identity initialization. This is consistent with our theory for deep equilibrium linear models and the previous work for ResNet (Zou et al., 2020b).
1000, where s is independently generated by the
consider random initialization and identity initialization (Zou et al.,
We repeated the same experiment by generating (I:), independently from the uniform distribution of the interval [-1, 1] instead for all i E {l....,,, n} and k E {1..... m² with n 1000 and m = 10. Figure 2 (b) shows the results of this experiment with the uniform distribution and confirm the global convergence of deep equilibrium linear models again with all initialization and random trials. In this case, linear ResNet with identity initialization also converged to the global minimum value. These observations are consistent with Corollary where deep equilibrium linear models are guaranteed to converge to the global minimum value without any condition the initialization.
We now consider the rate of the global convergence. In Corollary 1, we can set by- #
to get guarantee for the global linear convergence rate for all initializations in theory. However, in practice, this is pessimistic convergence rate and we may want to choose AT depending on initialization. To demonstrate this, using the same data as that in Figure 2 (a), Figure 2 (c) reports the numerical training trajectory along with theoretical upper bounds with Xy m(1+y)2 and initialization-dere dependent \y = infre(o.7) Amiu (Dz). As can be seen in Figure 2 (c), the theoretical bound with ditaliizationdddeendeee Xr demonstrates faster and more accurate convergence rate. A qualitatively same observation is reported for the logistic loss in Appendix D.2.
cetiiiatioonnndeeendde
DELM: idenkity
DELM: identity
trinning trajectory bound Ind-iid boundi optima
DELM: random
random
train loss
DELM: random 2
train loss
DELM: random 2
train loss
DELM: random 3
DELM: random 3
ResNet
global optima
global optima
iterations
of iteir ationns
of therations
(a) Gaussian data
(b) Uniform data
(c) Theoretical bounds
Figure 2: (a)- Convergence performances for deep equilibrium linear models (DELMs) with identity initialization and random initialization of three random trials, and linear ResNet with iden- tity initialization. (c) the numerical training trajectory of DELMs with random initialization along with theoretical upper bounds with pnttialiaatonninddeeeeeen X and imitialization-depende Ay.
8

---

Published as a conference paper at ICLR 2021
5 RELATED WORK
The theoretical study of gradient dynamics of deep networks with some linearized component is highly active area of research. Recently, Bartlett et al. (2019); Du & Hu (2019); Arora et al. (2019a); Zou et al. (2020b) analyzed gradient dynamics of deep linear networks and proved global convergence rates for the square loss under certain assumptions on the dataset, initialization, and network structures. For example, the dataset is assumed to be whitened (i.e., d = 4m or X Im_) and the initial loss is assumed to be smaller than the loss of any rank-deficient solution by Arora et al. (2019a): the input and output layers are assumed to represent special transformations and are fixed during training by Zou et al. (2020b).
Deep networks are also linearized implicitly in the neural tangent kernel (NTK) regime with signif- icant wwee-parameteriaation m > n (Yehudai & Shamir, 2019; Lee et al., 2019). By significantly increasing model parameters (or more concretely the width m), we can ensure deep features or cor- responding NTK to stay nearly the same during training. In other words, deep networks in this regime are approximately linear models with random features corresponding to the NTK at random initialization. Because of this implicit linearization, deep networks in the NTK regime are shown to achieve globally minimum training errors by interpolating all training data points (Zou et al., 2020a; Li & Liang. 2018; Jacot et al., 2018; Du et al., 2019; 2018; Chizat et al.. 2019; Arora et al., 2019b; Allen-Zhu et al., 2019; Lee al., 2019; Fang et al., 2020; Montanari & Zhong. 2020).
These previous studies have significantly advanced our theoretical understanding of deep learn- ing through the study of deep linear networks and implicitly linearized deep networks in the NTK regime. In this context, this paper is expected to contribute to the theoretical advancement through the study of a new and significantly different type of deep models deep equilibrium linear models. In deep equilibrium linear models, the function at each layer A + (x; xx, 0) is nonlinear due to the additional nonlinearity σ: A + h(2(111); 0) == 7o (A)2 linear networks, the function at each layer W + h (1 x. W (())((11)) is linear (it is linear also with skip connection). Furthermore, the nonlincarity o is not an element-wise func- tion, which poses an additional challenge in the mathematical analysis of deep equilibrium linear models. The nonlinearity a, the infinite depth, and weight tying in deep equilibrium linear models necessitated us to develop new approaches in our proofs. The differences in the models and proof' naturally led to qualitatively and quantitatively different results. For example, we do not require any of parametrriztaiiin m > n, interpolation of all training data points, and any assumptions mentioned above for deep linear networks.
((x). In contrast, for deep
Unlike previous papers, we also related the dynamics of deep equilibrium linear models to that of a trust region Newton method of shallow models with G,-quadratic norm. This suggested potential benefits of deep equilibrium linear models. Our theory consistent with our numerical observations.
6 CONCLUSION
For deep equilibrium linear models, despite the non-convexity, we have rigorously proven conver- gence of gradient dynamics to global minima, at a linear rate, for a general class of loss functions, including the square loss and logistic loss. Moreover, we have proven the relationship between the gradient dynamics of deep equilibrium linear models and that of the adaptive trust region method. These results apply models with any configuration on the width of hidden layers, the number of data points, and input/output dimensions, allowing rank-deficient covariance matrices as well as both paaaeettrratttiizatooonn and over-pavane puee-paramttrizztionn.
The crucial assumption for our analysis is the differentiabilityy of the function q + Elg. V:). which is satisfied by standard loss functions, such as the square loss, the logistic loss, and the smoothed hinge loss E(g. Dt) = (max{(0) 1 ):4))* with A > 2. However, it is not satisfied by the (non-smoothed) hinge loss E(g. y1) max (0, 1 y:g), the treatment of which is left to future work. Future work also includes corresponding theoretical analyses with stochastic gradient descent.
Our theoretical results (in Section 3) and numerical observations (in Section 2.2 and Appendix D) uncover the special properties of deep equilibrium linear models, providing basis of future work for theoretical studies of implicit bias and for further empirical investigations of deep equilibrium models. In our proofs, the treatments of the additional nonlinearity a,, the infinite depth, and weight tying are especially unique, and we expect our new proof techniques to be proven useful in further studies of gradient dynamics for deep models.
9

---

Published as a conference paper at ICLR 2021
REFERENCES
J Harold Ahlberg and Edwin N Nilson. Convergence properties of the spline fit. Journal of the Society for Industrial and Applied Mathematics, (1):955104, 1963.
Zeyuan Allen-Zhu, Yuanzhi Li, and Yingyu Liang. Learning and generalization in overparameter- ized neural networks, going beyond two layers. Advances neural information processing systems, pp. 6158-6169, 2019.
Sanjeev Arora, Nadav Cohen, and Elad Hazan. On the optimization of deep networks: Implicit acceleration by overparamcterization In International Conference on Machine Learning, 2018.
Sanjeev Arora, Nadav Cohen, Noah Golowich, and Wei Hu. A convergence analysis of gradient de- scent for deep linear neural networks. In International Conference on Learning Representations, 2019a.
Sanjeev Arora, Simon Du, Wei Hu, Zhiyuan Li, and Ruosong Wang. Fine-grained analysis of optimization and generalization for overparameterized two-layer neural networks. arXiv preprint arXiv:1901.08584 2019b.
Shaojie Bai, J Zico Kolter, and Vladlen Koltun. Deep equilibrium models. Advances in Neural Information Processing Systems, pp. 690-701, 2019a.
Shaojie Bai, J Zico Kolter, and Vladlen Koltun. Trellis networks for sequence modeling. In Inter- national Conference on Learning Representations, 2019b.
Randal J Barnes. Matrix differentiation. Springs Journal, pp- 1-9, 2006.
Peter L Bartlett, David P Helmbold, and Philip M Long. Gradient descent with identity initializa- tion efficiently learns positive-definite lincar transformations by deep residual networks. Neural computation, 31(3):477-502, 2019.
Bradley M Bell and James V Burke. Algorithmic differentiation of implicit functions and optimal values. In Advances in Automatic Differentiation, pp. 67-77. Springer, 2008
Lenaic Chizat, Edouard Oyallon, and Francis Bach. On lazy training in differentiable programming. In Advances in Neural Information Processing Systems, pp. 2937-2947, 2019.
Bruce Christianson. Reverse accumulation and attractive fixed points. Optimization Methods and Software, 3(4):311-326 1994.
Tarin Clanuwat, Mikel Bober-Irizar, Asanobu Kitamoto, Alex Lamb, Kazuaki Yamamoto, and David Ha. Deep learning for classical japanese literature. NeurlPS Creativity Workshop 2019, 2019.
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, pp. 6292- 2019.
Mostafa Dehghani, Stephan Gouws, Oriol Vinyals, Jakob Uszkoreit, and Lukasz Kaiser. Universal transformers. In International Conference on Learning Representations, 2019.
Simon Du and Wei Hu. Width provably matters in optimization for deep linear neural networks. In International Conference on Machine Learning, pp. 1655-1664. 2019.
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, pp. 1675- 1685, 2019.
Simon Du, Jason D Lee, Haochuan Li, Liwei Wang, and Xiyu Zhai. Gradient descent finds global minima of deep neural networks. arXiv preprint arXiv: /8/ 1.080044 2018.
Cong Fang, Jason D Lee, Pengkun Yang, and Tong Zhang. Modeling from features: a mean-fide field framework for pereptarmeteri deep neural networks. arXiv preprint 2007:20077.14522. 2020.
10

---

Published as a conference paper at ICLR 2021
Georg Frobenius. Über matrizen nicht negativen elementen. Sirzungsberichtt der Königlich Preussischen Akademie der Wissenschaften, pp. 456 456-477, 1912.
Andreas Griewank and Andrea Walther. Evaluating derivatives: principles and techniques of algo- rithmic differentiation. SIAM, 2008.
Suriya Gunasekar, Blake E Woodworth, Srinadh Bhojanapalli, Behnam Neyshabur, and Nati Srebro. Implicit regularization in matrix factorization. In Advances in Neural Information Processing Systems, pp. 6151-6159, 2017.
Suriya Gunasekar, Jason Lee, Daniel Soudry, and Nati Srebro. Implicit bias of gradient descent linear convolutional networks. In Advances in Neural Information Processing Systems, pp- 9461-9471, 2018.
Moritz Hardt and Tengyu Ma. Identity matters in deep learning. In International Conference on Learning Representations, 2017.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recog- nition. In Proceedings of the IEEE conference computer vision and pattern recognition, pp- 770-778, 2016.
Arthur Jacot, Franck Gabriel, and Clément Hongler. Neural tangent kernel: Convergence and gen- cralization in neural networks. Advances in neural information processing systems, pp. 8571- 8580, 2018.
Ameya D Jagtap, Kenji Kawaguchi, and George Em Karniadakis. Locally adaptive activation func- tions with slope recovery for deep and physics-infornet neural networks. Proceedings of the Royal Society A, 476(2239):20203344 2020a.
Ameya D Jagtap, Kenji Kawaguchi, and George Em Karniadakis. Adaptive activation functions accelerate convergence in deep and physics-informed neural networks. Journal of Computational Physics, 404: 109136. 2020b.
Ziwei Ji and Matus Telgarsky. Directional convergence and alignment in deep learning. arXiv preprint urXiv:2006.06657 2020.
Hamed Karimi, Julie Nutini, and Mark Schmidt. Linear convergence of gradient and proximal- gradient methods under the polyak-tojasiewin condition. Joint European Conference on Ma- chine Learning and Knowledge Discovery in Databases, pp. 795-811. Springer, 2016.
Kenji Kawaguchi. Deep learning without poor local minima. In Advances in Neural Information Processing Systems, pp. 586, 2016.
Kenji Kawaguchi and Yoshua Bengio. Depth with nonlincarity creates bad local minima in resnets. Neural Nenworks, 11::167-174, 2019.
Kenji Kawaguchi and Jiaoyang Huang. Gradient descent finds global minima for generalizable deep neural networks of practical sizes. In 2019 57th Annual Allerton Conference Communication, Control, and Computing (Allerton), pp. 92-99. IEEE, 2019.
Kenji Kawaguchi and Leslie Kaelbling. Elimination of all bad local minima in deep learning. In International Conference on Artificial Intelligence and Statistics, pp. 853-863. PMLR, 2020.
Kenji Kawaguchi, Leslie Pack Kaelbling, and Yoshua Bengio. Generalization in deep learning. arXiv preprint arXiv: 1w1..::710005466 2017.
Kenji Kawaguchi, Jiaoyang Huang, and Leslie Pack Kaelbling. Effect of depth and width local minima in deep learning. Neural computation, 31 31(7):1462-1499) 2019.
Alex Krizhevsky and Geoffrey Hinton. Learning multiple layers of features from tiny images. Tech- nical report, Citeseer, 2009.
Thomas Laurent and James Brecht. Deep linear networks with arbitrary loss: All local minima are global. In International conference machine learning, pp. 2902-2907. PMLR, 2018.
11

---

Published as a conference paper at ICLR 2021
Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11):2278-233 2324,, 1998.
Jachoon Lee, Lechao Xiao, Samuel Schoenholz, Yasaman Bahri, Roman Novak, Jascha Sohl- Dickstein, and Jeffrey Pennington. Wide neural networks of any depth evolve linear models under gradient descent. In Advances in neural information processing systems, pp. 8572-8583, 2019.
Yuanzhi Li and Yingyu Liang. Learning overparameterized neural networks via stochastic gradient descent on structured data. In Advances in Neural Information Processing Systems, pp. 8157- 8166, 2018.
Shiyu Liang, Ruoyu Sun, Jason D Lee, and R Srikant. Adding neuron eliminate all bad local minima. In Advances in Neural Information Processing Systems, 2018.
Andrea Montanari and Yiqiao Zhong. The interpolation phase transition in neural networks: mem- orization and generalization under lazy training. preprint rxXiv2200.112220 2020.
Nenad Morača. Bounds for norms of the matrix inverse and the smallest singular value. Linear algebra and its applications, 429(10))258) 002599922011. 2008.
Edward Moroshko, Suriya Gunasekar, Blake Woodworth, Jason D Lee, Nathan Srebro, and Daniel Soudry. Implicit bias in deep linear classification: Initialization scale vs training accuracy. arXiv preprint arXiv:2007.067: 06738, 2020.
Yuval Netzer, Tao Wang, Adam Coates, Alessandro Bissacco, Bo Wu, and Andrew Ng. Reading digits in natural images with unsupervised feature learning. In NIPS workshop on deep learning and unsupervised feature learning, 2011.
Quynh Nguyen. On connected sublevel sets in deep learning. In International Conference on Machine Learning, pp. 4790-4799, PMLR, 2019.
Quynh Nguyen. A note connectivity of sublevel sets in deep learning. arXiv:2101.08576, 2021.
arXiv preprint
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, high- performance deep learning library. Advances in neural information processing systems, pp- 8026-8037, 2019.
Oskar Perron. Zur theorie der matrices. Mathematische Annalen, 64(2):248-26) 1907.
Tomaso Poggio, Kenji Kawaguchi, Qianli Liao, Brando Miranda, Lorenzo Rosasco, Xavier Boix, Jack Hidary, and Hrushikesh Mhaskar. Theory deep learning iii: explaining the non-overfitiing puzzle. arXiv preprint (tti://800000170 2017.
Boris Teodorovich Polyak. Gradient methods for minimizing functionals. Zhurnal Vychislitel noi Matematiki Matematicheskoi Fiziki, 3(4):643-653, 1963.
Andrew M Saxe, James L McClelland, and Surya Ganguli. Exact solutions to the nonlinear dy- namics of learning in deep linear neural networks. In International Conference Learning Representations, 2014.
Ohad Shamir. Are ResNets provably better than linear predictors? Advances in Neural Informa- fion Processing Systems, appear, 2018.
Daniel Soudry, Elad Hoffer, Mor Shpigel Nacson, Suriya Gunasekar, and Nathan Srebro. The im- plicit bias gradient descent on separable data. The Journal of Machine Learning Research, 19 (1):2822-2878, 2018.
Tactile Srl and Italy Brescia. Semeion handwritten digit data set. Semeion Research Center of Sciences of Communication, Rome, Italy, 1994.
James Varah. A lower bound for the smallest singular value of a matrix. Linear Algebra and its Applications, 11(1):3-5, 1975.
12

---

Published as a conference paper at ICLR 2021
Vikas Verma, Meng Qu, Kenji Kawaguchi, Alex Lamb, Yoshua Bengio, Juho Kannala, and Jian Tang. Graphmix: Regularized training graph neural networks for semi-supervised learning. arXiv preprint drXiv:1909.117111 2019.
Blake Woodworth, Suriya Gunasekar, Jason D Lee, Edward Moroshko, Pedro Savarese, Itay Golan, Daniel Soudry, and Nathan Srebro. Kernel and rich regimes in overparametrized models. arXiv preprint arXiv:2002.09277 2020.
Han Xiao, Kashif Rasul, and Roland Vollgraf. Fashion-mnist a novel image dataset for benchmark- ing machine learning algorithms. arXiv preprint (rXiv:1708.0774:: 2017.
Gilad Yehudai and Ohad Shamir. On the power and limitations of random features for understanding neural networks. Advances in Neural Information Processing Systems, pp. 6598 6598-6608, 2019.
Difan Zou, Yuan Cao, Dongruo Zhou, and Quanquan Gu. Gradient descent optimizes over- parameterized deep ReLU networks. Machine Learning, 109(3):467-492, 2020a.
Difan Zou, Philip M Long, and Quanquan Gu. On the global convergence of training deep linear resnets, In International Conference on Learning Representations, 2020b.
13

---

Published as a conference paper at ICLR 2021
A PROOFS
In this appendix, we complete the proofs of our theoretical results. We present the proofs of Theorem Appendix A.I, Theorem in Appendix A.2, Corollary in Appendix A.3, Corollary 2 in Appendix A.4, and Proposition 1 in Appendix A.5. We also provide proof overview of Theorem in the beginning Appendix A.1.
Before starting proofs, first introduce additional notation used in the proofs and then dis- cuss alternative proofs using the implicit function theorem to avoid relying on the convergence of Neumann series.
Additional notation.
Given scalar-valued function R and matrix M Raxd" write
8M
aa OM
Rdxd'
E
0e 8Maa
aa üMaa'
where Mi represents the (i, entry of the matrix M. Given vector-valued function Rd and column vector E Rdd" we write
das 0b1
da1
Oa
E Rddxdd
bb
Baa
da.
ab
ūb
where b; represents the i-th entry of the column vector b. Similarly, given a vector-valued function a E Rd and a row vector E R xd' we write
da 8611
8a ab
0a Øb
E RRddd'
Dag
Bog
3011
0b,
where by represents the i-th entry of the row vector b. We use VAL to represent the map (A, B) + 81
8A (A, B) (without the usual transpose used in vector calculus). Given a matrix M and a function 4p, define Vario similarly the map M 0i (M). Our proofs also use the indicator function:
1 k) {1 0
ifi=k
if if
Finally, we recall the definition of the Kronecker product of two matrices: for matrices M E Idd xde and M E Rduxds,
Mn
Mid, M
M M
Tedarday xd
M
On alternative proofs using the implicit function theorem. In default proofs, utilize the Neumann series Ex= 7 (A)* when deriving the formula of the gradients with respect to A. Instead of using the Neumann series, we can alternatively use the implicit function theorem derive the formula of the gradients with respect to A. Specifically, in this alternative proof, we apply the implicit function theorem the function vi defined by
(<voo[A) A), 2) = 2 7a (A)= ((x),
where vec[A] and z R are independent variables of the function th: iee. vivve[[A], is allowed to be nonzero. On the other hand, the vector satisfying (b(ve[[A), is the fixed point
14

---

Published as a conference paper at ICLR 2021
2" tem 20) based on equation (3). Therefore, by applying the implicit function theorem to the function v. holds that if the the Jacobian matrix @ulive Al.e) 82
l== is invertible, then
O2
Oib (vb
(14)
i Vvec[A] [A]
O2
a vec[A]]
1==0
Since
(vec x)
= 1 - (A(( is invertible, it holds that
82
(vo(voc(A)). θ Ave[AA]
(I (A)))(()
(15)
vec[A] [A]
1==1
Moreover, since o(A)= E R is a column vector,
@e(vvc[[A] 2)
Oo (A)z
0 voc(o(A))
= θ vec[_A] A]
i @ve[[A] AA] ==
i vec[A]
s=
==
0|2 @ 1m] vec (ec) a vec[A]
1=2
->((>>) X 1m] θ @ecc(((A)) Ovec[A]
(16)
Combining equations (15) and (16), we have
Oz"
vec (
i vec[ (T - (A)) @ ml
(17)
8 vvoc[A]
In our proofs, whenever we require the gradients with respect to A. we can directly use equation (17), instead of relying on the convergence of the Neumann series. For example, equation (21) in the proof of Theorem is identical to equation (17) with additional multiplication Biu: i.c., for the left hand side,
Oz
0Bqez" O vec[A]
0. U-
B 2 vec[A]
0A
and for the right hand side,
a (e)(o(A))
y Bg (T - -A(( (=*)T [m Ovec[A] ) T
=> [ (MaLA)))) 8A.1
(Bn.U-1)T ((x) (U-T)a) ((U))()
(A(())) Ddem
) (Tx)U)))(((()(00
(x)))
A. PROOF OF THEOREM 1
We begin with proof overview of Theorem 1. We first compute the derivatives of the output of deep equilibrium linear models with respect to the parameters in Appendix A.1.1. Then using the derivatives, we rearrange the formula of V L such that is related the formula of Lo in Appendices A.1.1-A.1.3. Intuitively, we then want to understand V through the property of VLg, similarly to the landscape analyses of deep linear networks by Kawaguchi (2016). However, we note there that the additional nonlincarity a creates a complex interaction over the dimension m prevent from using such proof approach. Instead, using the proven relation of V and V Lo from Appendices A.1.1-A.1.3, we directly analyze the trajectories of the dynamics over time t in Appendices A.1.4-A.1.5 which results in a partial factorization of the iteration over the dimension Using such partial factorization, we derive the linear convergence rate in Appendices A. A.1.7 by using the PL inequality and the properties induced norms.
Before getting into the details of the proof, we now briefly discuss the property of our proof in terms of the tightness of a bound. In the condition of Bell < (1- (1 the statement of Theorem 1, the quantity (1 7)) comes from the proof in Appendix A.1.7: i.c.. it is the reciprocal of the quantity 1-y in the upper bound of [Tm (A)) & - Therefore, a natural question is whether or not we can improve this bound further. This bound turns out to be tight based on the following lower bound, The matrix Im (A) is Z-matrix since off-diagonal entries are less than equal to
15

---

Published as a conference paper at ICLR 2021
Furthermore, I.. (A) is M-matrix since eigenvalues of Im (A) are the eigenvalues of In -ya (o(A and the eigenvalues of - yo T are lower bounded by 1 11-- > This is because ({A) is stochastic matrix with the largest eigenvalue being one. Moreover, in the proof in Appendix A.1.7, we showed that / 1/-ra(Al (A))ss Eitj Il (A)lk = 1 1-y for all j. Therefore,
using the lower bound by Morača (2008), we have
- (()-)o()))=Il max, 70 E iØj 1 70 (A))1)
which matches with the upper bound of [I (o)) S - Therefore, cannot further improve the our bound on [Bll2 in general without making some additional assumption.
A.1.1 REARRANGING THE FORMULA OF V
We will use the following facts for matrix calculus (that can be derived by using definition of deriva- tives: e.g., see Barnes, 2006):
OM-1
da --M- -M
8.M da
M
Ou M
- ab M-T
0M
(M(M)
M aM
da
MMi da
0g(a)
0g(a) da 0M da 0M
Recall that U = ----70 (A). From the above facts, given a function we have
Ög(U)
m m dg(U) oug
AAN
aU; AAM
i=1 j=1
Ög(U) OUg, (A
aÜ (A)))
aAu
i=1 1=1
0g(U) (A
(18)
OUs dAu
Using the quotient rule,
(A)
exp(AA)
AAu
Dhn Ex exp(Ars)
Brxpl A,y)
)E,, exp(At))) - exp(AA;) BE, exp(Aes)
OA
(, exp(Ar)))
k(i=b)1lj= exp(eAApp An)) exp((xxp))))) (Aij)) - 1íj = expeepppppp exxp(A))))
(2 exp((AA) (1t))
1(i [i=b][[i=- = exp(An)
1 = () |eep(A (xp((As)exp((A)) exp(Ans
Exex((Ar) (1i))
exp(((Ar)))
exp(Aij) exp(Akj)
= = 1iii=)11( = ()0(A))) - 1 - E, (xx((Aes) E, exp(A19)
= 1 = 111 (A)( = -)))=)) (A)ki-
(19)
Thus,
0g(U) = - Σ
Ug(U) do(A)a - - (8g(U))
(Ao())) ER,
0AN
aUa aAH
0U11 dAw
16

---

Published as a conference paper at ICLR 2021
where Ug(U) 00.1
R and Ba(A)=
E R This yields
BAM
Øg(U)
(8g(U)
T (A)))
E IR 1x
0Aa
OAa
where An(A)ar BAat
E Rm
Now we want to set g(U) to be the output of deep equilibrium linear models as g(U)
Bq (lim-+>o (z, (2,A)))) A)) for each E (l...., my). To do this, we first simplify the formula of the output By (limi-200 (x, (x, using the following:
(Im- (A()
Σ ((()))
K=
= -m 70(A) (A) (r0( (o(1a)))-(((())-----((((((( (770
=1- - (>o())))
Therefore,
lim (x, (x, = lim (Im (A() Σ xth (a()())
(Im (o(A)
1:000
1+000
k==
Im
lim (AO( 6(x)
1000
6(x)
where the first line, the second line and the last line used the fact that (A)s 2 0.
[oo())l max Ll
and hence (A)]l for y (0. This shows that B (lim (x, Ax) BU -00-106 1(())
where the inverse U-1 exists as the corresponding Neumann series converges Σ ) (A)* since (A))l1 < 1. Therefore, we can now set g(U) = By (lim (0) (x,A)) A)) = 1 ((z).
Then, using
M-
- M M-T.
8M
0g(U) au
θB..U U ['(z) -Bgh (x(z) TU-T
aU
which implies that
a. U- - (B₂z) Txx=)=--T
-UU-T(Bpp
0U.1
-U-T((B,.) ((x) (U-T), R
(20)
Combining (18) and (20),
0g(U)
-
(g(UU)
T (o(())1
0A,
0A,i
aU., 0A,i
(ao(A))s
=7 -U- ((Br)) (B..)T((a TT(x)T (((T)((
AAal
(ao(A))L
=>((U-T) (MaUTT)a)T&a (
U-1
(
dA
17

---

Published as a conference paper at ICLR 2021
+=(U~")nb(2)BBBp0 (Do(A))
E IR lxm
AA.
where we used (+ ((U-1)TH ((U-))),)T and ((U-T)+) (((U-)).)))TT= (U-)) By taking transpose,
@Bp.U-ll(x)) =7 ( (oo(A)))) T (B 6(x)) (()))((UTT)) E Rmx!
0A
AAu
By rearranging this the matrix form,
(21)
0A
=> -r((2nAn) ((a(A)1)) JA.1
)
T
Or ( Jd..m
((x) (U (U-T)(n))))
mm
where @Rg o(x) 0A E Rmxm Each entry of this matrix represents the derivatives of the model
output with respect to the parameters A. We now use this to rearrange V VAL(A, B) = 2L(A.B) We set Big = B-, ((x) and 91 BU o(x) and define
0A
Oo( (A)aa aA
Jk
E Rm
and
(iii Ui)
T
Q =
d(I;) E Rm
Then, using the chain rule and the above formula of (B,
0A
(L(A,
0A
n my
al(gr, y.)
0ūgg
Bggg
0A
de(ii, Vi)
[ (B ((;)) ((:)))(UUTT)) UB T((x;) T )
aggg
B 1)
)o(x;) (U (U-T)a
J EB
)o(2;) (U ((U (U-T)
Bigg
Bise
(5 (BB
)o(a,) 1 (U (U-1)a1
J [""r ).g
ar )o(z,) (U- )am
8y
J. (BU-1) 0%
TT(2:) ((U-T)
086 YT J ( ov.
=> [ (Bold) (BU-))TQ(U-T)
T (BU-))TQ(U-)
(Al AA.m
Summarizing the above, have that
V AAL B) = AL(A, B)
0A 7 [] (BU(UU)))QUU-))
JT (BU-1)T QU--T)QUU))))))))) )mm]]
(22)
where (ATAA, B) E Rmxm.
18

---

Published as a conference paper at ICLR 2021
A.1.2 REARRANGING THE FORMULA OF Vo
In order to relate Lo to the gradient dynamics of L, we now rearrange the formula of Lg. We set Dig (xx) E R and D = W (xz) E "Ry for linear models. Then, by the chain rule,
@Lo(W) oW
dé(ii. yi) 0gig aw
=
Dili
Since
Ogjn
1 = g)o(x;)
OW
we have
al(ūr.v) Disk
allúr- 34) Dün
o(x;)
aW
Disk
ow
By rearranging this into the matrix form,
n 1=
(bT;))
Jūu
@L)
aW
ae(
(bxxa)
i=1
Oiimy
aH(U..99) .1))
(2x))
V6(g..ve) (1x))
Jylmyy
0e (8(().)))) M.)
((x,)
o(z;) E gm
Oyi
==1
where (.. Oy.
E 1xm, Thus,
VLo(W) = OL(W)
Oe( Ui) ) o(x;) E Rm
(23)
aw
0g.
A.1.3 COMBINING THE FORMULA OF V AND VLD
Combining (22) and (23) by resolving the different definitions of 0; yields that
VAL(A,B)
(24)
=y (BU-
JJ ((BU-))TV
Here, if there is no additional nonlinearity a, the matrices Jj (o(A)., become identity for all
8A.x
k. In that case, VAL(A, B) can be further simplified and factorize over m, which is desired for the analysis of gradient dynamics. However, due to the additional nonlinearity, we cannot factor- ize (AL B) over m. One of the key techniques in our analysis to keep this un- un-factorized V LAA, B) and find way to factorize it during the update of parameters (A,, B, in the gradient dynamics, as shown later this proof. To do so, we now start considering the dynamics over time
1.4 ANALYSING (lim) 20) (ax, At))
Now let us temporarily consider a gradient dynamics discretized by the Euler method as
A A ALA,, B.),
19

---

Published as a conference paper at ICLR 2021
with some step size 0. Then,
lim 2 (ctt)) = (4m (00(1++1) =0(x) 10000
= (Im -o - B.))(A,
where used (limp-so (x, U,1 x(x)) from Section A.1.1. By setting Piy(a) o(A aV aV,L(A, B) E I.
o alALaVaL(A B)), = Piy(a) = Piy(0)) + 2p a + O(a²).
da
By using the chain rule and setting M = - LAL ER**
@ppi(a)) da
m m
(M) dMR
OM
da
k=
1=1
Σ
Σ
1 (j = (i= - 1 1 ()) VAT(A, B)
Σ k=1
L = k)) J (AL
Therefore,
(0))(())
>
(A)))O(A)) ())(A,
da
R=1
(A) AL,
AAki
Oo A)ij
V (A, B).j ER.
OA.j
Recalling the definition of Jg = aA) BA.k
E R m
(0), (o(A).s
V VAL(A, = = B),, E mx
Oox
0A.j
Rearranging it into the matrix form,
0p(0)) JJV L(A, B).1
Jn V ALA, B)mm] gm
da
Putting the above equations together,
a - (A B)) = -(0) (p) 0(a²)
da
= =o(A) a A
Jn [mm
Thus,
Im a
- - L(A
Jm A B).m +
[Im 10--000 (o(A) +70 IV B),
V O(a3)]
[U ar
B).1 .1
Jm A ]m] + +(aa)]
By setting M [ A(A, B) Jn VA mm) and ((a) = (U + ay + -1
and by using aMd-1 M-
AM M-1
(A ((A, B)) ao +0
20

---

Published as a conference paper at ICLR 2021
p(a)
(0) + ag(0)
da O(a²)
U-! anu MU - 2a0(a) + O(a?)
- MU- +0
Summarizing above,
I - (A B))]
(25)
Jn V () O
A.1.5 PUTTING RESULTS TOGETIIER FOR INDUCED DYNAMICS
We now consider the dynamics of 2 = B,U," in my xm that is induced by the gradient dynamics of (Ae, B:):
0L (A,, B.),
AL (A. (A, B,), vt 2 0.
dt A 0A
B, = dB
dt
Continuing the previous subsection, we first consider the dynamics discretized by the Euler method:
Z++1 = B4+1U [B: - a [A,, B.)][m (a,, aVA,,A, B.))]
where a > 0. Then, substituting (25) into the right-hand side of this equation,
Zt+1
B,+1 UNU, [J,, AAL B:)a1
Jm A Bm).m + 0(a
B4+1U,1 ary B1+1 [J1, VALA,, B:).1
Jm, V ALL(A, Be) 0
Using Bi+1 = [B: a BL B.),. we have
UrB++1U, [J1, V (Ae, Bt)a1
Jmn, V A Be) Ui
=,y,,UTT [J1, AL(A, Bi).1
Jm,t V T. Be) am) Ur + +(a2)
= = [11, V AL(A, Be)+1
Jm, V mm)] U. + O(a2).
1
Since (+- ((U-1)T)+* -.
U-1
and V (A,
(L((A
=7JE (BU-1) (BU from (24), we have that )ak
0A
Z, [M,, ALLI,, B.).1 Jm V AL(i, Bm) U,!
=ZZ,
(ZO(Z)) ( )
JT, Z, Lo (Z))(((U "mm]] Ue
(U,
14
=7 [J z V VLo TM
Jm. 2, (LO(Z))(U 7mm))
)mm)
m
Z.Jk.J 2, (Lo(Z))(U (U,
Σ Z.Jae k. 2 (.,
On the other hand, using B = [Br aVnL(A,,B.) and VnL(A, B) BL(A.)) = 018
(BE()..)) o(x,) )(x))) -U-T = (IU(BU- we have Ogi
21

---

Published as a conference paper at ICLR 2021
B Z (Z))U U U- =x=1 (U )ak
Summarizing these
equations by noticing
(U )kk= yields that
2+1 Z - -VVol(Z))
Σ (UT (0,( )x- - -o?? Σ
VLo(Z))(T )oR(UT") )2. + +((a)) ZeJaNJ 2L
k-1
R-1
= 21 -a
VLO(ZNNU (o(ZNNUTT<<((U))) + 4 VLo )oa(U,)))))))))) ),. + +0(aP)
m
Z-
+ 2. )oRU;,)))))))))))))) + +0(a7)
By vectorizing both sides,
vec(Z++1)
=vec(Z,)
vec (Lmm 122 + +0(a²)
k=1
-vec(Z,)
(((TT))a(T )on(UT (Im, + 2, vec(VLolZ))++O0
By defining
D=:= Σ
@ (Im + Ze 2,)],
we have
vec(Z++1) vec(Zr) - aD velovvee
This implies that
vec (Z++1) - vec(Zr) Dr vec (L)( O(a):
where > O. By recalling the definition of the Euler method and defining Z(t) Z₁ we
rewrite this as
vec( ((t + a)) vec (tt) = D, vec vec (L,)Zz)) + O(a).
By taking the limit for a - 0 and going back to continuous-time dynamics, this implies that
dt vec(Z.) = D; voc V Lo())))
(26)
Here, we note that the complex interaction over m due the nonlinearity is factorized out into the matrix Dr. Furthermore, the interaction within the matrix D, has more structures when compared with that in the gradients themselves from (24). For example, unlike the gradients, the interaction over m even within D, can be factorized out in the case of my = 1 as:
D (U=T))>(U~\)e (mm, + 2
(I+ 2 ZJ LT) (U- )A )m+
1 22 ZT
U-T diag
1
U-1
2 Jn TT
Z2 ZT
= T Im + diag
:
U-1
>22Jm JT
Although do not assume my 1, this illustrates the additional structure well.
22

---

Published as a conference paper at ICLR 2021
A.1.6 ANALYSIS OF THE MATRIX De
From the definition of Dr, we have that
D2 [UT (,; @ (Imu t
(w-T))a(U~)))))))) [mm] +
(a- 2 Of 2,]
m
(U (0 (U= Am
)- 2 2
U TU-TU-1 Imy] +
Σ (U Jun(U kek 80² Z
(27)
Since U-TU-1 is positive definite, I,, positive definite, and Kronecker product of two positive definite matrices is positive definite (since the eigenvalues of Kronecker product are the products of eigenvalues of the two matrices), we have
[U-TU- 1m,,] O.
(28)
Since (+= is positive semidefinite, Ze Kronecker product of two positive semidefinite matrices is positive semidefinite (since the eigenval- ues of Kronecker product are the products of cigenvalues of the two matrices), we have
k, ZT is positive semidefinite, and
[(U (U- @ > J.JK,. Z ž 0.
Since sum of positive semidefinite matrices is positive semidefinite (from the definition of positive semi-definiteness: x M, 2 0 = x
(Σx Ma = Σ₁ x M 2 0),
m
ICU )OR (U Z. J
T 2
O.
(29)
Since a sum of positive definite matrix and positive semidefinite matrix is positive definite (from the definition of positive definiteness and positive definiteness: (x: Mi > 0 NxT M2x) = T (M1 +
M2)r x. Mi + x M2 > O),
m
D. =
)ak (U
ka )u.
(Im, + +22 22 J. (T,
[U TU-TU-L Imy] + > ((U & 222 2i]] 0.
Therefore, D, is a positive definite matrix for any and hence
AT inf xE(0,T]
Amiu (D.) >
(30)
A.1.7 CONVERGENCE RATE VIA POLYAK--OOISSEWIIC INEQUALITY AND NORM BOUNDS Let R E (0, oo] and T > 0 be arbitrary. By taking derivative of Lo Zi) LØ.R with respect to time t with Z, == 1
d (Lo(Zi) L6,R)
my
X92
d.Lo
d (2))
d
(Zi)
dW
dt
is
i3
dLo
dW (2i)
c
d
(Zi)
di
ij
ij
23

---

Published as a conference paper at ICLR 2021
where we used the chain rule and the fact that d (z,) = deo (Z)),
di L O. By using the vectorization notation with
d
dt (Lo(Zz) Lū,n)
vee vcc (VL)(2)))
vec di
(2i)
By using (26) for the equation of vec [at (Z+)].
d R [V+)(Z+))
dt (To(Ze)
D3 vec vec (Lo((Z))
S - Amin vee Lo l1
= -mmin VLo)
Using the condition that Vo satisfies the the Polyak-Lojasiewicz inequality with radius R, if 2 < R. then we have that for all t E [0, T).
d (Lo(Z))- - n..n) S -x (Di)(L - 1.)
dt
S - Lo (Zi) LigR).
By solving the differential equation, this implies that if || < R,
Lo(Zr) - Lū,n S (Lo(Zo) - Lo.n) e -2ndrT
Since L(Ac, B.) (Z())) if |Ze||1 < R. L(Ar, Br) S LØ (L(Ao, Bo) RØ 2
(31)
We complete the proof of the first part of the desired statement by showing that [| < (1- It implies 2 < R. With Z BÚ-1, since any induced operator norm is submultiplicative matrix norm,
| 11 [ - -A)(
111 (A))
We can then rewrite
[I - xo())
II((Im - (A()) = 11(4m - (AA) '1loo.
Here, the matrix 4m - yo (A) is strictly diagonally dominant: i.c., 1m (A) Ii: Ejdi Im (A)ley for any i. This can be shown as follows: for any j.
1>7+ 1> Σ (A)i
+ 11 1> 70( +
Σ (A)is
(
+ 1 - (A))s
Σ
yor (A)is
+ 4m (o(A))
-A)(A)
ifj
+ Im - ro(A))bs
Σ
Im (o(A))
i
+ Im (a() 185
Σ
I'm (A) 1/3
This calculation also shows that IT, - (A)b;-- Im - = 1 for all j. Thus, using the Bulberg-Nilson-Weee bound for the strictly diagonally dominant matrix (Ahlberg Nilson, 1963; Varah, 1975; Morača, 2008), we have
1 [1(m 7) miny (T (A())) - Σitz I' (A(A)))))
1
1 1-7
24

---

Published as a conference paper at ICLR 2021
By taking transpose,
1
I1 [(mm=) 1111 -
1 1-7
Summarizing above,
1
1| [I - (A() 3 ||B||
1-7
Therefore, if B|]|| < (), then ||2||1 = [B (A)) 111 ||B|| 1 < R, as desired.
Combining this with (31) implies that if [ < R(1 77)
L(A., Bd) S Lo.R + (L(Ao, Bo) - R) 2xATT
Recall that T inf [nfll To (W) and 1. inf AERmXm BEBB [(A, B) where BR = ( B
2T |B]||| < (1 R}. Here, B E BR implies that || | ||| S [B] 1| <
(1 R|I |11 S R. using the above upper bond ||U~~1|2 [I yo || 5 I-7'
Since
L(A, B)- = Lo(Z) with Z = BU-1, this implies that LōR S R and thus
L(A,, Bb) S Lô + (L(Ao: Ba)
R)
2 L + (I(Ao, Bo) LØ 2RATT"
This completes the first part of the desired statement of Theorem 1.
The remaining task is to lower bound λr, which completed as follows: for any (A, B).
Amin (ainnD) (D) = min T Du
oz [|||o||11
m
Σ I(U ((UT)))R(U @ (In, Z 2,,1 v
min
min v [U-TU-1 & mmm]o
e: ==||v||=1
Aminn (U
- Jm,))
>minnU U--1
Amii
1
= ouin (1-1))
(32)
1111 m
where the third line follows from (27)- the fifth line follows from the property of Kronecker product (the eigenvalues of Kronecker product of two matrices are the products of eigenvalues of the two matrices), and the last inequality follows from the relation between the spectral norm and the norm || 1111. We now compute ||6/|/ as: for any (A, B).
||U|| 114m (A)il
= max Σ 1 - (A))isl
max Σ El(Tm)ks - (As(s)
-max (4m)ss - (A))ss) +
Σ
(1m))s - (A)
)
i
max |1 1- (A))s)) + +Σ - (A(v)ll
itj
max1 max1-~o(AA); + Σ (A))s)
(AA)s [[o(Aba)
=max 177 +1 j
-AA)ss
<1+9
25

---

Published as a conference paper at ICLR 2021
By substituting this into (32), we have that for any (A, B) (and hence for any t),
1
Amin (D)
(33)
m (1 + 7)²
This completes the proof for both the first and second parts of the statement of Theorem 1.
A.2 PROOF OF THEOREM 2
We first show that with 6, > 0 sufticiently small, we have Ge > 0. Recall that
G = Uz (S, Fi) U,T = U.S, U, - S,U,F,U, Thus, with 6, > 0 sufficiently small, for any + O,
1 Gv = UU,,, U,
de U,F.U, e,
which dominated by the first term v U.S, T v if the matrix U.S U.T is positive definite. Since S. == Im + diag(e)) with of E R and (o) ) == [B, 113 for = 1, m, the matrix U, is positive definite. Thus, with 6; > 0 sufficiently small, e' Gv is dominated by the first term, which is positive (since U,S, U, is positive definite), and thus we have G. > O.
Then we observe that the output of argmin oclivlle, S6,ll illl # veo(B,U,)))le (v(v) of the following constrained optimization problem:
is the set of solutions
minimize (v) s.t. lella 82 2 (B,U, 0.
G,
Since this optimization problem is convex, one of the sufficient conditions for global optimality is the KKT condition with a multiplier μ E R:
VLg(v) 2pG,v
#
|1112 - o6 d veel B,U,-
2
O.
dt
G.
Therefore, the desired statement is obtained if the above KKT condition is satisfied by 5, ( vec(B,U,~)) with some multiplier μ.. The rest of this proof shows that the KKT condi- tion is satisfied by setting 6, &d vec( () and 25, With this choice, the last two conditions of the KKT condition hold, since
1 2 0,
28:
and
d
2
2 = 0. - 55 d vec vec(BBU,))
|ul12, - 88 vec (BU,,
vec
= 5
vee B.U,
dl
G.
G.
Gr
The remaining task is to show that VLi(v) 2uGiv f
with 6, (4 vee vec(B.U,))) and
25. From the definition of Thy-
Vlg(o) + 2pG V vec (vec (BU,-)) + I (vecl B,U, ( 2pGrr.
(34)
We now compute and VLyee and 2 Lyec. Since (L)(W) aw E aLo(W)
(Be(g).)))
x..
0g:
(E yi)
(
(al(0e vi) T
vec (Lo(W)) = Σ
vec
Imy
MI;)
= > (p(xz)) @ Im,,
Oğs
1=1
26

---

Published as a conference paper at ICLR 2021
where
0 : Wo(z;;) = [b(();)) @ Imy] vec[W]
Therefore,
"
(
(Delin 0g:
Vi)
T
V L (vo veo voc voc(VLo(WW)
Ei((s Imy]
For the Hessian,
i d ve)(W)
2 V2L 12 (vec(W))
V (vec(W))
0
(OE( vi)
Oi
Σ [o(xs) @ Imgg]
Dis
Oūs
a vee(W)
θ
@E M1)
Σ (xx) 60 [mg]]
(ob(az) s mmg)
Oij
==1
By defining (x) ((z, vi) and (=((()) =
0
(6,(2)
0z
0z
4 Lg (vec()))
Σ @ I,,
( ([ @ Imy]
(35)
Since we have that
1 (Im (A))
E
(()))))) =m - (A( + (A) (r0( + (ro(2))) - (A))3)) - (A)
= -
we can write:
1
(Im - (A())
lim 200 (x, A) 1+200
lim (Jm - (o()) 10:0
Σ {
Lm
lim (A)()))+ o(x) 1000
((x),
where we used the fact that (A)) 2 O, [[o(A]]1 max, E, (A(A)b) 1. and thus (A))l < 1 for any y E 10. This shows that (lim 200) (z, A)) = 2" (x. A) = U-1)(x), from which we have Φ(x;) = U2 (z,, A).
Substituting d(x;) Uz U2, A) into (35),
(vec(ve a
Uz (I A) 2 (VV @ I,yl.
Σ
U & Im, 112 A) & Le, (2,,
T @ Lmgl.
1=1
U0 @ 1
[2 (1,, A) 2 1
10 ② m. [mmg]
27

---

Published as a conference paper at ICLR 2021
In the case of my = since I,,
1, (vec(W)) is further simplified to:
L0 (v)c(WW) = U
e, (Woo(())) (rz, A)2 (2,, A) UT.
Therefore,
LO (v
where
F.
((s)) Az)=" A
E=1
By plugging μ 26,
and V2 N (vecl - UFFUU into (34).
1
VLI(U) + 2uGrv= (vec( U.FUU + U. (S 1 5,Fi) U,T
8;
1
(vec((() + U.S U,
8,
By using u = 5, (áe (oc(BBU,11)
(T) + 2pG,v= V.oo" U.S U, ( voo(BB);)))
By plugging (26) into d vec(B with Z. B,U
VTg(o) 2pG Lvec - U. Dr vo)(VIn(Z))).
Recall that
D. Σ (L & (Lm, + 2 )..
In the case of my = 1, the matrix D, can be simplified as:
D
[(U-T) )o-())= @ (1m, + + 21)
(1 + 7J TT) (U-T)ab (U )k
=
1 + 72 ZJJJJ zT
= -U- diag
U-1
+ 72 22 2T
- S
Plugging this into the above equation for Dr,
(u)()) + 2 = (v) v (vec
(vec 0.
Therefore, the constrained optimization problem at time / is solved by
v = 5, ( dt
veol veo(BBu;)))
which implies that
d
1
veol (BU = voo(VV), vec(V))
argmin me]le|la,Sa 1 vec(
Ib((),
dt
5c
28

---

yielding that Published as a conference paper at ICLR 2021
By multiplying <(x) @ Im each side of the equation, we have
d (8(x)*
x Im,) veol (B,,, = B,
lim
(x)
dt
1
1
[0(x) Imyl voc(V))
6e Viń(1x),
6,
lim (x. A.)) = Vio(z). dt 10000
d B,
1
8e
This proves the desired statement of Theorem 2.
A.3 PROOF OF COROLLARY
The assumption rank( () min(n, m) implies that m ( > 0. Moreover, the square loss l satisfies the assumption of the differentiabilityy Thus, Theorem 1 implies the statement of this corollary if Lo with the square loss satisfies the Polyak-Lojasiewicz inequality for any W E Xm with parameter K. = 20mia 2immin (0) This is to be shown in the rest of this proof. By setting yp Lo in Definition 1, we have || VVvvec(voc()))
(W)) With the square loss, we can write 11 [W 4 - Y where & E Rmxx and Y E Rm xn with
Lo(W) = Ei=1 | ((x;)
3 = (xr)) and Vj= (Wi),. Thus,
(W) 2( r) E 1
We first consider the case of m S n. this case, we consider the vectorization Ly (vec(W)) = Lo (W) and derive the gradient with respect to vec(W):
L (vec())) = 2 vec( (W Φ - r) T) = =2 @ [my] vec( 1 Y). Then, the Hessian can be easily computed
= 2 @ 4m, T @ Lmm,] = 2 T 8 Amyl
where my is the identity matrix of size my by my. Since the singular values of Kronecker product of the two matrices is the product of singular values of each matrix, we have
W)) 2 20 II
where we used the fact that m < n in this case. Since W arbitrary, this implics that Ly is strongly convex with parameter 20 (4) > 0 in Rm Since strongly convex function with some parameter satisfies the Polyak-L.ojasewiize inequality with the same parameter (Karimi al., 2016), this implies that L (and hence Lo) satisfies the Dojssk-Lojasiewiez inequality with parameter ma (+)² in xm
xm in the case of m S n.
We now consider the remaining case of m n. In this case, using the singular value decomposition of Φ- UVV
VLo(W)]1 21 T) -Y
2
21|UEVT (W - T
= 2/JEV (W - T
2 20min (2) (W ( T)
20 LoW)
20 ()2 (Tn(W) T,5*)
for any Li 2 O, where the first line uses ||q/l= = ||g 1|2. the second line uses the singular value decomposition, and the third and fourth line uses the fact that U and V are orthonormal matrices. The forth line uses the fact that m > n in this case. Therefore, since W is arbitrary, we have shown that La satisfies the Polyak-Lojasiees 20 2 in both cases of m > and m < n.
inequality for any W Xm with parameter
29

---

Published as a conference paper at ICLR 2021
A.4 PROOF OF COROLLARY
The assumption rank ④) m implies that min (4) > O. Moreover, the logistic loss l satisfies the assumption of the differentiability. Thus, Theorem 1 implies the statement of this corollary if Lo with the logistic loss satisfies the Polyak-Lojasiewice inequality with the given radius R E (0, oo] and the parameter
(2T + p oo (Φ)² where p (R) depends on R. Let
(0, oo] be given. 0 if R 20. If (R,
(oc, 0), then 2T + p(R) = for which the statement of this corollary trivially holds (since the bound does not decrease). Therefore, we focus on the remaining case of (R, r) * (oo, 0). Since (R, 7) * (oo, we have 2T + p(B) > O. We first compute the Hessian with respect to as:
Note that we have p Rt) > 0 if R < co, and p
12
LoWW) Σ P.(1 (1)(()))) 2T Σ (x;))((N))
where i. Ite (oo- Therefore,
v. L' = Σ pi(1 - p-)u ((x))((Z)) To 2 T
Σ (=ss)((a))
2 (() Σ U. T ((x))(()))) Tu u 2To
Σ (oss)((a))
= (2r + p(R)) TTT
> (2r +
0.
Therefore, Lo is strongly convex with parameter (2r J (Φ)² > Since strongly con- vex function with a parameter satisfies the Polyak-tooasiewice inequality with the same parameter (Karimi al., 2016), this implies that Lo satisfies the Polysk--tojssiewiez inequality with the given radius R E co] and parameter (2r + p( (4)² > O. Since R E (0, oo] is arbitrary, this implies the statement of this corollary by Theorem 1.
A.5 PROOF OF PROPOSITION 1
Let (z, A) be given. By repeatedly applying the definition of (x, A), we obtain
20 A) 70 (A(=))(11) ((x)
= (A) + +(x)) + 6(x)
(36)
k=0
where a (A)* represents the matrix multiplications of i copies of the matrix a (A) with (A)) I In general, if a is identity, this sequence does not converge. However, with our definition of have
|o(A)] -max Σ
(A) 11
Therefore, (A)]l (A)) = 7 for any 7 (0, 1). Since an induced matrix norm is sub-multiplicative, this implies that
(A)*th) II
In other words, each term (>)o(A)*1 in the series Ei= (A)**())))) is bounded by * Since the series " 7² converges in with y E (0, 1), this implies, by the comparison test, that the series Σe= 11 (A)**11 converges in R. Thus, the sequence (E * is Cauchy
30

---

Published as a conference paper at ICLR 2021
sequence. By defining 81 Ek=o y (A)*, we have ||s: |s Ex=r+1 7 or (A) I Ex=rt (A)*] for any 1' > by the triangle inequality of the (matrix) norm. Since Ei (A) is a Cauchy sequence, this inequality implies that (Ei-o 7 (A)*)) is a
Cauchy sequence (in Banach space (Rm*m, II II1), which is isometric to Im under |- (11 Thus, the sequence (Ek=o > o (A)²), converges. From (36), this implies that the sequence (z (9) (a, A)), converges.
B ON THE IMPLICIT BIAS
In this section, we show that Theorem 2 suggests an implicit bias towards a simple function as a result of infinite depth, whereas understanding this bias in more details is left as an open problem. This section focuses on the case of the square loss E(g. yt) = ||g y: I12 with my = 1. By solving vec (V)) E argmingey (v) in Theorem 2 for the direction of the Newton method, Theorem implies that
V vr- - (A()) E E
(37)
where ra E R is an error vector with each entry being a function of the residuals fo, (2:)) ys-
Since (A) is positive matrix and is a left stochastic matrix due to the nonlinearity o, the Perron- Frobenius theorem (Perron, 1907; Frobenius, 1912) ensures that the largest eigenvalue of (A,) is one, any other eigenvalue in absolute value is strictly smaller than one, and any left eigenvector corresponding the largest eigenvalue the vector 7 1/[1, 1......1 E R.m where Ç E is some scalar. Thus, the largest eigenvalue of the matrix (4m (A.)) is T-7. any other eigenvalue is in the form of 1- with Ave < 1. and any left eigenvector corresponding the largest eigenvalue is 1y E Rm. By decomposing the error vector Tı Pir: + (1 Pi)r:, this implies that (V,) = V, = - Pirt + 9, ((1 - P Pi)r.) E IR", where P, - 11 is a projection onto the column space of 1 = [1,1,...., 1] E RTm and 9, is function such that for any q in its domain, [1)(()) < e for all y E (0,1) with some constant e in 7/
In other words, vec(V)) in Theorem can be decomposed into the two terms: - Pir, (the pro- jection of the error vector onto the column space of 1) and gy((1 - Pi)r:) (a function of the pro- jection of the error vector onto the null space space of 1). Here, as 7 + 1, -1 Pirell + oo and (1g.((1 Pi)r,)ll < c. This implies that with y < sufficiently large, the first term Pire dominates the second term 9,(((-Pi)r))
Since (Pirt) (xx) = Es T ((x) with je = Exe E R, this implies that with y < suf- ficiently large, the dynamics of deep equilibrium linear models f. = +o+ Vio learns a simple shallow function år o(x) first before learning more complicated components of the functions through 9, Pi)re). where po = S p,dt E F. Here, Az
m
simple average model that averages over the features 113 T Φ(x) and multiplies it by a scaler f.c. Moreover, large 7 < means that we have large effective depth large weighting for deeper layers since we have a shallow model with = 0 and y is a discount factor of the infinite depth.
m 1 Tzz())
C ADDITIONAL DISCUSSION
On existence of the limit. When Bai et al. (2019a) introduced general deep equilibrium models, they hypothesized that the limit lim 2!!) exists for several choices of h, and provided numerical results to support this hypothesis. model parameters, the limit is not ensured to exists. For example, when h increases the norm of the output at every layer, then it is easy to see that the sequence diverge or explode. This is also true when we set h to be that of deep equilibrium models without the nonlinearity (or equivalently redefining a to be the identity function): the operator norms on A are not bounded by one, then the sequence diverge in general. In other words, in general, some trajectory of gradient dynamics may potentially violate the assumption of the existence of the limit when learning models. In this view, the class of deep equilibrium linear models is one instance of general deep equilibrium models where the limit is guaranteed to exist for any values of model parameters as stated in Proposition 1.
general deep equilibrium models, depending on the values of
On the PL inequality. With R sufficiently large, the definition of the PL inequality in this paper is simply rearrangement in the form where Lo is allowed to take matrices as its inputs through
31

---

Published as a conference paper at ICLR 2021
[vec()))))) Lo(-), where vec(M) represents the standard vectorization of a matrix M. See (Polyak, 1963; Karimi et al., 2016) for more detailed explamations of the PL inequality.
On the reditus R for the logistic loss. As shown in Section 3.1.3, we can use R = oo for the square loss and the logistic loss, in order to get prior guarantee for the global linear convergence in theory. In practice, for the logistic loss, we may want to choose R depending on the different scenarios, because of the following observation.Fo For the logistic loss, we would like set the radius R to be large so that the trajectory on B is bounded as B.ll < (1 R) for all E [0, T] and the global minimum value on the constrained domain to decrease: i.e., L3 L' R Co. However, unlike in the case of the squared loss, the convergence rate decreases as we increase R in the case of the logistic loss, because p( R) decreases R increases. This does not pose an issue because we can always pick R < oc so that for any t > 0 and T > O, we have p R)) > Cp for some constant Cp > O. Moreover, this tradeoff does not appear for the square loss: i.e., we can set R do for the square loss without decreasing the convergence rate. We can also avoid this tradeoff for the logistic loss by simply setting R oo and r > 0.
On previous work without implicit linearization. In Section 5, we discussed the previous work on deep neural networks with implicit linearization via significant preepparameterization. Kawaguchi & Huang (2019) observed that we can also use the implicit lincarization with mild wver-parameterization by controlling learning rates to guarantee global convergence and general- ization performances at the same time. On the other hand, there is another line of previous work where deep nonlinear neural networks are studied without any (implicit explicit) linearization and without any strong assumptions; e.g., see the previous work by Shamir (2018); Liang et al. (2018); Nguyen (2019); Kawaguchi & Bengio (2019); Kawaguchi & Kaelbling (2020); Nguyen (2021). Whereas the conclusions of these previous studies without strong assumptions can be directly appli- cable to practical settings, their conclusions are not strong as those of previous studies with strong assumptions (c.g.. implicit linearization via sigmificant weerparameeterizttoon as expected. The di- rect practical applicability, however, comes with the benefit of being able to assist the progress of practical methods (Verma et al., 2019; Jagtap et al., 2020b;a).
D EXPERIMENTS
The purpose of our experiments is to provide secondary motivation for our theoretical analyses, instead of claiming the immediate benefits of using deep equilibrium linear models.
D.1 EXPERIMENTAL SETUP
For data generation and all models, we set ((x) = x. Therefore, we have m = my-
Data. To generate datasets, we first drew uniformly random 200 input images from stan- dard image dataset - CIFAR-10 (Krizhevsky & Hinton, 2009), CIFAR-100 (Krizhevsky & Hinton, 2009) or Kuzushij-MNNIST (Clanuwat et al., 2019) - as pre input data points -/"Te E RTT Out of 200 images, 100 images to be used for training were drawn from train dataset and the 100 other images to be used for testing were drawn from the corresponding test dataset. Then, the input data pints x; E IR' with my = 150 were generated as z; R.." where each entry of matrix R. E Rmm,xmms was set to 5/ /mz with 5 gldd NNO0,) 1) and was fixed over the indices i. We then generated the targets as y; = B" (lim/++00 20) (x;; 55, E R with y = 0.8 where 8, lid NN(,, 1). Each entry of the true (unknown) matrices A R. xm and B" E R Rmmm was independently drawn from the standard normal distribution.
Model. For DNNs, we used ReLU activation and (1 R for 1 1.2...,,I (W(H) E RIxm). Each entry of the weight matrices (!) for DNNs was initialized to 8/ /m where 1id. (O 1) for all / = 1, 2..... H. Similarly, for deep equilibrium linear models, each entry of A and B was initialized to mm where 5 li A/ 1). Linear models were initialized to represent the exact same functions those of initial deep equilibrium linear models: i.c., Wo BoU We used 7 = 0.8 for deep equilibrium linear models.
Training. For each dataset, used stochastic gradient descent (SGD) to train linear models, deep equilibrium linear models, and fully-connected feedforward deep neural networks (DNNs). We used the square loss {(g, vz) ||g v:11.. We fixed the mini-batch size to be 64 and the momentum
32

---

Published as a conference paper at ICLR 2021
coefficient to be 0.8. Under this setting, linear models are known to find minimum norm solution (with extra elements from initialization) (Gunasekar et al., 2017: Poggio et al., 2017). Similarly, DNNs have been empirically observed to have implicit regularization effects (although the most well studied setting is with the loss functions with exponential tails) (e.g., see discussions in Poggio al., 2017; Moroshko et al., 2020; Woodworth et al., 2020). In order minimize the effect of learning rates on our conclusion, we conducted experiments with all the values of learning rates from the choices of 0.01, 0.005, 0.001, 0.0005, 0.0001 and 0.00005, and reported the results with both the worst cases and the best cases separately for each model (and each depth H for DNNs)
All experiments were implemented in PyTorch (Paszke al., 2019).
D.2 ADDITIONAL EXPERIMENTS
In this subsection, we report additional experimental results.
Additional datasets. We repeated the same experiments those for Figure with four additional datasets modified MNIST (LeCun et al., 1998), SVHN (Netzer et al., 2011), SEMEION (Srl & Brescia, 1994), and Fashion-MNIST (Xiao et al., 2017). We report the result of this experiment in Figure 4. As can be seen from Figures 4, we confirmed qualitatively the same observations as in Figure 1: i.e., all models preformed approximately the same at initial points, but deep equilibrium linear models outperformed both linear models and nonlinear DNNs in test errors after training.
DNNs without bias terms. In Figures 1-4, the results of DNNs are reported with bias terms. To consider the effect of discarding bias term, we also repeated the same experiments with DNNs without bias term and reported the results in Figure 5. As can be seen from Figures 5, we confirmed qualitatively the same observations: i.c., deep equilibrium linear models outperformed nonlinear DNNs in test errors.
DNNs with deeper networks. To consider the effect of deeper networks, we also repeated the same experiments with deeper DNNs with depth H 10, 100 and 200, and we reported the re- sults in Figures 6-7. As can be seen from Figures 6-7, we again confirmed qualitatively the same observations: i.e., deep equilibrium linear models outperformed nonlinear DNNs in test errors, al- though DNNs can reduce training errors faster than deep equilibrium linear models. We experienced gradient explosion and gradient vanishing for DNNs with depth H 100 and H 200.
Larger datasets. In Figures 1, we used only 200 data points so that we can observe the effect of inductive bias and overfitting phenomena under a small number of data points. If we use a large number of data points, it is expected that the benefit of the inductive bias with deep equilibrium linear models tends to become less noticeable because using large number of data points can reduce the degree of overfitting for all models, including linear models and DNNs. However, repeated the same experiments with all data points of each datasets: for example, we use 60000 training data points and 10000 test data points for MNIST. Figure 8 reports the results where the values are shown with the best learning rates for each model from the set of learning rates SLR = {0.01, 0.005, 0.001, 0.0005, 0.0001, 0.00005) (in terms of the final test errors epoch 100). As can be seen in the figure, deep equilibrium linear models outperformed both linear models and nonlinear DNNs in test errors.
Logistic loss and theoretical bounds. In Corollary 2, we can set Xy = m[l+7)² to get guarantee for the global linear convergence rate for any ini- tialization in theory. However, in practice, this is pessimistic convergence rate and we may want to choose Xr depending on initializations. To demon- strate this, Figure 3 reports the numerical train- ing trajectory along with theoretical upper bounds with initialization Inittalization-nndeeeednnnt X m(1+7)² and initialization -etirradtoon-deeeeeeett >y = inf Amia (D2). As can be seen in Figure 3, the theoretical upper bound with mnttializttion-deeeeeeen Ar demonstrates a faster convergence rate.
training trajectory hound: inet-ind A bound: inst-dopp A. optimal
train loss
of iterations
Figure 3: Logistic loss and theoretical bounds with cewpliiiatiooinneeenddeee Ar and initialization-deeendde År.
33

---

Published as a conference paper at ICLR 2021
(best)
199
DNN (H-
Linear (wors))
13
DWN (H-3)
5022 lOSS
DELM (best)
(best) DELM (wors)
(worst)
ein
200
(a) Modified MNIST: Linear v.s. DELM
(b) Modified MNIST: DNN v.s. DELM
Linear (best)
DNN [H+
Linear (worst)
DNN (H-3)
bo0 toze
DELM (best)
(H-4)
loss
DELM (worst)
tost
DELM (wers)
3
$200
(c) Modified MMISTT Linear v.s. DELM
(d) Modified Fashion-MNIST. DNN v.s. DELM
Linear (best)
DNN (H
Linear (worst)
DAN (H-
DELM (best)
DNN (H-4)
DELM (best)
DELM (worst)
DELM (woret)
1991
3500
0704
1080
283
1991
1200
1290
1066
2003
100
1093
281
550
1931
(c) Modified SEMEION: Linear v.s. DELM
(f) Modified SEMEION: DNN v.s. DELM
Linear (best)
(worst)
DNN (H-
cart (002
DELM (best)
U33 look
DELM (best)
DNN (H--4)
(worst)
1000
2300
.000
1000
2003
2002
00.0
$2.0
2003
2006
43.0
2000
3002
20.6
325
4603
$200
(g) Modified SVHNN: Linear v.s. DELM
(h) Modified SVHN: DNN v.s. DELM
Figure 4: Test and train losses (in log scales) versus the number of epochs for linear models, deep equilibrium linear models (DELMs), and deep neural networks with ReLU (DNNs). The plotted lines indicate the mean values over five random trials whereas the shaded regions represent error bars with one standard deviations for both test and train losses. The plots for linear models and DELMs are shown with the best and worst learning rates (in terms of the final test errors at epoch 5000). The plots for DNNs are shown with the best learning rates for each depth 2, 3, and (in terms of the final test errors at epoch 5000).
34

---

Published as a conference paper at ICLR 2021
DNN-N6 (H-2)
DNN-NB (H-3)
loss
DELM (bes)
DELM (worst)
1990
(a) Modified Kuzushij-MNIST: NN-NB V.S. DELM
DNN NB (H-2)
107
DNN NB (H-2)
361
DNN-NB (1-3)
DNN NB 01-3)
36
toist 10:00
DNN-N (1:4)
DNN-ND 01m4)
Ibest
DELM (worst!
DELM Iworst
39:12
Li0e
4230
2000
$200
(b) Modified CIFAR-I0: DNN-NB v.s. DELM
(c) Modified CIFAR- 100: DNN-NB v.s. DELM
DNN-NB (H-2)
DNN (H-2)
DNN NB (H-3)
DAN-NI 0H-1)
30"
taxst loss
DNN NB (H-4)
Fiin loss
DNN-NB (H-4)
DELM
DELM [woret]
twone)
50
12-4
23
1006
1003
1091
1203
682
1035
0903
200
100
wiee
1999
200
120
023
1991
(d) Modified Fashion-MNIST: DNN-NB v.s. DELM
(e) Modified MNIST: DNN-NB v.s. DELM
(H-2)
(1-2)
DNNANB (H-3)
DRN--N (H-3)
best 0055
DNNINE (H-4)
Lrein 1oss
DNN-NB (H-4)
DELM lbest)
DELM
1006
2006
1003
4830
3220
2306
1802
6822
$2.0
2003
2006
2005
3002
2006
37.0
4403
$2.0
(I) Modified SEMEION: DNN-NB v.s. DELM
(g) Modified SVHN: DNN-NB v.s. DELM
Figure 5: Test and train losses (in log scales) versus the number of epochs for deep equilibrium linear model (DELM) and deep neural network with no bias term (DNN-NB). The plotted lines indicate the mean values over five random trials whereas the shaded regions represent error bars with one standard deviations for both test and train losses. The plots for DELM are shown with the best and worst learning rates (in terms of the final test errors at epoch 5000). The plots for DNN-NB are shown with the best learning rates for each depth 2, 3, and 4 (in terms of the final test errors at epoch = 5000).
35

---

Published as a conference paper at ICLR 2021
DNN-N6 (H-10)
161
(H-100)
Lest loss
train l055
10
DELM (besi)
DNN NB (H:-200)
30
2000
4000
5007
3000
2006
3902
6003
3e00
spect
spach
(a) Modified Kuzushiji-Min MNIST: DNN-NB v.S. DELM
best look
Urain loss
Dost
2003
1806
1022
1230
1156
2003
2000
2002
286
4103
5230
(b) Modified CIFAR-10: DNN-NB v.s. DELM
(c) Modified CIFAR- 100: DNN-NB v.s. DELM
36"
text loss
Erain loss
Drain
L0C8
2000
2006
4000
2008
3000
200
3006
4669
$900
100
2006
2099
4900
2008
3000
2006
200
4043
5200
(d) Modified Fashion-MNiBn DNN-NB v.s. DELM
(e) Modified MNIST: DNN-NB v.s. DELM
best lo8s
Drain losS
107
1008
2000
2000
4000
5000
3090
300
3000
4009
5900
1008
2000
3000
4300
3009
3800
2006
3300
4009
5200
(f) Modified SEMEION: DNN-NB v.s. DELM
(g) Modified SVHN: DNN-NB v.s. DELM
Figure 6: Test and train losses (in log scales) versus the number of epochs for deep equilibrium linear model (DELM) and deeper neural network with no bias term (DNN-NB). The legend the same for all subplots and shown in subplot (a). The plotted lines indicate the mean values over five random trials whereas the shaded regions represent error bars with one standard deviations for both test and train losses. The plots for DELM are shown with the best and worst learning rates (in terms of the final test errors at epoch 5000). The plots DNN-NB are shown with the best learning rates for each depth H 10, 100, and 200 (in terms of the final test errors at epoch = 5000). The values for DNNs with depth H 100 and 200 coincide in the plots.
36

---

Published as a conference paper at ICLR 2021
18"
DNN (H-10)
(H=100)
test loss
train loss
10:7-
DNN (H-200)
(best)
20%
DELM (worst)
1006
2006
1000
6000
5992
1099
2010
1000
0000
502e
spach
(a) Modified Kuzushiji-MNiS MNNIT:: DNN v.s. DELM
307
Lest Lestloss
1
Urain 105
Lest lood
1pch
60127
(b) Modified CIFAR-10: DNN v.s. DELM
(c) Modified CIFAR-100: DNN v.s. DELM
307
28"
best loss
Lrain 1O0S
best
55
2000
$300
2000
2003
4890
1009
3006
4000
5350
2004
300
2000
(d) Modified Fashion -ANIIT:: DNN v.s. DELM
(c) Modified MNIST: DNN v.s. DELM
1
test 1
10
Lest loss
Trein 1000
1283
2003
63c
1330
201
1200
2000
V000
188
1203
120
1293
12:0
4940
V00
(1) Modified SEMEION: DNN v.s. DELM
(g) Modified SVHN: DNN v.s. DELM
Figure 7: Test and train losses (in log scales) versus the number of epochs for deep equilibrium linear model (DELM) and deeper neural networks with ReLU (DNNss. The legend is the same for all subplots and shown in subplot (a). The plotted lines indicate the mean values over five random trials whereas the shaded regions represent train losses. The plots for DELM are shown with the best and worst learning rates (in terms of the final errors epoch 5000). The plots for DNNs are shown with the best learning rates for each depth H = 10, 100, and 200 (in terms of the final test errors at epoch 5000). The values for DNNs with depth 100 and 200 coincide in the plots.
bars with standard deviations for both test and
37

---

Published as a conference paper at ICLR 2021
DNN (H==2)
18
13"
test loss
DNN (H-10)
20-1-
train
Linear
16-9
28-4
DELM
23
(a) Modified Kuzushiji-MNIST: NN-NB DELM
teet loss
train lOSS
12
30-
12
16-4
6002
1000
(b) Modified CIFAR-10: DNN-NB v.s. DELM
(c) Modified CIFAR-100: DNN-NB v.s. DELM
36"
30"
10*
bost loss
10.
train los
20-4
40
12.0
10:0
180
390
(d) Modified Fashion-MNIST: DNN-NB v.s. DELM
(e) Modified MNIST: DNN-NB v.s. DELM
bentt loos
train loss
5.
15
20-
62 12a
28
126
40
(l) Modified SEMEION: DNN-NB v.S. DELM
(g) Modified SVHN: DNN-NB v.s. DELM
Figure 8: Test and train losses (in log scales) versus the number of epochs for linear models, deep equilibrium linear models (DELMss), and deep neural networks with ReLU (DNNs). The legend the same for all subplots and shown in subplot (a). The plotted lines indicate the mean values over three random trials whereas the shaded regions represent error bars with one standard deviations.
38


Reviews:
[{"review_id": "qkiR3cyajSb", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper discusses the theory of deep equilibrium models with linear activations. The model weights are softmaxed to ensure that inference converges to a fixed point, a necessary condition for training deep equilibrium models. The paper then analyzes the gradient flow dynamics of such models. The main result is that linear-rate convergence is guaranteed for a class of loss functions, including quadratic and logistic losses, when training with gradient flow. This conclusion is supported by experiments conducted in a teacher-student-like setup, where the labels are generated by a teacher deep equilibrium model, showing that training does converge in practice.\n\nDeep equilibrium models represent a novel way to train neural networks, and not much is known about them yet theoretically. It is important that we understand the dynamics of such models better, and this paper is a good step in that direction.\n\nSuggested improvements:\n\n1. Definition 1 has a typo: On the right-hand side of the inequality there shouldn't be a gradient.\n\n2. The main results are presented clearly, and the paper is generally easy to read. The only exception is section 3.2 on the connection with trust region Newton methods, which I did not understand. I recommend clarifying the intuition behind Theorem 2, as well as the main message of this section.\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "qXeYyPfOWAU", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "> Summary: This work focuses on the study of (global) convergence and gradient dynamics of a recently proposed family of models, the deep equilibrium  (linear) models (DELM) under common classes of loss functions. Exploiting the Neumann series convergence and the PL inequality analysis, the authors proved convergence to global optima of DELM without prior assumption on the width m of the model (relative to the number of data n). \n\nHere is my general opinion:\n\nWhile deep linear models in general (as the authors acknowledged) has been widely studied, and despite the extreme simplicity of the DELM when compared to the original DEQ model that Bai et al. [1] studied, I found this work *interesting and important* in establishing a solid foundation for the theoretical study of this class of implicit-depth models. The authors managed to demonstrate that the gradient dynamics and convergence assumptions of DELM is __indeed__ different from typical \"stacked\" deep linear models that prior approaches study, such as deep linear ResNet. And throughout the arguments of the paper and the proof in the appendix, I can tell how the \"equilibrium\" property of DELM is making the story different, and think this paper sets a good starting point for future similar in this direction for general implicit-depth models. But still, the paper has a limited scope in terms of the structure it studies.\n\nPros:\n1. One of the first theoretical works on the gradient dynamics and convergence properties of the deep equilibrium models [1] (and implicit models in general [2,3]), which are quite different from conventional deep networks.\n2. Clear notation and theoretical insights, with proof relatively easy to follow. The proof seems overall correct (there are some that I didn't check closely though).\n3. Clear discussion of the relation, including (and especially) the differences of the prior analysis on deep linear neural networks.\n\nCons:\n1. The very definition of DELM, which the author provided a particular formulation of, is of limited scope (see my comment below that expands on this point).\n2. The empirical studies to validate the conclusion of the theoretical results could be strengthened.\n\n-------------------------------\n\nI have some comments/questions for the authors, detailed below:\n\n1. The major limitation that I found while reading the discussion and the proof of this paper is that while the authors claim to study deep equilibrium **linear** models, the insights mostly only apply to the models converging with Neumann series guarantee, and can be written in the form $BU^{-1} \\phi(x)$. I understand the motivation for fixating a provably convergent equilibrium model formulation. But:\n\n    a) A provably convergent deep equilibrium linear model doesn't need to be Neumann (for its Jacobian) for the implicit function theorem to work. In the simplest case, the fixed point of a function $h(x)$ on 2D can have a local derivative with absolute value > 1. This certainly implies that repeatedly unrolling the function $h(x)$ may not converge, and yet **there still is** a unique fixed point and one can reliably solve for it (see [4]). However, without the nice Neumann series form, which allows one to write $(I-\\gamma \\sigma(A))^{-1}$ as a closed-form representation for the \"infinite-depth\" network forward pass, I don't think the theorems will hold directly. Typically, the $(I-J)^{-1}$ term should only appear in the implicit function theorem, which is used for the backward pass. I expect the authors to clarify this further.\n\n    b) The very design of $\\sigma(A)$ in the model the authors study is a bit bizarre to me. Why applying a \"softmax\" on the weight? Is it just to ensure that proposition 1 holds (i.e., that you have a handy, provably-convergent linear model)? The authors stressed a few times that the $h(\\cdot)$ function is thus \"non-linear\" w.r.t. $\\sigma(A)$, but I fail to directly see why it matters so much as the model is still linear w.r.t. the input (it's really a one-linear layer; though the inverse from Neumann does make a difference on A), and in terms of the gradient dynamics, the major difference this makes will merely be $\\frac{\\partial \\sigma(A)}{\\partial A}$. I might have missed something here and would appreciate if the authors can clarify.\n\n2. I didn't quite get the specific point the authors were trying to make in Section 3.2 in terms of the implicit bias. Could you expand on that?\n\n3. Overall, I feel that the empirical support of the theoretical findings can be stronger, for instance, by inspecting different initialization $(A_0, B_0)$, or validating the radius discussion at the end of Section 3.1 for the logistic loss. Like in Zou et al. [5], some synthetic data could probably work just fine. What is the reason for only using 200 images from the MNIST/CIFAR datasets? Is it to keep the size of $\\Phi$ small (but I didn't see the authors report anything about it in Section 2.2). And since the primary purpose of Sec. 2.2. is to \"discuss whether the model would also make sense in practice\" (which I take to mean that you only want to compare the test accuracies coming out of these models), wouldn't a 200-sample version of MNIST/CIFAR too small to draw a robust conclusion on this?\n\n4. One that I think could be useful for further thought is the convergence property not just for GD, but also SGD, like Zou et al. provided in [5] (up to a probability).\n\n-------------------------------\n\nMinor things:\n\n  - i) Page 3: \"the outputs of the deep equilibrium linear models $f_\\theta(x)= \\cdots$ are nonlinear and non-multilinear in the optimization variables $(A,B)$.\" Non-linear in even $B$?\n  - ii) Page 14: $V_{q*}$ --> $B_{q*}$\n  - iii) Page 16: $\\nabla_F L(A,B)$ --> $\\nabla_A L(A,B)$\n\n-------------------------------\n\n- [1] https://arxiv.org/abs/1909.01377\n- [2] https://arxiv.org/abs/1908.06315\n- [3] https://arxiv.org/abs/2009.06211\n- [4] https://arxiv.org/abs/2006.08591\n- [5] https://arxiv.org/abs/2003.01094", "extra_scores": {}, "extra_text": {}}, {"review_id": "JKthYvQ0jlm", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This submission studies the dynamics and convergence properties of \"deep equilibrium models\", which are parametric fixed-point iterations corresponding to the infinite depth limit of \"weight-tied\" neural networks. As the authors point out, these networks differ from deep linear networks and networks in the NTK scaling in that the optimization remains nonlinear w/r/t the parameters. The authors prove two results: first, they establish linear convergence to the global minimum under the relatively strict assumption of a \"local\" PL-inequality; secondly, they show that the dynamics of the deep equilibrium models differs from gradient descent dynamics and, in fact, is related to a trust region Newton method.\n\nThe first theorem is nice and well-presented, but I think the issue of the radius over which the PL-inequality holds could be better discussed. I could not tell whether or not the convergence depends on starting within the locally smooth and quadratically bounded region of the loss. \n\nThe second theorem, regarding the nature of the dynamics, lacks a clear interpretation. Basically, all that is said is that the dynamics is distinct from what would be seen in a linear model, and it's not clear that this dynamics has anything to do with implicit regularization, as the authors suggest. I would recommend the authors clarify the discussion of this result. \n\nThe experiments are somewhat bizarre and I felt that they were a little misleading about the representative power / potential of these models, but perhaps I did not fully understand the set-up and intent. The authors randomly sample a deep equilibrium model and then use it to represent a conditional probability distribution where the conditioning is with image data. Then, the authors fit various models to this distribution and show that the deep equilibrium models (which is precisely the underlying function representing the distribution) has better performance than other classes of functions. I think the description in this section could be vastly improved, perhaps presenting this more akin to a student-teacher problem. \n\n ", "extra_scores": {}, "extra_text": {}}, {"review_id": "DyWSOpWjXTw", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "**Overview**\n\nThis paper purports to study training deep equilibrium models by studying the optimization dynamics of deep *linear* equilibrium models. The original deep equilibrium model is formalized as follows: Given a training dataset $(x_i,y_i)$ for i = 1,...,n where $x_i \\in \\mathcal{X}\\subseteq \\mathbb{R}^{m_x}$ and $y_i\\in \\mathcal{Y}\\subseteq \\mathbb{R}^{m_y}$ are the $i$-th input and output, respectively. The goal is to learn a predictor from a family $\\mathcal{H} = \\\\{\\{ f_\\theta : \\mathbb{R}^{m_x} \\rightarrow \\mathbb{R}^{m_y} | \\theta\\in\\Theta\\}\\\\}$. Then, instead of trying to map $x$ to $y$ using finite amount of layers, deep equilibrium models assume infinite number of layers, and the output $z^*$ of the last hidden layer is defined by\n\\begin{equation}\n    z^* = \\lim_{l\\rightarrow \\infty} z^{(l)} = \\lim_{l\\rightarrow \\infty} h(z^{(l-1)};x,\\theta) = h(z^*;x,\\theta)\n\\end{equation}\nwhere $h$ is some continuous function of choice. \n\nIn particular with deep equilibrium **linear** models, $h$ is constrained as follows:\n\\begin{equation}\n    h(z^{(l-1)};x,\\theta) = \\gamma\\sigma(A)z^{(l-1)}+\\phi(x)\n\\end{equation}\nwhere $\\phi(x)$ is a feature map of $x$ and transforms $x\\in\\mathbb{R}^{m_x}$ into $\\phi(x)\\in\\mathbb{R}^m$. $\\theta = (A,B)$ are two trainable matrices, where $A\\in\\mathbb{R}^{m\\times m}$ is for computing each hidden output and $B\\in\\mathbb{R}^{m_y\\times m}$ is for computing the final output of the network. $\\gamma \\in (0,1)$ is some positive real number, and $\\sigma$ is a nonlinear function to ensure the existence of the fixed point $z^*$. This model is linear in z. The objective function of this deep equilibrium linear model can be written as follows:\n\\begin{equation}\n    L(A,B) = \\sum_{i=1}^n \\ell\\left( B\\left( \\lim_{l\\rightarrow\\infty}z^{(l)}(x_i,A)\\right),y_i\\right)\n\\end{equation}\nwhere $\\ell$ is some choice of loss function.\n\nThen this paper provides some motivations behind studying the dynamics of these deep equilibrium linear models by presenting some interesting comparisons between deep equilibrium linear models and \"normal\" linear models and additional, less interesting comparisons with fully-connected feed-forward deep neural networks (FNN), using standard image datasets CIFAR-10, CIFAR-100 and Kuzushiji-MNIST. In their tests, the deep equilibrium linear models outperformed both linear models and FNNs. \n\nThe main results from this paper is uncovering the dynamics behind these deep equilibrium linear models. This paper provided a sequence of proofs that shows linear convergence of these models step by step, under the assumption that the loss functions $\\ell$ are differentiable, which is satisfied by some standard loss functions, such as square loss, losgistic loss, and smoothed hinge loss $\\ell(f_\\theta(x),y) = (\\max\\{0,1-f_\\theta(x)y\\})^k$ with $k\\geq 2$. Also, despite of non-convexity of the loss functions, it is shown in this paper that a global minimum $L^*$ always exists and under their assumptions, the deep equilibrium linear models will always converge to the global minimum linearly.\n\n**Strengths**\n\nThe main strength of this paper is that it brought  some interesting insights into deep equilibrium  models, and they showed their results rigorously. The definitions and propositions are clear enough for readers with some analysis and machine learning background to fully understand. \n\nSince the dynamics of deep learning models is an open field of research and isn't discovered fully, this paper will definitely contribute to the deep learning field in understanding the dynamics and convergence theory of these deep equilibrium linear models, and potentially benefiting researches on understanding more general deep learning models.\n\n**Weakness**\n\nAlthough this paper brought some nice ideas, some of the methodologies are not quite convincing. For example, in the experiments shown in this paper, they  compare performance of deep equilibrium linear models with linear models and deep neural networks. Especially in the comparison with the deep networks, it doesn't seem that the networks are deep enough for the comparison to be compelling. Also in the same experiments, they assumed the true data distribution is approximately given by a deep equilibrium linear model and generated data according to this model -- i.e. the data is only semi-\"real\". It would be better to show that the deep equilibrium linear model outperforms other models under a more general setting.\n\nThroughout the entire paper, it's unclear what are the contributions. To be specific, in the first two sections, it's unclear whether they want to show the trainability of the deep equilibrium linear models or learn the dynamics of these models. Also, some of the interesting aspects of this paper were missing some details. For example, one would be interested in seeing why exactly can the deep equilibrium linear models outperform linear models, since they are both linear and the only difference is how they are trained. In the last two sections of this paper, they also brought up *implicit bias*, which would be another interesting topic to dive into. It would be great seeing more comparisons on those aspects.", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper analyzes the gradient flow dynamics of deep equilibrium models with linear activations and establishes linear convergence for quadratic loss and logistic loss; several exciting results and connections, solid contribution, accept! 
Decision: Accept (Spotlight)