A) Samples for training

Let the training set be denoted by $S = \{ z(n) \}_{n=1}^N$, where each training example $z(n)$ consists of an input $x(n)$ and the corresponding target $y(n)$. An example $z(n) = ( x(n), y(n) )$ is independently identically chosen from a probability distribution $\mathcal {D}$ on the joint space $\mathcal {Z}$ of the input space $\mathcal {X}$ and the output space $\mathcal {Y}$. Note that the indices of the examples are omitted if they are clear from the context.

B) Training of deep neural network

The DNN $f: \mathcal {X} \times \Theta \rightarrow \mathcal {Y}$ with parameters $\theta \in \Theta$ predicts the corresponding target $y \in \mathcal {Y}$ for a given input $x \in \mathcal {X}$, where Θ is the parameter space. Its performance is measured on the basis of the expected risk,

(1)\begin{equation} R(\theta) = \mathbb{E}_{z} \left[ \ell (z, \theta) \right], \tag{1}\end{equation}

where $\ell (z, \theta ) = \frac {1}{2} \left\Vert f(x;\theta ) - y \right\Vert ^2$ is the squared loss. The parameters are trained by the gradient descent (GD),

(2)\begin{equation} \theta_{t+1} = \theta_t - \eta \nabla_{\theta} R_S(\theta_t), \tag{2}\end{equation}

to minimize the empirical risk,

(3)\begin{equation} R_S(\theta) = \frac{1}{\vert S \vert} \sum_{z \in S} \ell (z, \theta), \tag{3}\end{equation}

instead of the expected risk because the data distribution $\mathcal {D}$ is not known, where $\theta _t$ and η denote the output of the GD at the tth update and the learning rate, respectively. Note that the parameters $\theta _0$ are initialized according to the method specified in each subsequent analysis section.

C) Decomposition of expected risk

The expected risk is decomposed into two components,

(4)\begin{equation} R(\theta_t) = \underbrace{R(\theta_t) - R_S(\theta_t)}_{{\rm generalization\ gap}} + \underbrace{R_S(\theta_t)}_{{\rm empirical\ risk}}. \tag{4}\end{equation}

The generalization gap measures the difference between the expected risk and empirical risk, while the empirical risk expresses how fast the GD optimizes the parameters. Recent analytical techniques evaluate each of them as described next.

A) Formulation of ResNets

We evaluate the effectiveness of ResNets by deriving upper bounds of the generalization gaps of the following linear DNNs: $f: \mathbb {R}^D \times \Theta \rightarrow \mathbb {R}^D$, where θ denotes the parameters of each NN.

MLP:

(5)\begin{equation} f(x;\theta) = \prod_{l=1}^L W^l x, \tag{5}\end{equation}

ResNet1:

(6)\begin{equation} f(x; \theta) = \prod_{l=1}^L (W^l + I) x, \tag{6}\end{equation}

ResNet2:

(7)\begin{equation} f(x; \theta) = \prod_{l=1}^{L/2} (W^{2l} W^{2L-1} + I) x. \tag{7}\end{equation}

Although these DNNs are linear with respect to the input x and have the same expressive ability, they are nonconvex with respect to the parameter θ and have different parameter representations.

B) Algorithmic stability

A training algorithm $\mathcal {A}$ receives the training set S and outputs a trained model $\mathcal {A}(S)$. The algorithmic stability measures how much the removal of one example $z(n)$ from the training set S affects the trained model $\mathcal {A}(S^n)$ in terms of the expected loss, where $S^n = S {\setminus} z(n)$ (Fig. 1).

Fig. 1. Algorithmic stability. It measures how much the removal of one example $z(n)$ from the training set S affects the trained model $A(S^n)$.

Definition 1 (Definition 4 in [Reference Bousquet and Elisseeff18])

The training algorithm $\mathcal {A}$ is pointwise hypothesis stable if there exists $\epsilon _{stab}$ such that $\forall n \in [N],$

