GRAM-GAUSS-NEWTON METHOD: LEARNING OVER-
Year: 2019

Abstract:
First-order methods such as stochastic gradient descent (SGD) are currently the standard algorithm for training deep neural networks. Second-order methods, despite their better convergence rate, are rarely used in practice due to the prohibitive computational cost in calculating the second-order information. In this paper, we propose a novel Gram-Gauss-Newton (GGN) algorithm to train deep neural networks for regression problems with square loss. Our method draws inspiration from the connection between neural network optimization and kernel regression of neural tangent kernel (NTK). Different from typical second-order methods that have heavy computational cost in each iteration, GGN only has minor overhead compared to first-order methods such as SGD. We also give theoretical results to show that for sufficiently wide neural networks, the convergence rate of GGN is quadratic. Furthermore, we provide convergence guarantee for mini-batch GGN algorithm, which is, to our knowledge, the first convergence result for the mini-batch version of a second-order method on overparameterized neural networks. Preliminary experiments on regression tasks demonstrate that for training standard networks, our GGN algorithm converges much faster and achieves better performance than SGD.

heading: 1 INTRODUCTION
text: First-order methods such as Stochastic Gradient Descent (SGD) are currently the standard choice for training deep neural networks. The merit of first-order methods is obvious: they only calculate the gradient and therefore are computationally efficient. In addition to better computational efficiency, SGD has even more advantages among the first-order methods. At each iteration, SGD computes the gradient only on a mini-batch instead of all training data. Such randomness introduced by sampling the mini-batch can lead to better generalization (Hardt et al., 2015; Keskar et al., 2016; Masters & Luschi, 2018; Mou et al., 2017; Zhu et al., 2018) and better convergence (Ge et al., 2015; Jin et al., 2017a;b), which is crucial when the function class is highly overparameterized deep neural networks. Recently there is a huge body of works trying to develop more efficient first-order methods beyond SGD (Duchi et al., 2011; Kingma & Ba, 2014; Luo et al., 2019; Liu et al., 2019).
Second-order methods, despite their better convergence rate, are rarely used to train deep neural networks. At each iteration, the algorithm has to compute second order information, for example, the Hessian or its approximation, which is typically an m by m matrix where m is the number of parameters of the neural network. Moreover, the algorithm needs to compute the inverse of this matrix. The computational cost is prohibitive and usually it is not even possible to store such a matrix.
Recently, there is a series of works (Du et al., 2018b;a; Zou et al., 2018; Allen-Zhu et al., 2018a; Oymak & Soltanolkotabi, 2018; Arora et al., 2019b;a; Cao & Gu, 2019; Zou & Gu, 2019) considering the optimization of neural networks based on the idea of neural tangent kernel (NTK) (Jacot et al., 2018). Roughly speaking, the idea of NTK is to linearly approximate the output of a network w.r.t. the parameters at a local region, thus resulting in a kernel feature map, which is the gradient of the output w.r.t. the parameters. Jacot et al. (2018) shows that when the width of the network tends to infinity, the NTK tends to be unchanged during gradient flow since the network only needs o(1) change (goes to zero as the width tends to infinity) of parameter to fit the data, so the change of kernel
is also o(1). As for finite-width networks, however, the NTK changes and previous works (Allen-Zhu et al., 2018b; Arora et al., 2019a) show that for sufficiently wide neural networks the optimization dynamic of GD/SGD is equivalent to that of using GD/SGD to solve an NTK kernel regression problem where the kernel is slowly evolving (see Lemma 1 in Section 2 for a precise description). A natural question then arises:
Can we gain acceleration by directly solving kernel regression w.r.t. the NTK at each step?
In this paper, we give a positive answer to this question and reveal the connection between NTK regression and Gauss-Newton method. We propose a novel optimization method – the Gram-GaussNewton (GGN) method. Instead of doing gradient descent, GGN solves the kernel regression w.r.t. the NTK at each step of the optimization. Following this idea, we theoretically prove that for overparameterized networks, GGN enjoys quadratic convergence compared to the linear rate of gradient descent.
Besides theoretical fast convergence rate, GGN is also very efficient in practice. In fact, GGN is implicitly a reformulation of the Gauss-Newton method (see Section 3.1 for details) which is a classic second-order algorithm often used for solving nonlinear regression problems with square loss. In the Gauss-Newton method, one uses J>J as an approximation of the Hessian (see Section 2 for a formal description) where J is the Jacobian matrix. However, the original Gauss-Newton method faces challenges when used for training deep neural networks. Most seriously, the size of the approximate Hessian J>J is m by m, where m is the number of parameters of the neural network. Moreover, for overparameterized neural networks, J>J is not invertible, which may make the algorithm intractable for training commonly-used neural networks.
GGN bypasses the difficulty stated above as follows. Instead of using J>J as approximate Hessian and applying Newton-type method, each step of GGN only involves the Gram matrix JJ> whose size is n by n where n is the number of data. Furthermore, as already mentioned, to get better generalization performance, it is crucial to use mini-batch to introduce sampling noise when calculating derivatives. Therefore, like SGD, we also use mini-batch in GGN. In this case, the size of the Gram matrix further reduces to b by b, where b is the batch size. Though conventional wisdom may suggest that applying mini-batch scheme to second-order methods will introduce biased estimation of the accelerated gradient direction, we give the first convergence result for mini-batch second-order method on overparameterized networks. Regarding computational complexity, we show that at each iteration, the overhead of GGN is small compared to SGD: the extra computation of GGN is mainly the matrix product JJ> and the inverse of this matrix whose size is small for a mini-batch. Detailed analyses can be found in Section 3.3. We next conduct experiments on two regression tasks to study the effectiveness of the GGN algorithm. We demonstrate that in these two real applications, using a practical neural network (e.g., ResNet-32) with standard width, our proposed GGN algorithm can converge faster and achieve better performance than several baseline algorithms.

heading: 1.1 RELATED WORKS
text: Despite the prevalence of first-order methods for training deep neural networks, there have been continuing efforts in developing practical second-order methods (Becker et al., 1988; Pascanu & Bengio, 2013). We summarize some of these works below.
The main approach for these methods is to develop delicate approximations of the second-order information matrix so that the update direction can be computed as efficiently as possible. For example, Botev et al. (2017) proposed a recursive block-diagonal approximation of the Hessian. The blocks are Kronecker factored and can be efficiently computed and inverted. Grosse and Martens in a series of works developed the K-FAC method (Martens & Grosse, 2015; Grosse & Martens, 2016). The key idea is a Kronecker-factored approximation of the Fisher information matrix, which is used as the second-order matrix in natural gradient methods. These works received considerable attention and have been further improved (Wu et al., 2017; George et al., 2018; Martens et al., 2018). Bernacchia et al. (2018) derived an exact expression of the natural gradient update, but only works for linear networks. Different from all these works, our GGN algorithm does not try to approximate the second-order matrix whose size is inevitably huge. Instead, we present an easy-to-compute solution of the updating direction, reducing the computational cost significantly. One exceptional concurrent work Ren & Goldfarb (2019) also aims to use the exact Gauss-Newton update. They focus on reducing the complexity of inverting approximate Hessian by Sherman-Morrison-Woodbury
Formula and require subtle implementation tricks to use backpropagation. In contrast, GGN has simpler update rule and better guarantee for neural networks.
In a concurrent and independent work, Zhang et al. (2019a) showed that natural gradient method and K-FAC have a linear convergence rate for sufficiently wide networks in full-batch setting. In contrast, our method enjoys a higher-order (quadratic) convergence rate guarantee for overparameterized networks, and we focus on developing a practical and theoretically sound optimization method. We also reveal the relation between our method and NTK kernel regression, so using results based on NTK (Arora et al., 2019b), one can easily give generalization guarantee of our method. Another independent work (Achiam et al., 2019) proposed a preconditioned Q-learning algorithm which has similar form of our update rule. Unlike the methods considered in Zhang et al. (2019a); Achiam et al. (2019) which contain the learning rate that needed to be tuned, our derivation of GGN does not introduce a learning rate term (or understood as suggesting that the learning rate can be fixed to be 1 to get good performance which is verified in Figure 2 (c)).

