Title: PROVABLE BENEFIT OF ORTHOGONAL INITIALIZA-TION IN OPTIMIZING DEEP LINEAR NETWORKS

Abstract: The selection of initial parameter values for gradient-based optimization of deep neural networks is one of the most impactful hyperparameter choices in deep learning systems, affecting both convergence times and model performance. Yet despite significant empirical and theoretical analysis, relatively little has been proved about the concrete effects of different initialization schemes. In this work, we analyze the effect of initialization in deep linear networks, and provide for the first time a rigorous proof that drawing the initial weights from the orthogonal group speeds up convergence relative to the standard Gaussian initialization with iid weights. We show that for deep networks, the width needed for efficient convergence to a global minimum with orthogonal initializations is independent of the depth, whereas the width needed for efficient convergence with Gaussian initializations scales linearly in the depth. Our results demonstrate how the benefits of a good initialization can persist throughout learning, suggesting an explanation for the recent empirical successes found by initializing very deep non-linear networks according to the principle of dynamical isometry.

1. INTRODUCTION:

Through their myriad successful applications across a wide range of disciplines, it is now well established that deep neural networks possess an unprecedented ability to model complex real-world datasets, and in many cases they can do so with minimal overfitting. Indeed, the list of practical achievements of deep learning has grown at an astonishing rate, and includes models capable of human-level performance in tasks such as image recognition (Krizhevsky et al., 2012) , speech recognition (Hinton et al., 2012) , and machine translation (Wu et al., 2016 ).

Yet to each of these deep learning triumphs corresponds a large engineering effort to produce such a high-performing model. Part of the practical difficulty in designing good models stems from a proliferation of hyperparameters and a poor understanding of the general guidelines for their selection. Given a candidate network architecture, some of the most impactful hyperparameters are those governing the choice of the model's initial weights. Although considerable study has been devoted to the selection of initial weights, relatively little has been proved about how these choices affect important quantities such as rate of convergence of gradient descent.

In this work, we examine the effect of initialization on the rate of convergence of gradient descent in deep linear networks. We provide for the first time a rigorous proof that drawing the initial weights from the orthogonal group speeds up convergence relative to the standard Gaussian initialization with iid weights. In particular, we show that for deep networks, the width needed for efficient convergence for orthogonal initializations is independent of the depth, whereas the width needed for efficient convergence of Gaussian networks scales linearly in the depth.

Orthogonal weight initializations have been the subject of a significant amount of prior theoretical and empirical investigation. For example, in a line of work focusing on dynamical isometry, it was found that orthogonal weights can speed up convergence for deep linear networks (Saxe et al., 2014; Advani & Saxe, 2017) and for deep non-linear networks (Pennington et al., 2018; Xiao et al., 2018; Gilboa et al., 2019; Chen et al., 2018; Pennington et al., 2017; Tarnowski et al., 2019; Ling & Qiu, 2019) when they operate in the linear regime. In the context of recurrent neural networks, orthogonality can help improve the system's stability. A main limitation of prior work is that it has focused almost exclusively on model's properties at initialization. In contrast, our analysis focuses on the benefit of orthogonal initialization on the entire training process, thereby establishing a provable benefit for optimization.

The paper is organized as follows. After reviewing related work in Section 2 and establishing some preliminaries in Section 3, we present our main positive result on efficient convergence from orthogonal initialization in Section 4. In Section 5, we show that Gaussian initialization leads to exponentially long convergence time if the width is too small compared with the depth. In Section 6, we perform experiments to support our theoretical results.

2. RELATED WORK:

Deep linear networks. Despite the simplicity of their input-output maps, deep linear networks define high-dimensional non-convex optimization landscapes whose properties closely reflect those of their non-linear counterparts. For this reason, deep linear networks have been the subject of extensive theoretical analysis. A line of work (Kawaguchi, 2016; Hardt & Ma, 2016; Lu & Kawaguchi, 2017; Yun et al., 2017; Zhou & Liang, 2018; Laurent & von Brecht, 2018) studied the landscape properties of deep linear networks. Although it was established that all local minima are global under certain assumptions, these properties alone are still not sufficient to guarantee global convergence or to provide a concrete rate of convergence for gradient-based optimization algorithms.

Another line of work directly analyzed the trajectory taken by gradient descent and established conditions that guarantee convergence to global minimum (Bartlett et al., 2018; Arora et al., 2018; Du & Hu, 2019) . Most relevant to our work is the result of Du & Hu (2019) , which shows that if the width of hidden layers is larger than the depth, gradient descent with Gaussian initialization can efficiently converge to a global minimum. Our result establishes that for Gaussian initialization, this linear dependence between width and depth is necessary, while for orthogonal initialization, the width can be independent of depth. Our negative result for Gaussian initialization also significantly generalizes the result of Shamir (2018) , who proved a similar negative result for 1-dimensional linear networks.

Orthogonal weight initializations. Orthogonal weight initializations have also found significant success in non-linear networks. In the context of feedforward models, the spectral properties of a network's input-output Jacobian have been empirically linked to convergence speed (Saxe et al., 2014; Pennington et al., 2017; 2018; Xiao et al., 2018) . It was found that when this spectrum concentrates around 1 at initialization, a property dubbed dynamical isometry, convergence times improved by orders of magnitude. The conditions for attaining dynamical isometry in the infinitewidth limit were established by Pennington et al. (2017; 2018) and basically require that input-output map to be approximately linear and for the weight matrices to be orthogonal. Therefore the training time benefits of dynamical isometry are likely rooted in the benefits of orthogonality for deep linear networks, which we establish in this work.