(8)\begin{equation} \mathbb{E}_{\mathcal{A},S} \left[ \left\vert \ell (z(n), \mathcal{A}(S)) - \ell (z(n), \mathcal{A}(S^n) \right\vert \right] \leq \epsilon_{stab}, \tag{8}\end{equation}

where the expectation is taken with respect to the randomness of the algorithm $\mathcal {A}$ and the training set S.

A stable algorithm $\mathcal {A}$ with small $\epsilon _{stab}$ outputs a trained model with a small generalization gap in the framework of statistical learning theory.

Theorem 1 (Theorem 11 in [Reference Bousquet and Elisseeff18])

If the training algorithm $\mathcal {A}$ is pointwise hypothesis stable, the following holds with probability at least $1-\delta :$

(9)\begin{equation} R(A(S)) - R_S(A(S)) \leq \sqrt{\frac{M^2 + 12 MN \epsilon_{stab}}{2 N \delta}}, \tag{9}\end{equation}

where M is an upper bound of the loss function.

The algorithmic stability $\epsilon _{stab}$ of the GD depends on the flatness of the loss landscape around a global minimum.

Definition 2 (Definition 4 in [Reference Charles and Papailiopoulos19])

The empirical risk $R_S$ satisfies the Polyak–Lojasiewicz $($PL$)$ condition with a constant μ if the following holds:

(10)\begin{equation} \forall \theta_t \in \Theta, \quad R_S(\theta_t) - R_S(\theta^*) \leq \frac{1}{2 \mu} \left\Vert \nabla_{\theta} R_S(\theta_t) \right\Vert ^2, \tag{10}\end{equation}

where $\theta _*$ is a global minimum.

Here, the constant μ for the PL condition expresses the flatness of the loss landscape around a global minimum. If the empirical risk $R_S$ satisfies the PL condition with μ and is β-smooth, that is, the gradient is β-Lipschitz, then the following inequality holds:

(11)\begin{equation} R_S(\theta_t) - R_S(\theta_*) \leq \frac{\beta^2}{2\mu} \left\Vert \theta_t - \prod_{\Theta_*} (\theta_t) \right\Vert ^2, \tag{11}\end{equation}

where $\prod _{\Theta _*}(\theta _t)$ is the projection of $\theta _t$ on the set of global minima $\Theta _*$. This shows that an excess risk is smaller than a quadratic function of parameters and that the constant μ controls its flatness (Fig. 2).

Fig. 2. Excess risk is smaller than a quadratic function of parameters θ. The constant μ for the PL condition controls its flatness.

A training algorithm has better stability if it converges faster and its loss function has flatter minima.

Theorem 2 (Theorem 3 in [Reference Charles and Papailiopoulos19])

Suppose that the empirical risk $R_S$ satisfies the PL condition with μ and the loss function is α-Lipschitz. If the training algorithm $\mathcal {A}$ converges parameters to the global minima $\theta _*$ with $\left\Vert \theta _t - \theta _* \right\Vert \leq \epsilon _{t},$ it is pointwise hypothesis stable, that is,

(12)\begin{equation} \epsilon_{stab} \leq 2 \alpha \epsilon_{t} + \frac{2 \alpha^2}{\mu (N-1)}. \tag{12}\end{equation}

C) Upper bounds of the generalization errors

We applied the above stability analysis to MLP, ResNet1, and ResNet2 under Assumptions 1 and 2 [Reference Furusho, Liu and Ikeda24].

These assumptions are rather weak since a dataset satisfies Assumption 1 if it is preprocessed by principal component analysis (PCA) whitening. In addition, the PCA-whitened MNIST dataset satisfies Assumption 2.

Theorem 3 (Theorems 3 and 4 in [Reference Furusho, Liu and Ikeda24])

Initialize the linear DNNs by orthogonal initialization [Reference Saxe, McClelland and Ganguli27]. Then, under Assumptions 1 and 2, ResNet2 has flatter minima, as shown in Table 1, where $a_{\min }$ and $a_{\max }$ are the minimum and maximum singular values of the weights, during training, respectively, and C is a constant $($Fig. 3$)$. In addition, its parameters converge slower than the other DNNs, as shown in Table 2, where γ is the minimum singular value of the transform by the layer during training.

Table 1. PL condition of the L-layer linear DNNs.

Fig. 3. Constants μ for the PL-condition of the 10-layer linear DNNs.

Table 2. Convergence of the L-layer linear DNNs.

D) Numerical experiments

To confirm the validity of the above analyses and the applicability to the DNNs with the ReLU activation function, $\phi (\cdot ) = \max \{ 0, \cdot \}$, some numerical experiments were carried out. The dataset was the MNIST dataset [Reference LeCun, Bottou, Bengio and Haffner28] after PCA whitening, so that $\forall d \in [D], \mathbb {E} [ x_d ] = 0$ and ${\rm Var} (x_d)=1$, and projection into the principal subspace of 10 dimensions.

We initialized the DNNs with 10 hidden units in each layer by the orthogonal initialization and trained these by the GD. During the training, the training loss, the test loss, and the approximate value of the stability $\epsilon _{stab}$,

(13)\begin{equation} \frac{1}{N} \sum_{n=1}^N \vert \ell (z(n), A(S) ) - \ell ( z(n), A(S^n) ) \vert, \tag{13}\end{equation}