heading: 2 NEURAL TANGENT KERNEL AND THE CLASSIC GAUSS-NEWTON METHOD FOR NONLINEAR LEAST SQUARES REGRESSION
text: Nonlinear least squares regression problem is a general machine learning problem. Given data pairs {xi, yi}ni=1 and a class of nonlinear functions f , e.g. neural networks, parameterized by w, the nonlinear least squares regression aims to solve the optimization problem
min w∈Rm
L(w) = 1
2 n∑ i=1 (f(w,xi)− yi)2. (1)
In the seminal work (Jacot et al., 2018), the authors consider the case when f is a neural network with infinite width. They showed that optimization on this problem using gradient flow involves a special kernel which is called neural tangent kernel (NTK). The follow-up works further extended the relation between optimization and NTK which can be concluded in the following lemma: Lemma 1 (Lemma 3.1 in Arora et al. (2019a), see also Dou & Liang (2019); Mei et al. (2019)). Consider optimizing problem (1) by gradient descent with infinitesimally small learning rate: dwtdt = −∇L(wt). where wt is the parameters at time t. Let ft = (f(wt,xi))ni=1 ∈ Rn be the network outputs on all xi’s at time t, and y = (yi)ni=1 be the desired outputs. Then ft follows the following evolution:
dft dt
= −Gt · (ft − y), (2) where Gt is an n × n positive semidefinite matrix, i.e. the Gram matrix w.r.t. the NTK at time t, whose (i, j)-th entry is 〈∇wf(wt,xi),∇wf(wt,xj)〉.
The key idea of Jacot et al. (2018) and its extensions (Du et al., 2018b;a; Zou et al., 2018; Allen-Zhu et al., 2018a; Oymak & Soltanolkotabi, 2018; Lee et al., 2019; Yang, 2019; Arora et al., 2019b;a; Cao & Gu, 2019; Zou & Gu, 2019) is that when the network is sufficiently wide, the Gram matrix at initialization G0 is close to a fixed positive definite matrix defined by the infinite-width kernel and Gt is close to G0 during training for all t. Under this situation, Gt remains invertible, and the above dynamics is then identical to the dynamics of solving kernel regression with gradient flow w.r.t. the current kernel at time t. In fact, Arora et al. (2019a) rigorously proves that a fully-trained sufficiently wide ReLU neural network is equivalent to the kernel regression predictor.
As pointed out in Chizat & Bach (2018), the idea of NTK can be summarized as a linear approximation using first order Taylor expansion. We give an example of this idea on the NTK at initialization:
f(w,xi)− f(w0,xi) ≈ ∇wf(w0,xi) · (w −w0), (3) where∇wf(w0,x) can then be viewed as an explicit expression of feature map at x, w −w0 is the parameter in reproducing kernel Hilbert space (RKHS) induced by NTK and f(w,xi)− f(w0,xi) the target value.
The idea of linear approximation is also used in the classic Gauss-Newton method (Golub, 1965) to obtain an acceleration algorithm for solving nonlinear least squares problem (1). Concretely, at iteration t, Gauss-Newton method takes the following first-order approximation:
f(w,xi)− f(wt,xi) ≈ ∇wf(wt,xi) · (w −wt), (4)
where wt stands for the parameter at iteration t. We note that this is also the linear expansion for deriving NTK at time t. According to Eq. (1) and (4), to update the parameter, one can instead solve the following problem.
wt+1 = argmin w
1 2 ‖ft + Jt(w −wt)− y‖22 , (5)
where ft,y have the same meaning as in Lemma 1, and Jt = (∇wf(wt,x1), · · · ,∇wf(wt,xn))> ∈ Rn×m is the Jacobian matrix.
A necessary and sufficient condition for w to be the solution of Eq. (5) is
(J>t Jt) · (w −wt) = −J>t (ft − y). (6)
Below we will denote Ht := J>t Jt ∈ Rm×m. For under-parameterized model (i.e., the number of parameters m is less than the number of data n), Ht is invertible, and the update rule is
wt+1 = wt −H−1t J>t (ft − y). (7)
This can also be viewed as an approximate Newton’s method using Ht = J>t Jt to approximate the Hessian matrix. In fact, the exact Hessian matrix is
∇2w 1
2 n∑ i=1 (f(wt,xi)− yi)2 = J>t Jt + n∑ i=1 (f(wt,xi)− yi)∇2wf(wt,xi). (8)
In the case when f is only mildly nonlinear w.r.t. w at data point xi’s,∇2wf(wt,xi) ≈ 0, and Ht is close to the real Hessian. In this situation, the behavior of the Gauss-Newton method is similar to that of Newton’s method, and thus can achieve a superlinear convergence rate (Golub, 1965).

heading: 3 THE GRAM-GAUSS-NEWTON METHOD
text: The classic second-order methods using approximate Hessian such as Gauss-Newton method described in the previous section face obvious difficulties dealing with the intractable approximate Hessian matrix when the regression model is an overparameterized neural network.
In Section 3.1, we develop a Gram-Gauss-Newton (GGN) method which is inspired by NTK kernel regression and does not require the computation of the approximate Hessian. In Section 3.2, we show that for sufficiently wide neural networks, GGN has quadratic convergence rate. In Section 3.3, we show that the additional computational cost (per iteration) of GGN compared to SGD is small.

heading: 3.1 THE GRAM-GAUSS-NEWTON METHOD FOR OVERPARAMETERIZED NEURAL NETWORKS
text: We now describe our GGN method to learn overparameterized neural networks for regression problems. As mentioned in the previous sections, for sufficiently wide networks, using gradient descent for solving the regression problem (1) has similar dynamics as using gradient descent for solving NTK kernel regression (Lemma 1) w.r.t. NTK at each step. However, one can also solve the kernel regression problem w.r.t. the NTK at each step immediately using the explicit formula of kernel regression. By explicitly solving instead of using gradient descent to solve NTK kernel regression, one can expect the optimization to get accelerated. We propose our Gram-Gauss-Newton (GGN) method to directly solve the NTK kernel regression with Gram matrix Gt at each time step t. Note that the feature map of NTK at time t, based on the approximation in Eq. (4), can be expressed as x 7→ ∇wf(wt,x), and the linear parameters in RKHS are w −wt, also the target is f(w,xi)− f(wt,xi). Therefore, the kernel (ridgeless) regression solution (Mohri et al., 2018; Liang & Rakhlin, 2018) of Jt,S(w −wt) = (f(w,xi)− f(wt,xi)) w.r.t. (w −wt) gives the update
wt+1 = wt − J>t,SG−1t,S(ft,S − yS), (9)
where Jt,S is the matrix of features at iteration t computed on the training data set S which is equal to the Jacobian , ft,S and yS are the vectorized outputs of neural network and the corresponding targets on S respectively, and
Gt,S = Jt,SJ > t,S
is the Gram matrix of the NTK on S.
One may wonder what is the relation between our derivation from NTK kernel regression and the Gauss-Newton method. We point out that for overparameterized models, there are infinitely many solutions of Eq. (5) but our update rule (9) essentially uses the minimum norm solution. In other words, the GGN update rule re-derives the Gauss-Newton method with the minimum norm solution. This somewhat surprising connection is due to the fact that in kernel learning, people usually choose a kernel with powerful expressivity, i.e. the dimension of feature space is large or even infinite. However, by the representer theorem (Mohri et al., 2018), the solution of kernel (ridgeless) regression lies in the n-dimensional subspace of RKHS and minimizes the RKHS norm. We refer the readers to Chapter 11 of Mohri et al. (2018) for details.
As mentioned in Section 1, the design of learning algorithms should consider not only optimization but also generalization. It has been shown that using mini-batch instead of full batch to compute derivatives is crucial for the learned model to have good generalization ability (Hardt et al., 2015; Keskar et al., 2016; Masters & Luschi, 2018; Mou et al., 2017; Zhu et al., 2018). Therefore, we propose a mini-batch version of GGN. The update rule is the following:
wt+1 = wt − J>t,BtG −1 t,Bt (ft,Bt − yBt), (10) where Bt is the mini-batch used at iteration t, Jt,Bt and Gt,Bt are the Jacobian and the Gram matrix computed using the data of Bt respectively, and ft,Bt ,yBt are the vectorized outputs and the corresponding targets on Bt respectively. Gt,Bt = Jt,BtJ > t,Bt
is a very small matrix when using a typical batch size.
One difference between Eq. (10) and Eq. (7) is that our update rule only requires to compute the Gram matrix Gt,Bt and its inverse. Note that the size of Gt,Bt is equal to the size of the mini-batch and is typically very small. So this also greatly reduces the computational cost.
Using the idea of kernel ridge regression (which can also be viewed as Levenberg-Marquardt extension (Levenberg, 1944) of Gauss-Newton method), we introduce the following variant of GGN:
wt+1 = wt − J>t,Bt(λGt,Bt + αI) −1(ft,Bt − yBt), (11)
where λ > 0 is another hyper-parameter controlling the learning process. Our algorithm is formally described in Algorithm 1.
Algorithm 1 (Mini-batch) Gram-Gauss-Newton Method 1: Input: Training dataset S. Hyper-parameters λ and α. 2: Initialize the network parameter w0. Set t = 0. 3: for each iteration do 4: Fetch a mini-batch Bt from the dataset. 5: Calculate the Jacobian matrix Jt,Bt . 6: Calculate the Gram matrix Gt,Bt = Jt,BtJ > t,Bt
. 7: Update the parameter by wt+1 = wt − J>t,Bt(λGt,Bt + αI)
−1(ft,Bt − yBt). 8: t = t+ 1. 9: end for