Orthogonal matrices are also frequently used in the context of recurrent neural networks, for which the stability of the state-to-state transition operator is determined by the spectrum of its Jacobian (Haber & Ruthotto, 2017; Laurent & von Brecht, 2016) . Orthogonal matrices can improve the conditioning, leading to an ability to learn over long time horizons (Le et al., 2015; Henaff et al., 2016; Chen et al., 2018; Gilboa et al., 2019) . While the benefits of orthogonality can be quite large at initialization, little is known about whether or in what contexts these benefits persist during training, a scenario that has lead to the development of efficient methods of constraining the optimization to the orthogonal group (Wisdom et al., 2016; Vorontsov et al., 2017; Mhammedi et al., 2017) . Although we do not study the recurrent setting in this work, an extension of our analysis might help determine when orthogonality is beneficial in that setting.

3. PRELIMINARIES:

3.1 NOTATION Let [n] = {1, 2, . . . , n}. Denote by I d the d × d identity matrix, and by I an identity matrix when its dimension is clear from context. Denote by N (µ, σ 2 ) the Gaussian distribution with mean µ and variance σ 2 , and by χ 2 k the chi-squared distribution with k degrees of freedom. Denote by • the 2 norm of a vector or the spectral norm of a matrix. Denote by • F the Frobenius norm of a matrix. For a symmetric matrix A, let λ max (A) and λ min (A) be its maximum and minimum eigenvalues, and let λ i (A) be its i-th largest eigenvalue. For a matrix B ∈ R m×n , let σ i (B) be its i-th largest singular value (i = 1, 2, . . . , min{m, n}), and let σ max (B) = σ 1 (B), σ min (B) = σ min{m,n} (B). Denote by vec (A) be the vectorization of a matrix A in column-first order. The Kronecker product between two matrices A ∈ R m1×n1 and B ∈ R m2×n2 is defined as

A ⊗ B =    a 1,1 B • • • a 1,n1 B . . . . . . . . . a m1,1 B • • • a m1,n1 B    ∈ R m1m2×n1n2 ,

where a i,j is the element in the (i, j)-th entry of A.

We use the standard O(•), Ω(•) and Θ(•) notation to hide universal constant factors. We also use C to represent a sufficiently large universal constant whose specific value can differ from line to line.

3.2. PROBLEM SETUP:

Suppose that there are

n training examples {(x k , y k )} n k=1 ⊂ R dx × R dy .

Denote by X = (x 1 , . . . , x n ) ∈ R dx×n the input data matrix and by Y = (y 1 , . . . , y n ) ∈ R dy×n the target matrix. Consider an L-layer linear neural network with weight matrices W 1 , . . . , W L , which given an input x ∈ R dx computes

EQUATION

where

W i ∈ R di×di-1 (i = 1, . . . , L), d 0 = d x , d L = d y ,

and α is a normalization constant which will be specified later according to the initialization scheme. We study the problem of training the deep linear network by minimizing the 2 loss over training data:

EQUATION

The algorithm we consider to minimize the objective (2) is gradient descent with random initialization, which first randomly samples the initial weight matrices {W i (0)} L i=1 from a certain distribution, and then updates the weights using gradient descent: for time t = 0, 1, 2, . . .,

EQUATION

where η > 0 is the learning rate.

For convenience, we denote

W j:i = W j W j-1 • • • W i (1 ≤ i ≤ j ≤ L) and W i-1:i = I (i ∈ [L]).

The time index t is used on any variable that depends on W 1 , . . . , W L to represent its value at time t, e.g.,

W j:i (t) = W j (t) • • • W i (t), (t) = (W 1 (t), . . . , W L (t)), etc.

4. EFFICIENT CONVERGENCE USING ORTHOGONAL INITIALIZATION:

In this section we present our main positive result for orthogonal initialization. We show that orthogonal initialization enables efficient convergence of gradient descent to a global minimum provided that the hidden width is not too small.

In order to properly define orthogonal weights, we let the widths of all hidden layers be equal:

d 1 = d 2 = • • • = d L-1 = m, and let m ≥ max{d x , d y }. Note that all intermediate matrices W 2 , . . . , W L-1 are m × m square matrices, and W 1 ∈ R m×dx , W L ∈ R dy×m .

We sample each initial weight matrix W i (0) independently from a uniform distribution over scaled orthogonal matrices satisfying

EQUATION

In accordance with such initialization, the scaling factor α in ( 1) is set as

α = 1 √ m L-1 dy , which ensures E f (x; W L (0), . . . , W 1 (0)) 2 = x 2 for any x ∈ R dx . 1

The same scaling factor was adopted in Du & Hu (2019) , which preserves the expectation of the squared 2 norm of any input.

Let W * ∈ arg min W ∈R dy ×dx W X -Y F and * = 1 2 W * X -Y 2 F .

Then * is the minimum value for the objective (2). Denote r = rank(X), κ = λmax(X X) λr(X X) , and r =

X 2 F X 2 . 2

Our main theorem in this section is the following:

EQUATION

for some δ ∈ (0, 1) and a sufficiently large universal constant C > 0. Set the learning rate η ≤ dy 2L X2 . Then with probability at least 1 -δ over the random initialization, we have