were calculated every five updates (Figs. 4 and 5). The results show that ResNet2 had a greater stability and a smaller generalization gap than the other DNNs and that the analyses are valid even for the DNNs with the ReLU activation.

Fig. 4. Training loss (solid lines) and test loss (dotted lines) of the 10-layer DNNs with the ReLU activation.

Fig. 5. Approximation of the stability $\epsilon _{stab}$ of the 10-layer DNNs with the ReLU activation.

A) Formulation of batch normalization

We evaluate the effectiveness of BN by deriving the empirical risk of the following DNNs with the ReLU activation $\phi (\cdot ) = \max \{ 0, \cdot \}$:

ResNet:

(14)\begin{equation} \begin{aligned} & u^l_i = \phi (h^{l-1}_i),\quad h^l_i = \sum_{j=1}^D W^l_{i,j} u^l_j + h^{l-1}_i, \\ & \hat{y} = 1^T h^L. \end{aligned} \tag{14}\end{equation}

ResNet with BN:

(15)\begin{align} & u^l_i = \phi ( {\rm BN} (h^{l-1}_i) ),\quad {\rm BN} (h^l_i) = \frac{h^l_i - \mathbb{E}_{x} \left[ h^l_i \right]}{\sqrt{{\rm Var} \left( h^l_i \right)}}, \notag\\ & h^l_i = \sum_{j=1}^D W^l_{i,j} u^l_j + h^{l-1}_i,\quad \hat{y} = 1^T h^L. \tag{15}\end{align}

Here, $h^0 = W^0 x$ is the projection of the input, and the expectation of the BN is taken with respect to the input in the batch of the GD. Without the loss of generality, the projection matrix is a square matrix initialized by Xavier initialization [Reference Glorot and Bengio29], and the inputs in the training set are normalized to $\mathbb {E}_{x} [ x_i ] = 0$ and ${\rm Var} ( x_i ) = 1$.

B) Hessian and Fisher information matrices

The empirical risk $R_S(\theta _t)$ is approximated by the second-order Taylor expansion around the minima $\theta _*$,

(16)\begin{equation} R_S(\theta_t) = R_S(\theta_*) + \frac{1}{2} (\theta_t - \theta_*)^T H(\theta_*) (\theta_t - \theta_*), \tag{16}\end{equation}

where $H(\theta _*) = \nabla _{\theta } \nabla _{\theta } R_S(\theta _*)$ is the Hessian matrix. The Hessian matrix is decomposed as $H(\theta _*) = U \Lambda U^T$, where U and Λ are a unitary square matrix comprising the eigenvectors and a diagonal matrix filled with the eigenvalues, respectively, which simplifies the empirical risk to

(17)\begin{equation} R_S(\theta_t) = R_S(\theta_*) + \frac{1}{2} v_t^T \Lambda v_t, \tag{17}\end{equation}

where $v_t = U^T (\theta _t - \theta _*)$, and the GD to

(18)\begin{equation} v_{t+1} = \left( I - \eta \Lambda \right) v_t. \tag{18}\end{equation}

Let $\lambda _{\rm min}$ and $\lambda _{\rm max}$ be the minimum and maximum eigenvalues of $H(\theta _*)$, respectively. The GD converges when the learning rate is $\eta =2/\lambda _{\rm max}$ and it converges fastest when the learning rate is $\eta = 1/\lambda _{\rm max}$. The fastest convergence rate is the reciprocal of the condition number, $\lambda _{\rm max} / \lambda _{\rm min}$ [Reference LeCun, Bottou, Orr and Muller21], as is well known in adaptive filtering theory [Reference Haykin30]. However, the Hessian matrix and its eigenvalues are difficult to calculate owing to the complicated structure of DNNs.

Recently, the Fisher information matrix (FIM) of $p(x,y;\theta )$,

(19)\begin{equation} F(\theta) = \mathbb{E}_{x} \left[ \nabla_{\theta} \log p(x;\theta) \, \nabla_{\theta} \log p(x;\theta)^T \right], \tag{19}\end{equation}

has been found to approximate the Hessian matrix of a DNN, $f(x;\theta )$, where $p(x,y;\theta ) = p(x) p(y \vert x;\theta )$, $p(y \vert x; \theta ) = \mathcal {N} ( f(x;\theta ), 1 )$, and $p(x)$ is the probability of the input. In this case, the FIM is rewritten as

(20)\begin{equation} F(\theta) = \mathbb{E}_{x} \left[ \nabla_{\theta} f(x;\theta) \, \nabla_{\theta} f(x;\theta)^T \right] \tag{20}\end{equation}

and the following holds:

(21)\begin{equation} H(\theta) = F(\theta) - \mathbb{E}_{x} \left[ (y-f(x;\theta)) \cdot \nabla_{\theta} \nabla_{\theta} f(x;\theta) \right], \tag{21}\end{equation}

the second term of which is negligible when the error is small. In addition, the eigenvalues of the FIM of a sufficiently wide neural network do not change during training [Reference Jacot, Gabriel and Hongler31,Reference Karakida, Akaho and Amari32].

C) Bounds of the empirical risk

We calculated the eigenvalues of the FIMs of the naive ResNet and the ResNet with BN averaged over the random He initialization [Reference He, Zhang, Ren and Sun33] (expected FIM) under Assumptions 3 and 4 [Reference Furusho and Ikeda25].

Although Assumption 3 is rather unrealistic, some theorems have been derived on the basis of Assumption 3 and their results were in agreement with those of numerical experiments [Reference Karakida, Akaho and Amari22,Reference Poole, Lahiri, Sohl-Dickstein and Ganguli34,Reference Yang and Schoenholz35]. In addition, the binary class PCA-whitened MNIST dataset satisfies Assumption 4 (Fig. 6).

Fig. 6. Activation rate of the hidden units in each layer.

Theorem 4 (Modification of Table 1 in [Reference Furusho and Ikeda25])

Under Assumptions 3 and 4, the maximum eigenvalue $\lambda _{\rm max}$ of the expected FIM of the ResNet grows exponentially with the depth,

(22)\begin{equation} \begin{aligned} & m_{\lambda} \leq \lambda_{\rm max} \leq (L+1) D^2 m_{\lambda}, \\ & m_{\lambda} = \frac{L+4}{4L + 4} \cdot 2^L, \end{aligned} \tag{22}\end{equation}

where $m_{\lambda }$ is the mean of all the eigenvalues $\{ \lambda _i \}_{i=1}^{(L+1) D^2}$.

Theorem 5 (Modification of Table 1 in [Reference Furusho and Ikeda25])

Under Assumptions 3 and 4, BN relaxes the exponential growth of the eigenvalue to $L \log L$ order at most,

(23)\begin{align} & m_{\lambda} \leq \lambda_{\rm max} \leq (L+1) D^2 m_{\lambda}, \tag{23}\end{align}

(24)\begin{align} & m_{\lambda} = \frac{ H_{L+1} + 1}{2}, \tag{24}\end{align}

where $H_L = \sum _{k=1}^{L} \frac {1}{k}$ is the harmonic number.

Note that our discussion is focused on the minima $\theta _*$. This is justified by the fact that the GD and the stochastic GD make the parameters into the minima under some conditions [Reference Jacot, Gabriel and Hongler31,Reference Ge, Huang, Jin and Yuan36,Reference Lee, Jason, Scimochowitz, Jordan and Recht37].

D) Numerical experiments

To confirm the validity of the above analyses, some numerical experiments were carried out. The dataset was a subset of the MNIST dataset [Reference LeCun, Bottou, Bengio and Haffner28] with class labels of 0 and 1 after the PCA whitening so that $\forall d \in [D], \mathbb {E} [x_d ] = 0$, and ${\rm Var} (x_d) = 1$, with projection into the principal subspace of 50 dimensions.

We initialized the ResNet and ResNet with BN, which have 50 hidden units in each layer, by the He initialization, calculated the mean eigenvalues and maximum eigenvalues of the expected FIMs of these DNNs, and found that the mean eigenvalues were in agreement with the theoretical values and that the maximum eigenvalues were bounded by the theoretical upper and lower bounds (Fig. 7).

Fig. 7. Mean eigenvalues and maximum eigenvalues of the expected FIM.

In addition, the convergence properties of the ResNet and ResNet with BN were numerically examined. Each algorithm with various numbers of layers L and learning rates η updated the parameters 50 times for each run and the training loss was averaged over five runs (Fig. 8). It was found that the algorithms converged if the learning rate was less than the lower bounds of the stable convergence, that is, $2 / ( \mbox {upper bound of } \lambda _{\rm max})$.

Fig. 8. Training loss under various settings. Read lines: theoretical lower-bounds of the maximum learning rates. White color: divergence (the loss ${>}1000$). (a) ResNet. (b) ResNet with BN.

From a more practical viewpoint, we evaluated the training loss and test loss of the ResNet and ResNet with BN at each update, where each algorithm used an optimal learning rate, $\eta = 1 / ( \mbox {upper bound of } \lambda _{\rm max})$ (Fig. 9). The result shows that the BN accelerates the convergence and increases the stability.

Fig. 9. Loss of the 4-layer DNNs (solid line: average, shadowed area: within one s.d.)