heading: 3.2 CONVERGENCE ANALYSIS FOR OVERPARAMETERIZED NEURAL NETWORKS
text: In this subsection, we show that for two-layer neural networks, if the width is sufficiently large, then: (1) Full-batch GGN converges with quadratic convergence rate. (2) Mini-batch GGN converges linearly. (For clarity, here we only present a proof for two-layer neural networks, but we believe that it is not hard for the conclusion to be extended to deep neural networks using the techniques developed in Du et al. (2018a); Zou & Gu (2019)).
As we explained through the lens of NTK, the result is a consequence of the fact that for wide enough neural networks, if the weights are initialized according to a suitable probability distribution, then with high probability the output of the network is close to a linear function w.r.t. the parameters (but nonlinear w.r.t. the input of the network) in a neighborhood containing the initialization point and a global optimum. Although the neural networks used in practice are far from that wide, this still motivates us to design the GGN algorithm.
Neural network structure. We use the following two-layer network
f(w,x) = 1√ M a>σ(W>x) = 1√ M M∑ r=1 arσ(wrx), (12)
where x ∈ Rd is the input, M is the network width, W = (w>1 , · · · ,w>M )>and σ(·) is the activation function. Each entry of W is i.i.d. initialized with the standard Gaussian distribution wr ∼ N (0, Id) and each entry of a is initialized from the uniform distribution on {±1}. Similar to Du et al. (2018b), we only train the network on parameter W just for the clarity of the proof. We also assume the activation function σ(·) is `-Lipschitz and β-smooth, and ` and β are regarded as O(1) absolute constants.
The key finding, as pointed out in Jacot et al. (2018); Du et al. (2018b;a), is that under such initialization, the Gram matrix G has an asymptotic limit, which is, under mild conditions (e.g. input data not being degenerate etc., see Lemma F.2 of Du et al. (2018a)), a positive definite matrix
K(xi,xj) = Ew∼N (0,I) [ x>i xjσ ′(wxi)σ ′(wxj) ] . (13)
Assumption 1 (Least Eigenvalue of the Limit of Gram Matrix). We assume the matrix K defined in (13) above is positive definite, and denote its least eigenvalue as
λ0 = λmin(K) > 0.
Now we are ready to state our theorem of full-batch GGN: Theorem 1 (Quadratic Convergence of Full-batch GGN on Overparameterized Neural Networks). Assume Assumption 1 holds. Assume the scale of the data is ‖xi‖2 = O(1), |yi| = O(1) for i ∈ {1, · · · , n}. If the network width
M = Ω ( max ( n4
λ40 , n2d log(16n/δ) λ20
)) ,
then with probability 1 − δ over the random initialization, the full-batch version of GGN whose update rule is given in Eq. (9) satisfies the following:
1) The Gram matrix Gt,S at each iteration is invertible;

heading: 2) The loss converges to zero in a way that
text: ‖ft+1 − y‖2 ≤ C√ M ‖ft − y‖22 (14)
for some C that is independent of M , which is a second-order convergence.
For the mini-batch version of GGN, by the analysis of its NTK limit, the algorithm is essentially doing serial subspace correction (Xu, 2001) on subspaces induced by mini-batch. So mini-batch GGN is similar to the Gauss-Siedel method (Golub & Van Loan, 1996) applied to solving systems of linear equations, as shown in the proof of the following theorem. Similar to the full batch situation, GGN takes the exact solution of the “kernel regression problem on the subspace” which is faster than just doing a gradient step to optimize on the subspace. Moreover, we note that existing results of the convergence of SGD on overparameterized networks usually use the idea that when the step size is bounded by a quantity related to smoothness, the SGD can be reduced to GD. However, our analysis takes a different way from the analysis of GD, thus does not rely on small step size.
In the following, we denote G0 ∈ Rn×n as the initial Gram matrix. Let n = bk, where b is the batch size and k is the number of batches, and let
G0,ij := G0((i− 1)b+ 1 : ib, (j − 1)b+ 1 : jb) be the (i, j)-th b× b block of G0. We define the iteration matrix A = L>(D− L)−1 ∈ Rn×n, (15) where
D =  G0,11 0 · · · 0 0 G0,22 · · · 0 ... ... . . . ...
0 0 . . . G0,kk
 , L = −  0 0 · · · 0 G0,21 0 · · · 0 ... ... . . . ...
G0,k1 G0,k2 . . . 0
 ,
represents the block-diagonal and block-lower-triangular parts of G0. We will show that the convergence of mini-batch GGN is highly related to the spectral radius of A. To simplify the proof, we make the following mild assumption on A: Assumption 2 (Assumption on the Iteration Matrix). Assume the matrix A defined in (15) above is diagonalizable. So we choose an arbitary diagonalization of A as A = P−1QP and denote
µ := ‖P‖2 ∥∥P−1∥∥ 2 .
We note that Assumption 2 is only for the sake of simplicity. Even if it does not hold, an infinitesimally small perturbation can make any matrix diagonalizable, and it will not affect the proof.
Now we are ready to state the theorem for mini-batch GGN. Theorem 2 (Convergence of Mini-batch GGN on Overparameterized Neural Networks). Assume Assumption 1 and 2 hold. Assume the scale of the data is ‖xi‖2 = O(1), |yi| = O(1) for i ∈ {1, · · · , n}. We use the mini-batch version of GGN whose update rule is given in Eq.(10), and the batch Bt is chosen sequentially and cyclically with a fixed batch size b and k = n/b updates per epoch. If the network width
M = max ( Ω ( µ2n18
λ160
) ,Ω ( n2d log(16n/δ)
λ20
)) ,
then with probability 1− δ over the random initialization, we have the following: 1) The Gram matrix Gt,Bt at each iteration is invertible;
2) The loss converges to zero in a way that after T epochs, we have
‖fTk − y‖2 ≤ µ √ n ( 1− Ω ( λ20 n2 ))T . (16)
Proof sketch for Theorem 1 and 2. Denote Jt = J(Wt) and
Jt,t+1 = ∫ 1 0 J((1− s)Wt + sWt+1)ds.
For the full-batch version, we have
‖ft+1 − y‖2 = ‖ft − y + Jt,t+1(vec(Wt+1)− vec(Wt))‖2 = ∥∥ft − y − Jt,t+1J>t G−1t (ft − y)∥∥2
= ∥∥(Jt − Jt,t+1)J>t G−1t (ft − y)∥∥2
≤ ‖Jt − Jt,t+1‖2 ∥∥J>t ∥∥2 ∥∥G−1t ∥∥2 ‖ft − y‖2 . (17)
Then we control the first term in Eq. (17) in a way similar to the following:
‖Jt − Jt,t+1‖2 ≤ Ĉ√ M ‖Wt −Wt+1‖2 ≤ Ĉ√ M ∥∥J>t ∥∥2 ∥∥G−1t ∥∥2 ‖ft − y‖2 , and if
∥∥J>t ∥∥2 ∥∥G−1t ∥∥2 can be upper bounded, we get our result Eq. (14). For the mini-batch version, similarly we have
ft+1 − y = ft − y − Jt,t+1J>t,BtG −1 t,Bt (ft,Bt − yBt)
= ft − y − Jt,t+1J>t,BtG̃ −1 t,Bt (ft − y), (18)
where the subscript Bt denotes the sub-matrix/vector corresponding to the batch, and G̃−1t,Bt ∈ R b×n is a zero-padded version of G−1t,Bt to make Eq. (18) hold. Therefore, after one epoch (from the (tk + 1)-th to the ((t+ 1)k)-th update), we have
f(t+1)k − y =  tk∏ t′=(t+1)k−1 ( I− Jt′,t′+1J>t′,Bt′ G̃ −1 t′,Bt′ ) (ftk − y) =: At(ftk − y).
We will see that the matrix At is close to the matrix A defined in (15), so it boils down to analyzing the spectral properties of A.
For both theorems, we can compute that as M increases, the norm of the update ‖Wt −Wt+1‖F does not increase with M , so the update is small compared to the Gaussian initialization where ‖W0‖2 = Θ( √ M). From this we can derive that the matrices J,G etc. remain close to their initialization, which makes bounding their norms possible. The full proof is in the appendix.
In conclusion, the accelerated convergence is related to the local linearity and the stability of the Jacobian and Gram matrix. We emphasize that our theorems serve more as a motivation than a justification of our GGN algorithm, because we expect that GGN works in practice, even under milder situations when M is not as large as the theorem demands or for deep networks with different architectures, and that GGN would still perform much better than first-order methods.