(0) - * ≤ O 1 + log(r/δ) d y + W * 2 X 2 F , (t) - * ≤ 1 - 1 2 ηLλ r (X X)/d y t ( (0) - * ), t = 0, 1, 2, . . . ,

where (t) is the objective value at iteration t.

Notably, in Theorem 4.1, the width m need not depend on the depth L. This is in sharp contrast with the result of Du & Hu (2019) for Gaussian initialization, which requires m ≥ Ω(Lrκ 3 d y ). It turns out that a near-linear dependence between m and L is necessary for Gaussian initialization to have efficient convergence, as we will show in Section 5. Therefore the requirement in Du & Hu (2019) is nearly tight in terms of the dependence on L. These results together rigorously establish the benefit of orthogonal initialization in optimizing very deep linear networks.

If we set the learning rate optimally according to Theorem 4.1 to η = Θ( dy L X 2 ), we obtain that (t) - * decreases by a ratio of 1 -Θ(κ -1 ) after every iteration. This matches the convergence rate of gradient descent on the (1-layer) linear regression problem min

W ∈R dy ×dx 1 2 W X -Y 2 F .

4.1. PROOF OF THEOREM 4.1:

The proof uses the high-level framework from Du & Hu (2019) , which tracks the evolution of the network's output during optimization. This evolution is closely related to a time-varying positive semidefinite (PSD) matrix (defined in ( 7)), and the proof relies on carefully upper and lower bounding the eigenvalues of this matrix throughout training, which in turn implies the desired convergence result.

First, we can make the following simplifying assumption without loss of generality. See Appendix B in Du & Hu (2019) for justification. Assumption 4.1. (Without loss of generality) X ∈ R dx×r , rank(X) = r, Y = W * X, and * = 0. Now we briefly review Du & Hu (2019)'s framework. The key idea is to look at the network's output, defined as

U = αW L:1 X ∈ R dy×n .

We also write U (t) = αW L:1 (t)X as the output at time t. Note that (t) = 1 2 U (t) -Y 2 F . According to the gradient descent update rule, we write

W L:1 (t + 1) = i W i (t) -η ∂ ∂W i (t) = W L:1 (t) - L i=1 ηW L:i+1 (t) ∂ ∂W i (t)W i-1:1 (t) + E(t),

where E(t) contains all the high-order terms (i.e., those with η 2 or higher). With this definition, the evolution of U (t) can be written as the following equation:

EQUATION

where

EQUATION

Notice that P (t) is always PSD since it is the sum of L PSD matrices. Therefore, in order to establish convergence, we only need to (i) show that the higher-order term E(t) is small and (ii) prove upper and lower bounds on P (t)'s eigenvalues. For the second task, it suffices to control the singular values of W i-1:1 (t) and W L:i+1 (t) (i ∈ [L]).3 Under orthogonal initialization, these matrices are perfectly isometric at initialization, and we will show that they stay close to isometry during training, thus enabling efficient convergence.

The following lemma summarizes some properties at initialization. Lemma 4.2. At initialization, we have

EQUATION

Furthermore, with probability at least 1 -δ, the loss at initialization satisfies

EQUATION

Proof sketch. The spectral property (8) follows directly from (4).

To prove (9), we essentially need to upper bound the magnitude of the network's initial output. This turns out to be equivalent to studying the magnitude of the projection of a vector onto a random lowdimensional subspace, which we can bound using standard concentration inequalities. The details are given in Appendix A.1. Now we proceed to prove Theorem 4.1. We define

B = O 1 + log(r/δ) dy + W * 2 X 2

F which is the upper bound on (0) from ( 9). Conditioned on (9) being satisfied, we will use induction on t to prove the following three properties A(t), B(t) and C(t) for all t = 0, 1, . . .:

• A(t): (t) ≤ 1 -1 2 ηLσ 2 min (X)/d y t (0) ≤ 1 -1 2 ηLσ 2 min (X)/d y t B. • B(t): σ max (W j:i (t)) ≤ 1.1m j-i+1 2 , σ min (W j:i (t)) ≥ 0.9m j-i+1 2 , ∀1 ≤ i ≤ j ≤ L, (i, j) = (1, L). • C(t): W i (t) -W i (0) F ≤ 8 √ Bdy X Lσ 2 min (X) , ∀1 ≤ i ≤ L.

A(0) and B(0) are true according to Lemma 4.2, and C(0) is trivially true. In order to prove A(t), B(t) and C(t) for all t, we will prove the following claims for all t ≥ 0:

Claim 4.3. A(0), . . . , A(t), B(0), . . . , B(t) =⇒ C(t + 1). Claim 4.4. C(t) =⇒ B(t). Claim 4.5. A(t), B(t) =⇒ A(t + 1).

The proofs of these claims are given in Appendix A. Notice that we finish the proof of Theorem 4.1 once we prove A(t) for all t ≥ 0.

5. EXPONENTIAL CURSE OF GAUSSIAN INITIALIZATION:

In this section, we show that gradient descent with Gaussian random initialization necessarily suffers from a running time that scales exponentially with the depth of the network, unless the width becomes nearly linear in the depth. Since we mostly focus on the dependence of width and running time on depth, we will assume the depth L to be sufficiently large.

Recall that we want to minimize the objective

EQUATION

F by gradient descent. We assume Y = W * X for some W * ∈ R dy×dx , so that the optimal objective value is 0. For convenience, we assume X F = Θ(1) and Y F = Θ(1).

Suppose that at layer:

i ∈ [L], every entry of W i (0) is sampled from N (0, σ 2 i )

, and all weights in the network are independent. We set the scaling factor α such that the initial output of the network does not blow up exponentially (in expectation): 1) .

EQUATION

) Note that E f (x; W 1 (0), . . . , W L (0)) 2 = α 2 L i=1 (d i σ 2 i ) x 2 . Thus (10) means α 2 L i=1 (d i σ 2 i ) ≤ L O(

We also assume that the magnitude of initialization at each layer cannot vanish with depth:

EQUATION

Note that the assumptions (10) and ( 11) are just sanity checks to rule out the obvious pathological cases -they are easily satisfied by all the commonly used initialization schemes in practice.

Now we formally state our main theorem in this section. Theorem 5.1. Suppose max{d 0 , d 1 , . . . , d L } ≤ O(L 1-γ ) for some universal constant 0 < γ ≤ 1.

Then there exists a universal constant c > 0 such that, if gradient descent is run with learning rate η ≤ e cL γ , then with probability at least 0.9 over the random initialization, for the first e Ω(L γ ) iterations, the objective value is stuck between 0.4

Y 2 F and 0.6 Y 2 F .

Theorem 5.1 establishes that efficient convergence from Gaussian initialization is impossible for large depth unless the width becomes nearly linear in depth. This nearly linear dependence is the best we can hope for, since Du & Hu (2019) proved a positive result when the width is larger than linear in depth. Therefore, a phase transition from untrainable to trainable happens at the point when the width and depth has a nearly linear relation. Furthermore, Theorem 5.1 generalizes the result of Shamir (2018) , which only treats the special case of

d 0 = • • • = d L = 1.

5.1 PROOF OF THEOREM 5.1

For convenience, we define a scaled version of W i : let 1) and αW L:1 = βA L:1 , where

A i = W i /( √ d i σ i ) and β = α L i=1 ( √ d i σ i ). Then we know β ≤ L O(

A j:i = A j • • • A i .

We first give a simple upper bound on A j:i (0) for all 1 ≤ i ≤ j ≤ L.

Lemma 5.2. With probability at least 1 -δ, we have

A j:i (0) ≤ O L 3 δ for all 1 ≤ i ≤ j ≤ L.

The proof of Lemma 5.2 is given in Appendix B.1. It simply uses Markov inequality and union bound.

Furthermore, a key property at initialization is that if j -i is large enough, A j:i (0) will become exponentially small. Lemma 5.3. With probability at least 1 -e -Ω(L γ ) , for all 1 ≤ i ≤ j ≤ L such that j -i ≥ L 10 , we have A j:i (0) ≤ e -Ω(L γ ) .

Proof. We first consider a fixed pair (i, j) such that j -i ≥ L 10 . In order to bound A j:i (0) , we first take an arbitrary unit vector v ∈ R di-1 and bound A j:i (0)v . We can write

A j:i (0)v 2 = j k=i Z k , where Z k = A k:i (0)v 2 A k-1:i (0)v 2 . Note that for any nonzero v ∈ R d k-1 independent of A k (0), the distribution of d k • A k (0)v 2 v 2 is χ 2 d k .

Therefore, Z i , . . . , Z j are independent, and

d k Z k ∼ χ 2 d k (k = i, i + 1, . . . , j).

Recall the expression for the moments of chi-squared random variables:

E Z λ k = 2 λ Γ(d k /2+λ) d λ k Γ(d k /2) (∀λ > 0)

. Taking λ = 1 2 and using the bound

Γ(a+ 1 2 ) Γ(a) ≤ √ a -0.1 (∀a ≥ 1 2 ) (Qi & Luo, 2012), we get E √ Z k ≤ 2(d k /2-0.1) d k = 1 -0.2 d k ≤ 1 -0.1 d k . Therefore we have E j k=i Z k ≤ j k=i 1 - 0.1 d k ≤ 1 - 0.1 O(L 1-γ ) j-i+1 ≤ 1 -Ω(L γ-1 ) L 10 = e -Ω(L γ ) .

Choose a sufficiently small constant c > 0.

By Markov inequality we have

Pr j k=i Z k > e -c L γ ≤ e c L γ E j k=i Z k ≤ e c L γ e -Ω(L γ ) = e -Ω(L γ )

. Therefore we have shown that for any fixed unit vector v ∈ R di-1 , with probability at least 1 -e -Ω(L γ ) we have

A j:i (0)v ≤ e -Ω(L γ ) .

Next, we use this to bound A j:i (0) via an -net argument. We partition the index set

[d i-1 ] into [d i-1 ] = S 1 ∪ S 2 ∪ • • • ∪ S q such that |S l | ≤ L γ/2 (∀l ∈ [q]

) and q = O( di-1 L γ/2 ). For each l ∈ [q], let N l be a 1 2 -net for all the unit vectors in R di-1 with support in S l . Note that we can choose N l such that |N l | = e O(|S l |) = e O(L γ/2 ) . Taking a union bound over ∪ q l=1 N l , we know that A j:i (0)v ≤ e -Ω(L γ ) v simultaneously for all v ∈ ∪ q l=1 N l with probability at least 1 -(

q l=1 |N l |) e -Ω(L γ ) ≥ 1 -q • e O(L γ/

2 ) e -Ω(L γ ) = 1 -e -Ω(L γ ) . Now, for any u ∈ R di-1 , we write it as u = q l=1 a l u l where a l is a scalar and u l is a unit vector supported on S l . By the definition of 1 2 -net, for each l ∈ [q] there exists v l ∈ N l such that v l -u l ≤ 1 2 . We know that A j:i (0)v l ≤ e -Ω(L γ ) v l for all l ∈ [q]. Let v = q l=1 a l v l . We have

A j:i (0)v ≤ q l=1 |a l | • A j:i (0)v l ≤ q l=1 |a l | • e -Ω(L γ ) v l ≤ e -Ω(L γ ) q • q l=1 a 2 l v l 2 = √ qe -Ω(L γ ) v = e -Ω(L γ ) v . Note that u -v = q l=1 a l (u l -v l ) = q l=1 a 2 l u l -v l 2 ≤ 1 4 q l=1 a 2 l = 1 2 u , which implies v ≤ 3 2 u . Therefore we have A j:i (0)u ≤ A j:i (0)v + A j:i (0)(u -v) ≤ e -Ω(L γ ) v + A j:i (0) • u -v ≤ e -Ω(L γ ) • 3 2 u + A j:i (0) • 1 2 u = e -Ω(L γ ) u + A j:i (0) • 1 2 u .

The above inequality is valid for any u ∈ R di-1 . Thus we can take the unit vector u that maximizes A j:i (0)u . This gives us A j:i (0) ≤ e -Ω(L γ ) + 1 2 A j:i (0) , which implies A j:i (0) ≤ e -Ω(L γ ) . Finally, we take a union bound over all possible (i, j). The failure probaility is at most L 2 e -Ω(L γ ) = e -Ω(L γ ) .

The following lemma shows that the properties in Lemmas 5.2 and 5.3 are still to some extent preserved after applying small perturbations on all the weight matrices. Lemma 5.4. Suppose that the initial weights satisfy A j:i (0) ≤ O(L 3 ) for all 1 ≤ i ≤ j ≤ L, and A j:i (0) ≤ e -c1L γ if j -i ≥ L 10 , where c 1 > 0 is a universal constant. Then for another set of matrices A 1 , . . . , A L satisfying A i -A i (0) ≤ e -0.6c1L γ for all i ∈ [L], we must have

EQUATION

)