heading: 3.3 ANALYSIS OF THE PER-ITERATION COMPUTATIONAL COMPLEXITY
text: We have proved that for sufficiently overparametrized deep neural networks, full-batch GGN has quadratic convergence rate. In this subsection, we analyze the per-iteration computational cost of GGN, and compare it to that of SGD.
For every mini-batch (i.e., iteration), there are two major steps of computation in GGN:
• (A). Forward, and then backpropagate for computing the Jacobian matrix J. • (B). Use J to compute the update J>(λG + αI)−1(f − y)
We show that the computational complexity of (A) is the same as that of SGD with the same batch size; and the computational complexity of (B) is small compared to (A) for typical networks and batch sizes. Thus, the per-iteration computation overhead of GGN is very small compared to SGD. Overall, in terms of training time, GGN can be much faster than SGD.
For the computation in step (A), the forward part is just the same as that of SGD. For the backward part, for every input data, GGN keeps track of the output’s derivative for the nodes in the middle of the computational graph. This part is just the same as backpropagation in SGD. What is different is that GGN also, for every input data, keeps track of the output’s derivative for the parameters; while in SGD the derivatives for the parameters are averaged over a batch of data. However, it is not difficult to see the computational costs of GGN and SGD are the same.
For the computation in step (B), observe that the size of the Jacobian is b × m where b is the batch size and m is the number of parameters. The Gram matrix Gt,Bt = Jt,BtJ > t,Bt
in our GramGauss-Newton method is of size b× b and it only requires O(b2m+ b3) for computing Gt,Bt and a matrix inverse. Multiplying the two matrices to f − y requires even less computation. Overall, the computational cost in step (B) is small compared to that of step (A).

heading: 4 EXPERIMENTS
text: Given the theoretical findings above, in this section, we compare our proposed GGN algorithm with several baseline algorithms in real applications. In particular, we mainly study two regression tasks, AFAD-LITE (Niu et al., 2016) and RSNA Bone Age (rsn, 2017).

heading: 4.1 EXPERIMENTAL SETTING
text: AFAD-LITE task is to predict the age of human from the facial information. The training data of the AFAD-LITE task contains 60k facial images and the corresponding age for each image. We choose ResNet-32 (He et al., 2016) as the base model architecture. During training, all input images are resized to 64 ∗ 64. We study two variants of the ResNet-32 architecture: ResNet-32 with batch normalization layer (referred to as ResNetBN), and ResNet-32 with Fixup initialization (Zhang et al., 2019b) (referred to as ResNetFixup). In both settings, we use SGD as our baseline algorithm. In particular, we follow Qian (1999) to use its momentum variant and set the hyper-parameters lr=0.003 and momentum=0.9 determined by selecting the best optimization performance using grid search. Since batch normalization is computed over all samples within a mini-batch, it is not consistent with
our assumption in Section 2 that the regression function has the form of f(w,x), which only depends on w and a single input datum x. For this reason, the GGN algorithm does not directly apply to ResNetBN, and we test our proposed algorithm on ResNetFixup only. We set λ = 1 and α = 0.3 for GGN. We follow the common practice to set the batch size to 128 for our proposed method and all baseline algorithms. Mean square loss is used for training.
RSNA Bone Age task is a part of the 2017 Pediatric Bone Age Challenge organized by the Radiological Society of North America (RSNA). It contains 12,611 labeled images. Each image in this dataset is a radiograph of a left hand labeled with the corresponding bone age. During training, all input images are resized to 64 ∗ 64. We also choose ResNetBN and ResNetFixup for this experiment, and use ResNetBN and ResNetFixup trained in the first task as warm-start initialization. We use lr= 0.01 and momentum= 0.9 for SGD, and use λ = 1 and α = 0.1 for GGN. Batch size is set to 128 in these experiments, and mean square loss is used for training.

heading: 4.2 EXPERIMENTAL RESULTS
text: Convergence. The training loss curves of different optimization algorithms for AFAD-LITE and RSNA Bone Age tasks are shown in Figure 1. On both tasks, our proposed method converges much faster than the baselines. We can see from Figure 1a and Figure 1b that, on the AFAD-LITE task, the loss using our GGN method quickly decreases to nearly zero in 30 epochs. On the contrary, for both baselines using SGD, the loss decays much slower than our method in terms of wall clock time and epochs. Similar advantage of GGN can also be observed on the RSNA bone age task.
Generalization performance and different hyper-parameters. We can see that our proposed method trains much faster than other baselines. However, as a machine learning model, generalization performance also needs to be evaluated. Due to space limitation, we only provide the test curve for the RSNA Bone Age task in Figure 2a. From the figure, we can see that the test loss of our proposed method also decreases faster than the baseline methods. Furthermore, the loss of our GGN algorithm is lower than those of the baselines. These results show that the GGN algorithm can not only accelerate the whole training process, but also learn better models.
We then study the effect of hyper-parameters used in the GGN algorithm. We try different λ and α on the RSNA Bone Age task and report the training loss of all experiments at the 10th epoch. All results are plotted in Figure 2c. In the figure, the x-axis is the value of λ and the y-axis is the value of α. The gray value of each point corresponds to the loss, the lighter the color, the higher the loss. We can see that the model converges faster when λ is close to 1. In GGN, α can be considered as the inverse value of the learning rate in SGD. Empirically, we find that the convergence speed of training loss is not that sensitive to α given a proper λ, such as λ = 1. Some training loss curves of different hyper-parameter configurations are shown in Figure 2b.

heading: 5 CONCLUSION AND DISCUSSIONS
text: We propose a novel Gram-Gauss-Newton (GGN) method for solving regression problems with square loss using overparameterized neural networks. Despite being a second-order method, the computation overhead of the GGN algorithm at each iteration is small compared to SGD. We also prove that if the neural network is sufficiently wide, GGN algorithm enjoys a quadratic convergence rate. Experimental results on two regression tasks demonstrate that GGN compares favorably to SGD on these data sets with standard network architectures. Our work illustrates that second-order methods have the potential to compete with first-order methods for learning deep neural networks with huge number of parameters.
In this paper, we mainly focus on the regression task, but our method can be easily generalized to other tasks such as classification as well. Consider the k-category classification problem, the neural network outputs a vector with k entries. Although this will increase the computational complexity of getting the Jacobian whose size increases k times, i.e., J ∈ R(bk)×m, each row of J can be still computed in parallel, which means the extra cost only comes from parallel computation overhead when we calculate in a fully parallel setting. While most first-order methods for training neural networks can hardly make use of the computational resource in parallel or distributed settings to accelerate training, our GGN method can exploit this ability. For first-order methods, basically extra computational resource can only be used to calculate more gradients at a time by increasing batch size, which harms generalization a lot. But for GGN, more resource can be used to refine the gradients and achieve accelerated convergence speed with the help of second-order information. It is an important future work to study the application of GGN to classification problems.