Proof. It suffices to show that the difference

A j:i -A j:i (0) is tiny. Let ∆ i = A i -A i (0). We have A j:i = (A j (0) + ∆ j ) • • • (A i+1 (0) + ∆ i+1 )(A i (0) + ∆ i ).

Expanding this product, except for the one term corresponding to A j:i (0), every other term has the form

A j:(ks+1) (0) • ∆ ks • A (ks-1):(ks-1+1) (0) • ∆ ks-1 • • • ∆ k1 • A (k1-1):i (0), where i ≤ k 1 < • • • < k s ≤ j.

By assumption, each ∆ k has spectral norm e -0.6c1L γ , and each A j :i (0) has spectral norm O(L 3 ), so we have

A j:(ks+1) (0) • ∆ ks • A (ks-1):(ks-1+1) (0) • ∆ ks-1 • • • ∆ k1 • A (k1-1):i (0) ≤ e -0.6c1L γ s O(L 3 )

s+1 . Therefore we have

A j:i -A j:i (0) ≤ j-i+1 s=1 j -i + 1 s e -0.6c1L γ s O(L 3 ) s+1 ≤ j-i+1 s=1 L s e -0.6c1L γ s O(L 3 ) s+1 ≤ O(L 3 ) ∞ s=1 O(L 4 )e -0.6c1L γ s ≤ O(L 3 ) ∞ s=1 (1/2) s = O(L 3 ), which implies A j:i ≤ O(L 3 ) for all 1 ≤ i ≤ j ≤ L.

The proof of the second part of the lemma is postponed to Appendix B.2.

As a consequence of Lemma 5.4, we can control the objective value and the gradient at any point sufficiently close to the random initialization. Lemma 5.5. For a set of weight matrices W 1 , . . . , W L with A i = W i /( √ d i σ i ) that satisfy (12), the objective and the gradient satisfy

0.4 Y 2 F < (W 1 , . . . , W L ) < 0.6 Y 2 F , ∇ Wi (W 1 , . . . , W L ) ≤ ( d i σ i ) -1 e -0.9c1L γ , ∀i ∈ [L].

The proof of Lemma 5.5 is given in Appendix B.3.

Finally, we can finish the proof of Theorem 5.1 using the above lemmas.

Proof of Theorem 5.1. From Lemmas 5.2 and 5.3, we know that with probability at least 0.9, we have (i) A j:i (0) ≤ O(L 3 ) for all 1 ≤ i ≤ j ≤ L, and (ii)

A j:i (0) ≤ e -c1L γ if (i, j) further satisfies j -i ≥ L 10 .

Here c 1 > 0 is a universal constant. From now on we are conditioned on these properties being satisfied. We suppose that the learning rate η is at most e 0.2c1L γ .

We say that a set of weight matrices W 1 , . . . , W L are in the "initial neighborhood" if

A i -A i (0) ≤ e -0.6c1L γ for all i ∈ [L].

From Lemmas 5.4 and 5.5 we know that in the "initial neighborhood" the objective value is always between 0.4 Y 2 F and 0.6 Y 2 F . Therefore we have to escape the "initial neighborhood" in order to get the objective value out of this interval. Now we calculate how many iterations are necessary to escape the "initial neighborhood." According to Lemma 5.5, inside the "initial neighborhood" each W i can move at most η( √ d i σ i ) -1 e -0.9c1L γ in one iteration by definition of the gradient descent algorithm. In order to leave the "initial neighborhood," some

W i must satisfy W i -W i (0) = √ d i σ i A i -A i (0) > √ d i σ i e -0.6c1L γ .

In order to move this amount, the number of iterations has to be at least

√ d i σ i e -0.6c1L γ η( √ d i σ i ) -1 e -0.9c1L γ = d i σ 2 i e 0.3c1L γ η ≥ 1 L O(1) • e 0.

3c1L γ e 0.2c1L γ ≥ e Ω(L γ ) . This finishes the proof. 

6. EXPERIMENTS:

In this section, we provide empirical evidence to support the results in Sections 4 and 5. To study how depth and width affect convergence speed of gradient descent under orthogonal and Gaussian initialization schemes, we train a family of linear networks with their widths ranging from 10 to 1000 and depths from 1 to 700, on a fixed synthetic dataset (X, Y ).4 Each network is trained using gradient descent staring from both Gaussian and orthogonal initializations. In Figure 1 , We lay out the logarithm of the relative training loss (t) (0) , using heap-maps, at steps t = 1258 and t = 10000. In each heat-map, each point represents the relative training loss of one experiment; the darker the color, the smaller the loss. Figure 1 clearly demonstrates a sharp transition from untrainable to trainable (i.e., from red to black) when we increase the width of the network:

• for Gaussian initialization, this transition occurs across a contour characterized by a linear relation between width and depth;

• for orthogonal initialization, the transition occurs at a width that is approximately independent of the depth.

These observations excellently verify our theory developed in Sections 4 and 5.

To have a closer look into the training dynamics, we also plot "relative loss v.s. training time" for a variety of depth-width configurations. See Figure 2 . There again we can clearly see that orthogonal initialization enables fast training at small width (independent of depth), and that the required width for Gaussian initialization depends on depth. 

7. CONCLUSION:

In this work, we studied the effect of the initialization parameter values of deep linear neural networks on the convergence time of gradient descent. We found that when the initial weights are iid Gaussian, the convergence time grows exponentially in the depth unless the width is at least as large as the depth. In contrast, when the initial weight matrices are drawn from the orthogonal group, the width needed to guarantee efficient convergence is in fact independent of the depth. These results establish for the first time a concrete proof that orthogonal initialization is superior to Gaussian initialization in terms of convergence time.

A PROOFS FOR SECTION 4

A.1 PROOF OF LEMMA 4.2

Proof of Lemma 4.2. We only need to prove (9). We first upper bound the magnitude of the network's initial output on any given input  and f (x; W 1 (0 

x ∈ R dx . Let z = 1 √ m L-1 W L-1:1 (0) • x ∈ R m . Then we have z = x ,

), . . . , W L (0)) = 1 √ dy W L (0) • z = m dy • 1 √ m W L (0) • z. Note that 1 √ m W L (0) • z is the (signed) projection of z onto a random subspace in R m of dimension d y . Therefore 1 √ m W L (0) • z 2 / z

2 has the same distribution as

g 2 1 +•••+g 2 dy g 2 1 +•••+g 2 m

, where g 1 , . . . , g m are i.i.d. samples from N (0, 1). By the standard tail bounds for chi-squared distributions we have

Pr g 2 1 + • • • + g 2 dy ≤ d y + 2 d y log(1/δ ) + 2 log(1/δ ) ≥ 1 -δ , Pr g 2 1 + • • • + g 2 m ≥ m -2 m log(1/δ ) ≥ 1 -δ .

Let δ = δ 2r . Note that m > C • log(r/δ). We know that with probability at least 1 -δ r we have

EQUATION

Finally, taking a union bound, we know that with probability at least 1 -δ, the inequality (13) holds for every x ∈ {x 1 , . . . , x r }, which implies

(0) = 1 2 r k=1 f (x k ; W 1 (0), . . . , W L (0)) -y k 2 ≤ r k=1 f (x k ; W 1 (0), . . . , W L (0)) 2 + y k 2 ≤ O 1 + log(r/δ) d y r k=1 x k 2 + r k=1 y k 2 = O 1 + log(r/δ) d y X 2 F + Y 2 F ≤ O 1 + log(r/δ) d y + W * 2 X 2 F . A.2 PROOF OF CLAIM 4.3 Proof of Claim 4.3. Let γ = 1 2 Lσ 2 min (X)/d y . From A(0), . . . , A(t) we have (s) ≤ (1-ηγ) s B for all 0 ≤ s ≤ t. The gradient of the objective function (2) is ∂ ∂Wi = αW L:i+1 (U -Y ) (W i-1:1 X) .

Thus we can bound the gradient norm as follows for all 0 ≤ s ≤ t and all i ∈ [L]:

EQUATION

where we have used B(s). Then for all i ∈ [L] we have:

W i (t + 1) -W i (0) F ≤ t s=0 W i (s + 1) -W i (s) F = t s=0 η ∂ ∂W i (s) F ≤ η t s=0 2 (1 -ηγ) s B d y X ≤ 2η √ B d y X t-1 s=0 (1 -ηγ/2) s ≤ 2η √ B d y X • 2 ηγ = 8 Bd y X Lσ 2 min (X)