heading: A SOME PREPARATIONS FOR THE PROOF OF THEOREMS
text: Notations. We use the following notations for the rest of the sections.
• Let [n] = {1, · · · , n}.
• JW,x ∈ RM×d denotes the gradient ∂f(W,x)∂W , which is of the same size of W. • The bold JW or J(W) denotes the Jacobian with regard to all n data, with each gradient
for W vectorized, i.e.
JW = vec(JW,x1) >
... vec(JW,xn) >  ∈ Rn×(Md). (19) • wr denotes the r-th row of W, which is the incoming weights of the r-th neuron. • W0 stands for the parameters at initialization. • dW,x := σ′(Wx) ∈ RM×1 denotes the (entry-wise) derivative of the activation function. • We use 〈·, ·〉 to denote the inner product, ‖·‖2 to denote the Euclidean norm for vectors or
the spectral norm for matrices, and ‖·‖F to denote the Frobenius norm for matrices.
We can easily derive the formula for J as
JW,x = 1√ M ((dW,x ◦ a)x>), (20)
where ◦ is the point-wise product. So we can also easily solve G as
Gij = 〈Jxi , Jxj 〉 = 1
M M∑ r=1 x>i xjσ ′(wrxi)σ ′(wrxj). (21)
Our analysis is based on the fact that G stays not too far from its infinite-width limit K(xi,xj) = Ew∼N (0,Id) ( x>i xjσ ′(wxi)σ ′(wxj) ) ,
which is a positive definite matrix with least eigenvalue denoted λ0, and we assume λ0 > 0. λ0 is a small data-dependent constant, and without loss of generality we assume λ0 ≤ 1, or else we can just take λ0 = 1 if λmin(K) > 1.
The first lemma is about the estimation of relevant norms at initialization. Lemma 2 (Bounds on Norms at Initialization). If M = Ω (d log(16n/δ)), then with probability at least 1− δ/2 the following holds (a). ‖W0‖2 = O( √ M).
(b). f(W0,xi) = O(1), for i ∈ [n]. (c). ‖JW0,xi‖F = O(1), for i ∈ [n].
Proof. (a). The lemma is a well-known result concerning the estimation of singular values of Gaussian random matrices (see Corollary 5.35 of Vershynin (2010)). Notice that W0 ∈ RM×d is a Gaussian random matrix, the Corollary states that with probability 1− 2e−t2/2 one has
‖W0‖2 ≤ √ M + √ d+ t.
By choosing M = max(d, √ 2 log(8/δ)), we obtain ‖W0‖2 ≤ 3 √ M with probability 1− δ/4.
(b). First, ar, r ∈ [M ] are Rademacher variables, thereby 1-sub-Gaussian, so with probability 1− 2e−Mt2/2 we have 1M ∑M r=1 ar ≤ t. This means if we take M = Ω (log(16/δ)),
Pr[ 1√ M M∑ r=1 ar = O(1)] ≥ 1− δ 8 . (22)
Next, the vector vi = 1‖xi‖2W0xi ∈ R M×1 is a standard Gaussian vector. Suppose the activation σ(·) is l-Lipschitz and l is O(1) by our assumption, with the vector a fixed, the function
φ : RM → R,vi 7→ 1√ M a>σ(‖xi‖2 vi) = f(W0,xi)
has a Lipschitz parameter of l ‖xi‖2 / √ M = O(1/ √ M). According to the classic result on the concentration of a Lipschitz function over Gaussian variables (see Theorem 2.26 of Wainwright (2019)), we have
Pr[|φ(vi)− EW0(φ(vi))| ≥ t] ≤ 2 exp
( − Mt 2
2l2 ‖xi‖2
) ,
which means if M = Ω(log(16n/δ)),∣∣∣∣∣ 1√M a>σ(W0xi)− 1√M ( M∑ r=1 ar)Ew∼N (0,Id)[σ(wxi)] ∣∣∣∣∣ = O(1) (23) holds jointly for all i ∈ [n] with probability 1− δ/8. Note that∣∣Ew∼N (0,Id)[σ(wxi)]∣∣ ≤ |σ(0)|+ l · Eξ∼N (0,‖xi‖22))[|ξ|] = O(1). (24)
Plugging in (22) and (24) into (23), we see that as long as M = Ω(log(16n/δ)), then with probability 1− δ/4, for all i ∈ [n],
f(W0,xi) = 1√ M a>σ(W0xi) = O(1).
(c). Since σ is O(1)-Lipschitz, we have ‖dW,xi‖∞ = O(1). According to (20) we can easily know that ‖JW,xi‖F ≤ 1√ M ‖Diag(d)‖2 ‖a‖2 ‖x‖2 = O(1).
The next lemma is about the least eigenvalue of the Gram matrix G at initialization. It shows that when M is large, GW0 is close to K and has a lower bounded least eigenvalue. It is the same as Lemma 3.1 in Du et al. (2018b), but here we restate it and its proof for the reader’s convenience.
Lemma 3 (Bound on the Least Eigenvalue of the Gram Matrix at Initialization). If the width M = Ω
( n2 log(2n/δ)
λ20
) , then with probability at least 1− δ/2 over random initialization, we have
λmin(GW0) ≥ 3
4 λ0.
Proof. Because σ is Lipschitz, σ′(wxi)σ′(wxj) is bounded by O(1). For every fixed (i, j) pair, at initialization Gij is an average of independent random variables, and by Hoeffding’s inequality, applying union bound for all n2 of (i, j) pairs, with probability 1− δ/2 at initialization we have
|Gij −Kij | ≤ O
(√ log(2n/δ)
M ) and then
‖GW0 −K‖ 2 2 ≤ ‖GW0 −K‖ 2 F = O
( n2 log(2n/δ)
M
) .
Thus ifM = Ω ( n2 log(2n/δ)
λ20
) we can have ‖GW0 −K‖2 ≤ 1 4λ0 and thus λmin(GW0) ≥ 3 4λ0.
Next, in Lemma 4 and 5 we will bound the relevant norms and the least eigenvalue of G inside some scope of W that covers the whole optimization trajectory starting from W0. Specifically, we consider the range
B(R) , {W : ‖W −W0‖F ≤ R} ,
where R is determined later to make sure that the optimization trajectory remains inside B(R). The idea of the whole convergence theorem is that when the width M is large, R is very small compared to its initialization scale: ‖W0‖2 = O( √ M). This way, neither the Jacobian nor the Gram matrix changes much during optimization.
Lemma 4 (Bounds on Norms in the Optimization Scope). Suppose the events in Lemma 2 hold. There exists a constant C > 0 such that if M ≥ CR2, we have the following: (a) For any W ∈ B(R), we have
‖W‖2 = O( √ M). (25)
(b) For any W1,W2 ∈ B(R), if ‖W1 −W2‖F ≤ R′, then we have
‖JW1,xi − JW2,xi‖F = O ( R′√ M ) , ∀i, and ‖JW1 − JW2‖2 = O ( R′ √ n M ) . (26)
Also, for any W ∈ B(R), we have
‖JW‖F = O(1) and ‖JW‖2 = O( √ n). (27)
Proof. (a). This is straightforward from Lemma 2(a), the definition of B(R), and M = Ω(R2).
(b). According to the O(1)-smoothness of the activation, we have
‖dW1,xi − dW2,xi‖2 ≤ O(1) · ‖W1xi −W2xi‖2 = O(R ′),
so we can bound
‖JW1,xi − JW2,xi‖F
= 1√ M ∥∥Diag(a)(dW1,xi − dW2,xi)x>∥∥F ≤ 1√
M ‖Diag(a)‖2 ‖dW1,xi − dW2,xi‖2 ‖xi‖2
≤O( R ′
√ M ).
And according to (19), we have
‖JW1 − JW2‖2 ≤ ‖JW1 − JW2‖F
= √√√√ n∑ i=1 ‖JW1,xi − JW2,xi‖ 2 F
= O ( R′ √ n
M
) .
Also, taking W1 = W and W2 = W0, combining with Lemma 2(c), we see there exists C such that for M ≥ CR2 we have ‖JW‖F = O(1), and naturally ‖JW‖2 = O( √ n). Note: The constants hidden in the O(·) notation are irrelevant with C.
The next Lemma deals with the Gram matrix within B(R), ensuring a lower bound on the least eigenvalue throughout the optimization process (Du et al., 2018a). Lemma 5 (Least Eigenvalue in the Optimization Scope). For W ∈ B(R), suppose the events in Lemma 2 and 3 hold, there exists a constant C ′ such that for M ≥ C ′n2R2
λ20 , we have
‖GW −GW0‖2 ≤ λ0 4 ,
and thus combined with Lemma 3, we know that GW remains invertible when W ∈ B(R) and satisfies ∥∥G−1W∥∥2 ≤ 2λ0 . Proof. Based on the results in Lemma 4(b), we have
‖GW −GW0‖2 = ∥∥JWJ>W − JW0J>W0∥∥2 (28)
≤ ∥∥(JW − JW0)J>W)∥∥2 + ∥∥JW0(J>W − J>W0)∥∥2 (29)
≤ O( nR√ M ). (30)
To make the above less than λ04 , choosing M greater than some C′n2R2
λ20 suffices, and the lemma is
thus proved. Again the constants hidden in the O(·) notation are irrelevant with C ′.