.

This proves C(t + 1).

A.3 PROOF OF CLAIM 4.4 Proof of Claim 4.4. Let R = 8 √ Bdy X Lσ 2 min (X) and ∆ i = W i (t) -W i (0) (i ∈ [L]). Then C(t) means ∆ i F ≤ R (∀i ∈ [L]).

For 1 ≤ i ≤ j ≤ L, we have

W j:i (t) = (W j (0) + ∆ j ) • • • (W i (0) + ∆ i ) .

Expanding this product, each term except W j:i (0) has the form:

EQUATION

where i ≤ k 1 < • • • < k s ≤ j are locations where terms like ∆ k l are taken out. Note that every factor in (15) of the form W j :i (0) satisfies W j :i (0) = m j -i +1 2

according to (8). Thus, we can bound the sum of all terms of the form (15) as

W j:i (t) -W j:i (0) ≤ j-i+1 s=1 j -i + 1 s R s m j-i+1-s 2 = ( √ m + R) j-i+1 -( √ m) j-i+1 = ( √ m) j-i+1 1 + R/ √ m j-i+1 -1 ≤ ( √ m) j-i+1 1 + R/ √ m L -1 ≤ 0.1( √ m) j-i+1 .

Here the last step uses m > C(LR) 2 which is implied by (5). Combined with (8), this proves B(t).

A.4 PROOF OF CLAIM 4.5

Proof of Claim 4.5. Recall that we have the dynamics (6) for U (t). In order to establish convergence from (6) we need to prove upper and lower bounds on the eigenvalues of P (t), as well as show that the high-order term E(t) is small. We will prove these using B(t).

Using the definition (7) and property B(t), we have

λ max (P (t)) ≤ α 2 L i=1 λ max (W i-1:1 (t)X) (W i-1:1 (t)X) • λ max W L:i+1 (t)W L:i+1 (t) ≤ 1 m L-1 d y L i=1 1.1m i-1 2 σ max (X) 2 1.1m L-i 2 2 ≤ 2Lσ 2 max (X)/d y , λ min (P (t)) ≥ α 2 L i=1 λ min (W i-1:1 (t)X) (W i-1:1 (t)X) • λ min W L:i+1 (t)W L:i+1 (t) ≥ 1 m L-1 d y L i=1 0.9m i-1 2 σ min (X) 2 0.9m L-i 2 2 ≥ 3 5 Lσ 2 min (X)/d y .

In the lower bound above, we make use of the following relation on dimensions: m ≥ d x ≥ r, which enables the inequality λ min (W i-1:

1 (t)X) (W i-1:1 (t)X) = σ 2 min (W i-1:1 (t)X) ≥ σ 2 min (W i-1:1 (t)) • σ 2 min (X)

. Next, we will prove the following bound on the high-order term E(t):

1 m L-1 d y E(t)X F ≤ 1 6 ηλ min (P t ) U (t) -Y F .

Recall that E(t) is the sum of all high-order terms in the product

W L:1 (t + 1) = i W i (t) -η ∂ ∂W i (t) .

Same as ( 14), we have ∂ ∂Wi (t)

F ≤ 2 √ (t) X √ dy (∀i ∈ [L]

). Then we have Proof of Lemma 5.2. Notice that for any 1 ≤ i ≤ j ≤ L we have E A j:i (0) 2 F = d i-1 . Then by Markov inequality we have Pr A j:i (0) 2 F ≥ di-1 δ/L 2 ≤ δ/L 2 . Taking a union bound, we know that with probability at least 1 -δ, for all 1 ≤ i ≤ j ≤ L simultaneously we have A j:i (0) ≤ A j:i (0) F ≤ di-1 δ/L 2 ≤ O(L 3 /δ) (note that d i-1 ≤ O(L 1-γ ) = O(L)).

EQUATION

B.2 PROOF OF LEMMA 5.4

Proof of Lemma 5.4 (continued). For the second part of the lemma (j -i ≥ L 4 ), we need to bound the terms of the form A j:k+1 (0) • ∆ k • A k-1:i (0) more carefully. In fact, if j -i ≥ L 4 , then max{j -k -1, k -1 -i} ≥ L 10 , which by assumption means either A j:k+1 (0) or A k-1:i (0) has spectral norm bounded by e -c1L γ . This implies A j:k+1 (0) • ∆ k • A k-1:i (0) ≤ e -c1L γ e -0.6c1L γ • O(L 3 ) = e -1.6c1L γ • O(L 3 ). Therefore we have A j:i -A j:i (0) ≤ (j -i + 1)e -1.6c1L γ

• O(L 3 ) + j-i+1 s=2

j -i + 1 s e -0.6c1L γ s O(L 3 ) s+1 ≤ e -c1L γ + ∞ s=2

L s e -0.6c1L γ s O(L 3 ) s+1 ≤ e -c1L γ + ∞ s=2 e -0.5c1L γ s = O e -c1L γ .

This implies A j:i ≤ O e -c1L γ .

B.3 PROOF OF LEMMA 5.5

Proof of Lemma 5.5. We can bound the network's output as αW L:1 (0)X F = βA L:1 (0)X F ≤ L O(1) • e -Ω(L γ ) X F = e -Ω(L γ ) .