heading: B PROOF OF THEOREM 1
text: Proof idea. In this section, we use Wt, t ∈ {0, 1, · · · } to represent the parameter W after t iterations. For convenience, Jt,Gt, ft is short for JWt ,GWt , fWt respectively. We introduce
Jt,t+1 = ∫ 1 0 J((1− s)Wt + sWt+1)ds. (31)
For each iteration, if Gt is invertible, we have
ft+1 − y = ft − y + Jt,t+1 vec(Wt+1 −Wt) = ft − y − Jt,t+1J>t G−1t (ft − y) = (Jt − Jt,t+1)J>t G−1t (ft − y),
hence
‖ft+1 − y‖2 ≤ ‖Jt − Jt,t+1‖2 ∥∥J>t ∥∥2 ∥∥G−1t ∥∥2 ‖ft − y‖2 . (32)
Then we control the first term of the right hand side based on Lemma 4 in the following form
‖Jt − Jt,t+1‖2 ≤ Ĉ√ M ‖Wt −Wt+1‖2
= Ĉ√ M ∥∥J>t G−1t (ft − y)∥∥2 ≤ Ĉ√
M ∥∥J>t ∥∥2 ∥∥G−1t ∥∥2 ‖ft − y‖2 , and plugging into (32) along with norm bounds on J and G we obtain a second-order convergence. Formal proof. Let Rt = ‖Wt −Wt+1‖F for t ∈ {0, 1, · · · }. We take R = Θ ( n λ0 ) in Lemma 4 and 5 (the constant is chosen to make the right hand side of (34) hold). We prove that there exists an M = Ω ( max ( n4
λ40 , n 2d log(16n/δ) λ20
)) (with enough constant) that suffices. First we can easily verify
that all the requirements for M in Lemma 2-5 can be satisfied . Hence, with probability at least 1− δ all the events in Lemma 2-5 hold. Under this situation, we do induction on t to show the following:
• (a). Wt ∈ B(R).
• (b). If t > 0, then ‖ft − y‖2 ≤ n3/2 λ20 √ M ‖ft−1 − y‖22
As long as (b) is true for all t, then choosing M large enough to make sure the series {‖ft − y‖2}∞t=0 converges to zero, we obtain the second-order convergence property.
For t = 0, (a) and (b) hold by definition. Suppose the proposition holds for t = 0, · · · , T . Then for t = 0, · · · , T , Gt is invertible. Recall that the update rule is vec(Wt+1) = vec(Wt)−J>t G−1t (ft− y), we have
Rt = ‖Wt −Wt+1‖F = ‖vec(Wt)− vec(Wt+1)‖2 ≤ ∥∥J>t ∥∥2 ∥∥G−1t ∥∥2 ‖ft − y‖2
≤ O (√ n
λ0 ‖ft − y‖2
) . (Lemma 4 and 5) (33)
According to Lemma 2(b) and the assumption that the target label yi = O(1), we have ‖f0 − y‖22 = O(n). When T > 1, the decay ratio at the first step is bounded as
‖f1 − y‖2 ‖f0 − y‖2 ≤ O
( n 3 2
λ20 √ M
) ‖f0 − y‖2 ≤ O ( n2
λ20 √ M
) =: r,
and taking M = Ω(n 4
λ40 ) with enough constant can make sure r is a constant less than 1, in which case
the second-order convergence property (b) will ensure a faster ratio of decay at each subsequent step in ‖ft − y‖Tt=0. Combining (33), we have
T∑ t=0 Rt ≤ O ( T∑ t=0 √ n λ0 ‖ft − y‖2 ) ≤ O ( n λ0 ∞∑ t=1 rt−1 ) = O ( n λ0 ) ≤ R. (34)
Therefore, WT+1 ∈ B(R), and (34) also holds when T = 0, so (a) is true. Since B(R) is convex, this means we also have sWT + (1− s)WT+1 ∈ B(R) for s ∈ [0, 1]. Hence we can bound the difference of the Jacobian as
‖JT,T+1 − JT ‖2 ≤ ∫ 1 0 ‖J(sWT + (1− s)WT+1)− J(WT )‖2 ds
≤ O (∫ 1
0
sRt
√ n
M ds
) (Lemma 4(b))
= O
( n
λ0 √ M ‖fT − y‖2
) . (Using (33)) (35)
So we use the bound of (32) and obtain ‖fT+1 − y‖2 ≤ ‖JT − JT,T+1‖2 ∥∥J>T ∥∥2 ∥∥G−1T ∥∥2 ‖fT − y‖
≤ n 3/2
λ20 √ M ‖fT − y‖22 , (Using (35) and Lemma 4, 5)
which proves (b). This concludes our proof.

heading: C PROOF OF THEOREM 2
text: In this section, we give the proof of the convergence of our mini-batch GGN algorithm (Theorem 2).
Some additional notations. For the convenience of our proof, we will use slightly different notations than that in Section 3.2. Let n = bk, where b is the batch size, k is the number of batches, or equivalently, number of updates in each epoch. For epoch t ∈ {1, 2, · · · } and batch i ∈ [k], we will use Wti ∈ Rd×M to denote the parameters before the i-th update in epoch t. Let Jti = J(Wti) ∈ Rn×Md, Gti = G(Wti) ∈ Rn×n, and fti = f(Wti) ∈ Rn×1. (All of W0,W1,W11,G0, etc. can represent initial values in the proof.) Also, let Wt = Wt1,Jt = Gt1,Gt = Jt1, ft = ft1, and let Wt(k+1) = W(t+1)1,Jt(k+1) = J(t+1)1, etc. In the i-th batch, we use the data (xl, yl)ibl=(i−1)b+1. To specify the matrix related to a batch, we use the following indexing method: For i′, i′′ ∈ [k], let
Jti,i′ = Jti((i ′ − 1)b+ 1 : i′b, 1 : Md) ∈ Rb×(Md),
which is the b rows of Jti that represent the Jacobian of the i-th batch. Similarly, we have
Gti,i′i′′ = Jti,i′J > ti,i′′ ∈ Rb×b,
fti,i′ = fti((i ′ − 1)b+ 1 : ib) ∈ Rb×1,
yi′ = y((i′−1)b+1:ib) ∈ Rb×1.
Similar to (31) in the proof of Theorem 1, we make use of the notation
Jt(i,i+1) = ∫ 1 0 J(sWti + (1− s)Wt(i+1))ds, (36)
and similarly define Jt(i,i+1),i′ ∈ Rb×Md, etc. In addition, we make use of the matrix
G̃ti = [ 0b×(i−1)b G −1 ti,ii 0b×(k−i)b ] ∈ Rb×n
in the formula of the update.
Similar to the proof of Theorem 1, the idea of our proof is that due to the fact that the change of W is small compared to the scale of its initialization, the matrix J and G remains stable. Informally, this makes the algorithm close to solving kernel regression iteratively by batches, which is equivalent to solving a system of linear equations
JJ>r = f(W0)− y, (37)
for variables r ∈ Rn×1 (where vec(W) = vec(W0) + J>r), using the Gauss-Siedel method, which means solving
JJ>G  r(old)(1 : (i− 1)b)r(new)((i− 1)b+ 1 : ib) r(old)((ib+ 1 : kb))  = (f(W0 − y))(i−1)b+1:ib for the i-th batch in every epoch. Therefore, it is natural that the matrix A = L>(D− L)−1 in (15) is introduced. We will show later that the real update follows
ft+1 − yt = At(ft − y)
for some At ≈ A. In order to prove the theorem, we need some additional lemmas as follows. Lemma 6 (Formula for the Update). If the Gram matrix at each step is invertible, then:
(a) The update of W is
vec(Wt(i+1)) = vec(Wti)− J>ti,iG−1ti,ii(fti,i − yi) (38)
= vec(Wti)− J>ti,iG̃ti(ft − y) (39)
(b) The formula for f − y is
ft(i+1) − y = (In − Jt(i,i+1)J>ti,iG̃ti)(fti − y) (40)
(c) The update of f − y satisfies
ft(i+1) − y = Uti(Dt − Lt)−1(ft − y), (41)
where
Dt =  Gt1,11 0 · · · 0 0 Gt2,22 · · · 0 ... ... . . .
... 0 0 · · · Gtk,kk
 , Lt = −  0 0 · · · 0 Jt(1,2),2J > t1,1 0 · · · 0 ... ... . . . ...
Jt(1,2),kJ > t1,1 Jt(2,3),kJ > t2,2 · · · 0  , and
Uti = Diag −  Jt(1,2),1J > t1,1 −Gt1,11 Jt(2,3),1J>t2,2 · · · Jt(i,i+1),1J>ti,i 0 Jt(2,3),2J > t2,2 −Gt2,22 · · · Jt(i,i+1),2J>ti,i ... ... . . . ...
0 0 · · · Jt(i,i+1),iJ>ti,i −Gti,ii  ib×ib ,
 Gt(i+1),(i+1)(i+1) 0 · · · 0 Jt(i+1,i+2),i+2J > t(i+1),i+1 Gt(i+2),(i+2)(i+2) · · · 0 ... ... . . . ...
Jt(i+1,i+2),kJ > t(i+1),i+1 Jt(i+2,i+3),kJ > t(i+2),i+2 · · · Gtk,kk  (k−i)b×(k−i)b  . Or, if we denote
Uti =
[ Ũti (∈ Rib×ib) 0
0 Uti(∈ R(k−i)b×(k−i)b)
] ,
Dti − Lti = [ D̃ti − L̃ti (∈ Rib×ib) 0
D̂ti − L̂ti Dti − Lti(∈ R(k−i)b×(k−i)b)
] ,
then we have
ft(i+1) − y = [
Ũti(D̃ti − L̃ti)−1 0 −(D̂ti − L̂ti)(D̃ti − L̃ti)−1 I(k−i)b
] (ft − y) (42)
Specifically, we have ft+1 − y = Ut(Dt − Lt)−1(ft − y) =: At(ft − y), (43)
where
Ut = −  Jt(1,2),1J > t1,1 −Gt1,11 Jt(2,3),1J>t2,2 · · · Jt(k,k+1),1J>tk,k 0 Jt(2,3),2J > t2,2 −Gt2,22 · · · Jt(k,k+1),2J>tk,k
... ...
. . . ...
0 0 · · · Jt(k,k+1),kJ>tk,k −Gtk,kk  . Proof. (a) For (38), this is exactly the update formula (10) for the i-th batch where the Jacobian and Gram matrix are Jti,i and Gti,ii respectively. Note that
G−1ti,ii(fti,i − yi) = [ 0b×(i−1)b G −1 ti,ii 0b×(k−i)b ] (fti − y),
we then obtain (39) from (38).
(b). We have (ft(i+1) − y)− (fti − y)
= ft(i+1) − fti = Jt(i,i+1)(vec(Wt(i+1))− vec(Wti)) = Jt(i,i+1)J > ti,iG̃ti(fti − y) (By formula (39)),
we obtain (40).
(c). Based on (40) we know that
ft(i+1) − y = 1∏ i′=i (In − Jt(i′,i′+1)J>ti′,i′G̃ti′)(fti − y),
where the index goes in decreasing order from left to right. So in order to prove (41) we only need to prove that [
1∏ i′=i (In − Jt(i′,i′+1)J>ti′,i′G̃ti′)
] (Dt − Lt) = Uti, (44)
which we will prove by induction on i. For i = 0 it is trivial that Dt − Lt = Ut0 by definition. Suppose (44) holds for i, then[
1∏ i′=i+1 (In − Jt(i′,i′+1)J>ti′,i′G̃ti′)
] (Dt − Lt)
= (In − Jt(i+1,i+2)J>t(i+1),i+1G̃t(i+1))Uti = Uti − Jt(i+1,i+2)J>t(i+1),i+1G −1 t(i+1),(i+1)(i+1)Uti((ib+ 1 : (i+ 1)b, 1 : n))
= Uti −  Jt(i+1,i+2),1J > t(i+1),i+1 Jt(i+1,i+2),2J > t(i+1),i+1
... Jt(i+1,i+2),kJ > t(i+1),i+1  [0b×(ib) Ib×b 0b×(k−i−2)b] = Ut(i+1),
which proves (41). Note that by the definition, we have Uti = Dti − Lti, we can then obtain (42) by Uti(Dt − Lt)−1 = [ Ũti 0 0 Uti ] [ D̃ti − L̃ti 0 D̂ti − L̂ti Dti − Lti ]−1 = [ Ũti 0 0 Uti ] [ (D̃ti − L̃ti)−1 0 −(Dti − Lti)−1(D̂ti − L̂ti)(D̃ti − L̃ti)−1 (Dti − Lti)−1
] = [ Ũti(D̃ti − L̃ti)−1 0
−(D̂ti − L̂ti)(D̃ti − L̃ti)−1 I
] .
By (43), we can see that the iteration matrix At is close to the matrix A defined in (15). The convergence of the algorithm is much related to the eigenvalues of A. In the next two lemmas, we bound the spectral radius of A and provide an auxiliary result on the convergence on perturbed matrices based on the spectral radius.
Lemma 7. Suppose the least eigenvalue of the initial Gram matrix G0 satisfies λmin(G0) ≥ 34λ0, (which is true with probability 1− δ if M = Ω
( n2 log(2n/δ)
λ20
) , according to Lemma 3). Also assume
‖JW0,xl‖F = O(1) for any l ∈ [n] (which is true with probability 1 − δ, according to Lemma 2). Then the spectral radius A, or equivalently, maximum norm of eigenvalues of A, denoted ρ(A), satisfies
ρ(A) ≤ 1− Ω ( λ20 n2 ) . (45)
Proof. For any eigenvalue λ ∈ C of A, it is an eigenvalue of A>, so there exists a corresponding unit vector v ∈ Cn×1 such that λv = A>v = (D− L>)−1Lv, which means
λ(D− L>)v − Lv = 0. (46)
Since G0 = D− L− L> is positive definite with eigenvalues at least 34λ0, we have v ∗(D− L− L>)v ≥ 34λ0, which means
d− 2 Re(l) ≥ 3 4 λ0, (47)
where d = v∗Dv, l = v∗Lv, and v∗L>v = l. Let vi = v((i − 1)b + 1 : ib) ∈ Rb×1 be the i-th batch of the vector v, and ṽi = [01×(i−1)b,v>i ,01×(k−i)b] > ∈ Rn×1. It is not hard to see that
d = v∗Dv = k∑ i=1 v∗iG0,iivi = k∑ i=1 ṽ∗iG0ṽi ≥ k∑ i=1 3 4 λ0 ‖ṽi‖22 = 3 4 λ0. (48)
Also, since by our assumption each entry of L (or of G0) is at most O(1), we have
|l| = v∗Lv ≤ ‖L‖2 ≤ ‖L‖F = O(n). (49)
Now we use take an inner product of v with (46) and get
0 = v∗(λ(Dv − L>v)− Lv) = λ(d− l) + l, (50)
and therefore
|λ| = ∣∣∣∣ ld− l ∣∣∣∣ (By solving (50)) = √ |l|2
(d− l)(d− l) =
√ 1
d(d−2Re(l)) |l|2 + 1 ≤ √√√√ 1
Ω ( λ20 n2 ) + 1
(Using (47), (48), (49))
= 1− Ω ( λ20 n2 ) .
This concludes the proof for this lemma. Lemma 8. Denote ρ(A) = ρ0 ≤ 1. Let A be diagonalized as A = P−1QP and µ = ‖P‖2 ∥∥P−1∥∥ 2 (see Assumption 2). Suppose we have ‖At −A0‖2 ≤ δ for t ≤ T , then∥∥∥∥∥ T∏ i=1 At ∥∥∥∥∥ 2 ≤ µ(ρ+ µδ)T
Proof. We have∥∥∥∥∥ T∏ i=1 At ∥∥∥∥∥ 2 = ∥∥∥∥∥ T∏ i=1 (A + (At −A0)) ∥∥∥∥∥ 2
= ∥∥∥∥∥ T∏ i=1 P−1 ( Q + P(At −A0)P−1 ) P ∥∥∥∥∥ 2
= ∥∥∥∥∥P ( T∏ i=1 ( Q + P(At −A0)P−1 )) P ∥∥∥∥∥ 2
≤ ‖P‖2 ( T∏ t=1 (‖Q‖2 + ‖P‖2 ‖At −A0‖2 ‖P‖2) )∥∥P−1∥∥ 2
≤ µ(ρ+ µδ)T .
In addition to the bounds used in Theorem 1, we provide the next lemma with useful bounds of the norms and eigenvalues of relevant matrices in the optimization scope
B(R) = {W : ‖W −W0‖F ≤ R} . Lemma 9 (Relevant Bounds for the Matrices in the Optimization Scope). We assume the events in Lemma 2 and 3 hold, and let M = Cµ 2n8R2
λ80 for some large enough constant C, then:
For any (t, i)-pair, we assume Wt′i′ ∈ B(R) for all i′ ∈ [k] when t′ ≤ t and i′ ∈ [i] when t′ = t in the following propositions.
(a). minv∈Rn,‖v‖2=1 ‖(D− L)v‖2 ≥ 3 8λ0.
(b). Suppose up to t, Wt is in B(R) (which means for i′ ∈ [k + 1], and t′ ∈ [t− 1], Wti ∈ B(R)). Then ‖At −A‖2 ≤ 1−ρ(A) 2µ .
Proof. (a). Because for ‖v‖2 = 1,
‖(D− L)v‖2 ≥ v >(D− L)v
= 1
2 v>Dv +
1 2 v>(D− L− L>)v
≥ 0 + 1 2 λmin(G0) ≥ 3 8 λ0. (By the positive-definiteness of G0 and Lemma 3)
(b). By Lemma 4 we know that within B(R), the size of the Jacobian ‖Jxl‖F w.r.t. data xl is O(1), and the difference ‖JW1,xl − JW2,xl‖F within O(R) is bounded by O ( R√ M ) , this can be applied to each entry of D,L,U, etc., including those Jt′(i′,i′+1) terms by the convexity of B(R), and we
can see that each entry of these matrices has size at most O(1) and varies inside an O (
R√ M
) range.
Therefore we get
‖Ut‖F = O(n), ‖Dt − Lt‖F = O(n),
and
‖Ut −U‖2 = O ( Rn√ M ) ,
‖(Dt − Lt)− (D− L)‖2 = O ( Rn√ M ) .
By our chosen M , O ( Rn√ M ) ≤ 18λ0 can definitely hold, so along with (a) we have∥∥(D− L)−1∥∥
2 = O
( 1
λ0 ) , ∥∥(Dt − Lt)−1∥∥2 = O( 1λ0 ) .
Hence, ‖At −A‖2 = ∥∥(Ut(Dt − Lt)−1 −U(D− L)−1∥∥2
= ∥∥(Ut −U)(Dt − Lt)−1 + U((Dt − Lt)−1((Dt − Lt)− (D− L))(D− L)−1∥∥2
≤ O ( Rn√ M ) O ( 1 λ0 ) +O (n)O ( 1 λ0 ) O ( Rn√ M ) O ( 1 λ0 ) = O ( n2R
λ20 √ M
) ,
and along with Lemma 7, in order to have ‖At −A‖2 ≤ 1−ρ(A) 2µ , all we need is
O
( n2R
λ20 √ M
) ≤ Ω ( λ20 n2 ) ,
so choosing some M = Cµ 2n8R2
λ80 suffices.
With all the preparations, now we are ready to prove Theorem 2. The logic of the proof is just the same as what we did in the formal proof of Theorem 1, where we then used induction on t and now we use induction on the pair (t, i). The key, is still selecting some R so that in each step Wti remains in B(R). Combined with the previous Lemmas, in B(R) we have At being close to A, and then convergence is guaranteed. Formal proof of Theorem 2. Let Rti = ∥∥Wt(i+1) −Wti∥∥F . We take R = Θ ( n5
λ40 ) in the range B(R) (where the constant is chosen to make the right hand side of (52) hold). We prove that there exists an
M = max ( Ω ( µ2n18
λ160
) ,Ω ( n2d log(16n/δ)
λ20 )) with a large enough constant that suffices. First we can easily verify that all the requirements for M in Lemma 2-9 (most importantly, Lemma 9) can be satisfied. Hence, with probability at least 1− δ all the events in Lemma 2-9 hold. Under this situation, we do induction on (t, i) (t ∈ {1, 2, · · · }, i ∈ [k], by dictionary order) to show that:
• Wti ∈ B(R).
For (t, i) = (1, 1), it holds by definition.
Suppose the proposition holds up to (t, i). Then since Wti ∈ B(R), by Lemma 4 we know that λmin(Gti) ≥ λ02 . This naturally gives us λmin(Gti,ii) ≥ λ0 2 , which means Gti,ii is invertible.
Similar to the argument in the proof of Lemma 9, we know that each entry of Jti,Jt(i,i+1),Dti,Lti,Uti, etc., whose index of (t′, i′) only contains i′ with i′ ≤ i, is of scale O(1) and varies at most O ( R√ M ) , which is a result of Lemma 4. Then we know
‖Jti,i‖F = O( √ b) = O( √ n),∥∥∥D̃t(i−1) − L̃t(i−1)∥∥∥
F = O(n),∥∥∥D̂t(i−1) − L̂t(i−1)∥∥∥
F = O(n),
etc. We then also know∥∥∥(D̃t(i−1) − L̃t(i−1))− (Dt − Lt)(1:(i−1)b,1:(i−1)b)∥∥∥ F ≤ O ( nR√ M ) .
Since we can also have ∥∥((Dt − Lt)− (D− L))(1:(i−1)b,1:(i−1)b)∥∥F ≤ O ( nR√M ), and given our choice of M and R the right hand side is less than 18λ0, along with Lemma (9) (which says
min v∈Rn,‖v‖2=1
‖(D− L)v‖2 ≥ 3
8 λ0
and naturally, since this matrix is block-lower-triangular,
min v∈R(i−1)b,‖v‖2=1 ∥∥(D− L)1:(i−1)b,1:(i−1)bv∥∥2 ≥ 38λ0 holds) we know that D̃t − L̃t is invertible and has the bound ∥∥∥(D̃t − L̃t)−1∥∥∥ 2 ≤ 4λ0 .
Based on the update rule (38), we have Rti = ∥∥Wti −Wt(i+1)∥∥F
= ‖vec(Wt)− vec(Wt+1)‖2 ≤ ∥∥J>ti,iG−1ti,ii(fti,i − yi)∥∥2
≤ O (√ n
λ0
) ‖fti − y‖2
(∥∥J>ti,i∥∥2 = O(√n),∥∥G−1ti,ii∥∥2 = O ( 1
λ0 )) = O (√ n
λ0 )∥∥∥∥[ Ũt(i−1)(D̃t(i−1) − L̃t(i−1))−1 0−(D̂t(i−1) − L̂t(i−1))(D̃t(i−1) − L̃t(i−1))−1 I(k−i+1)b ]∥∥∥∥
2
‖ft − y‖2
(By formula (42)) ≤ O (√ n
λ0
)(∥∥∥Ũt(i−1)(D̃t(i−1) − L̃t(i−1))−1∥∥∥ F
+ ∥∥∥(D̂t(i−1) − L̂t(i−1))(D̃t(i−1) − L̃t(i−1))−1∥∥∥
F + ∥∥I(k−i+1)b∥∥F ) ‖ft − y‖2 ≤ O ( n3/2
λ20
) ‖ft − y‖2 (51)
≤ O ( n3/2
λ20 )∥∥∥∥∥ t−1∏ i=1 At ∥∥∥∥∥ 2 ‖f1 − y‖2 (By Eq. (43))
≤ O ( n3/2
λ20
) · µ ( 1− Ω ( λ20 n2 ))(t−1) ‖(f1 − y)‖2 (Lemma 7, 8, and 9(b).)
≤ O ( µn2
λ20
)( 1− Ω ( λ20 n2 ))(t−1) , (Lemma 2)
where Eq. (51) uses the following bounds proved above:∥∥∥Ũt(i−1)(D̃t(i−1) − L̃t(i−1))−1∥∥∥ F ≤ ∥∥∥Ũt(i−1)∥∥∥ F ∥∥∥(D̃t(i−1) − L̃t(i−1))−1∥∥∥ 2 = O ( n
λ0 ) ,∥∥∥(D̂t(i−1) − L̂t(i−1))(D̃t(i−1) − L̃t(i−1))−1∥∥∥ F ≤ ∥∥∥(D̂t(i−1) − L̂t(i−1))∥∥∥ F ∥∥∥(D̃t(i−1) − L̃t(i−1))−1∥∥∥ 2 = O ( n
λ0 ) ,∥∥I(k−i+1)b∥∥F = O(n).
Since this also hold for previous (t, i) pairs, we have
∑ (t′,i′)≤(t,i) Rt′i′ ≤ O ( µkn2 λ20 ) t∑ t′=1 ( 1− Ω ( λ20 n2 ))(t′−1) ≤ O ( µkn2
λ20
) O ( n2
λ20 ) ≤ O ( n5
λ40
) ≤ R, (52)
which is the reason why we need to take R = Θ ( n5
λ40
) . This means that Wti ∈ B(R) holds. And by
induction, we have proved that W remains in B(R) throughout the optimization process.
The last thing to do is to bound ft − y. By the same logic from above, we have
‖ft − y‖2 ≤ ∥∥∥∥∥ t−1∏ i=1 At ∥∥∥∥∥ 2 ‖f1 − y‖2
≤ µ √ n ( 1− Ω ( λ20 n2 ))(t−1) ,
which proves our theorem.

heading: D ADDITIONAL EXPERIMENTAL RESULTS
text: In this section, we give test performance curve of AFAD-LITE dataset in Figure 3 under the same setting with Section 4.
In addition, we provide more baseline results, e.g. Adam (Kingma & Ba, 2014) and K-FAC (Martens & Grosse, 2015) on RSNA Bone Age dataset. Since we find that, as another alternation of BN, Group Normalization (GN) (Wu & He, 2018) can largely improve the performance of Adam, we also implement the GN layer for our GGN method. We use grid search to obtain approximately best hyper-parameters for every experiment. All experiments are performed with batch size 128, input size 64*64 and weight decay 10−4. We set the number of groups to 8 for ResNetGN. Other hyper-parameters are listed below.
• SGD+ResNetBN: learning rate 0.01, momentum 0.9. • Adam+ResNetBN: learning rate 0.001. • SGD+ResNetGN: learning rate 0.002, momentum 0.9. • Adam+ResNetGN: learning rate 0.0005. • K-FAC+ResNet: learning rate 0.02, momentum 0.9, = 0.1, update frequency 100. • GGN+ResNetGN: λ = 1, α = 0.075.
The convergence results are summarized in Figure 4. Note to make comparison clearer, we use logarithmic scale for training curves.

Review summaries:
The authors propose a scalable second order method for optimization using a quadratic loss. The method is inspired by the Neural Tangent kernel approach, which also allows them to provide global convergence rates for GD and batch SGD. The algorithm has a computational complexity that is linear in the number of parameters and requires to solve a system of the size of the minibatch.  They also show experimentally the advantage of using their proposed methods over SGD. 

Authors propose minimizing neural network using kernel ridge regression. (Formula 9 and Algorithm 1). Main difference of this method is compared to Gauss-Newton, is that it uses JJ' as curvature, which has dimensions b-by-by (batch size b), instead of J'J as curvature, which has dimensions m-by-m (number of parameters m).