Thus the objective value (W 1 , . . . , W L ) = 1 2 αW L:1 (0)X -Y 2 F must be extremely close to 

X -Y ) X W i-1:1 = β/( d i σ i ) • A L:i+1 (αW L:1 X -Y ) X A i-1:1 ≤ L O(1) √ d i σ i A L:i+1 • O(1) • A i-1:1 .

Using (12), and noting that either L -i -1 or i -1 is greater than L 4 , we have

∇ Wi (W 1 , . . . , W L ) ≤ σ -1 i L O(1) • O e -c1L γ • O(L 3 ) ≤ ( d i σ i )

-1 e -0.9c1L γ .

:

We have E f (x; WL(0), . . . , W1(0))

= α 2 E x W 1 (0) • • • W L (0)WL(0) • • • W1(0)x . Note that by our choice (4) we have E W L (0)WL(0) = dyIm and W i (0)Wi(0) = mI (1 ≤ i ≤ L -1), so we have E f (x; WL(0), . . . , W1(0)) 2 = α 2 m L-1 dy x 2 = x 2 .2 r is known as the stable rank of X, which is always no more than the rank.

Note that for symmetric matrices A and B, the set of eigenvalues of A ⊗ B is the set of products of an eigenvalue of A and an eigenvalue of B.

We choose X ∈ R 1024×16 and W * ∈ R 10×1024 , and set Y = W * X. Entries in X and W * are drawn i.i.d. from N (0, 1).

Figures and Tables:

Figure 1: log (t) (0) at t = 1258 and t = 10000, for different depth-width configurations and different initialization schemes. Darker color means smaller loss.

Figure 2: Relative loss v.s. training time. For each plot, we vary width from 50 (yellow) to 1200 (purple). Solid and dashed lines represent Gaussian (GS) and orthogonal (OT) initializations.

for large L, so 0.4 Y 2 < (W 1 , . . . , W L ) < 0.6 Y 2 F .As for the gradient, for any i ∈ [L] we have ∇ Wi (W 1 , . . . , W L ) = αW L:i+1 (αW L:1

It suffices to show that the above bound is at most 1 6 ηλ min (P t ) U (t) -Y F = 1 6 ηλ min (P t ) 2 (t). Combined with A(t), this proves A(t + 1).

Table content: <table><tr><td>t) X md y √ dy• (t) √ m X . Note that m &gt; C • s X 2ηL (t) X md y 2 L-2 s=2 ≤ 2 m d y X • 2ηL (t) X dyB X 2 ≥ C • 0.5 s-2 md y • d y • (t) dy (t) X 2 . Thus √ m X min (X)/d y , it suffices to have 2L X 2 , we have dy 2ηL √ (t) X √ mdy ≤ 1 m L-1 d y E(t)X F ≤ m d y ≤ 2 From η ≤ we have m d y X 2ηL (t) X md y 2 = 4ηL X • (t) md y . Since λ min (P t ) ≥ 3 5 Lσ 2 4ηL X • (t) md y ≤ 1 6 η • 3Lσ 2 min (X) 2 (t) 5d y , which is true since m &gt; C • dyB X 2 σ 4 min (X) ≥ C • dy (t) X 2 σ 4 min (X) . Finally, from (6) and η ≤ dy 2L X 2 ≤ 1 λmax(Pt) we have U (t + 1) -Y F = vec (U (t + 1) -Y ) = (I -ηP (t)) • vec (U (t) -Y ) + 1 m L-1 d y vec (E(t)X) ≤ (1 -ηλ min (P (t))) vec (U (t) -Y ) + 1 m L-1 d y E(t)X F ≤ (1 -ηλ min (P (t))) U (t) -Y F + 1 6 ηλ min (P t ) U (t) -Y F = 1 -5 6 ηλ min (P (t)) U (t) -Y F ≤ 1 -1 2 ηLσ 2 min (X)/d B PROOFS FOR SECTION 5 B.1 PROOF OF LEMMA 5.2</td></tr></table>

Peer review (iclr2020):
This paper studies the convergence of deep linear networks under orthogonal initialization. Most of the problem setup, analysis techniques and proof roadmap are adapted from Du & Hu (2019). The difference is the orthogonal initialization compared with the Gaussian initialization in Du & Hu (2019). Since the product of orthogonal weight matrices is identity, the new initialization can remove the dependency of the number of nodes on the depth of the network. The authors also proved a lower bound of the loss function trained by Gaussian initialization within certain iterations. This justifies the disadvantage of Gaussian initialization. My biggest concern is that this paper seems to be very similar to Du & Hu (2019) in many places. The whole Sections 3 & 4 are almost the same as Sections 3,4,5 & 7 in their paper. The contribution in this paper is too incremental given previous work. What is the advantage of using orthogonal matrices in (4) compared with just using identity matrices as initialization? Can we prove the same result with just identity initialization? In this case, what would we lose by restricting us to this special case? In the proof of Lemma 4.2, it seems that only the randomness of the last layer is used. Why do we need all the layers to be uniformly sampled? It should be explained in more details that is drawn from a uniform distribution over orthogonal matrices in space. Then is not distributed on the whole space of -sphere. The argument in the proof of Lemma 4.2 thus needs more justification. It would be interesting to see an empirical comparison of the proposed initialization and the Gaussian initialization. Due to the lower bound proved in this paper, the experiments are expected to show distinct difference between these two.