Universal Statistics of Fisher Information in Deep Neural Networks: Mean Field Approach

Ryo Karakida, Shotaro Akaho, Shun-ichi Amari

Introduction

Deep learning has succeeded in making hierarchical neural networks perform excellently in various practical applications . To proceed further, it would be beneficial to give more theoretical elucidation as to why and how deep neural networks (DNNs) work well in practice. In particular, it would be useful to not only clarify the individual models and phenomena but also explore various unified theoretical frameworks that could be applied to a wide class of deep networks. One widely used approach for this purpose is to consider deep networks with random connectivity and a large width limit . For instance, Poole et al. proposed a useful indicator to explain the expressivity of DNNs. Regarding the trainability of DNNs, Schoenholz et al. extended this theory to backpropagation and found that the vanishing and explosive gradients obey a universal law. These studies are powerful in the sense that they do not depend on particular model architectures, such as the number of layers or activation functions.

Unfortunately, such universal frameworks have not yet been established in many other topics. One is the geometric structure of the parameter space. For instance, the loss landscape without spurious local minima is important for easier optimization and theoretically guaranteed in single-layer models , shallow piecewise linear ones , and extremely wide deep networks with the number of training samples smaller than the width . Flat global minima have been reported to be related to generalization ability through empirical experiments showing that networks with such minima give better generalization performance . However, theoretical analysis of the flat landscape has been limited in shallow rectified linear unit (ReLU) networks . Thus, a residual subject of interest is to theoretically reveal the geometric structure of the parameter space truly common among various deep networks.

To establish the foundation of the universal perspective of the parameter space, this study analytically investigates the Fisher information matrix (FIM). As is overviewed in Section 2.1, the FIM plays an essential role in the geometry of the parameter space and is a fundamental quantity in both statistics and machine learning.

This study analyzes the FIM of deep networks with random weights and biases, which are widely used settings to analyze the phenomena of DNNs . First, we analytically obtain novel statistics of the FIM, namely, the mean (Theorem 1), variance (Theorem 1), and maximum of eigenvalues (Theorem 4). These are universal among a wide class of shallow and deep networks with various activation functions. These quantities can be obtained from simple iterative computations of macroscopic variables. To our surprise, the mean of the eigenvalues asymptotically decreases with an order of O(1/M)O(1/M) in the limit of a large network width MM, while the variance takes a value of O(1)O(1), and the maximum eigenvalue takes a huge value of O(M)O(M) by using the O()O(\cdot) order notation. Since the eigenvalues are non-negative, these results mean that most of the eigenvalues are close to zero, but the edge of the eigenvalue distribution takes a huge value. Because the FIM defines the Riemannian metric of the parameter space, the derived statistics imply that the space is locally flat in most dimensions, but strongly distorted in others. In addition, because the FIM also determines the local shape of a loss landscape, the landscape is also expected to be locally flat while strongly distorted.

Furthermore, to confirm the potential usage of the derived statistics, we show some exercises. One is on the Fisher-Rao norm (Theorem 5). This norm was originally proposed to connect the flatness of a parameter space to the capacity measure of generalization ability. We evaluate the Fisher-Rao norm by using an indicator of the small eigenvalues, κ1\kappa_{1} in Theorem 1. Another exercise is related to the more practical issue of determining the size of the learning rate necessary for the steepest descent gradient to converge. We demonstrate that an indicator of the huge eigenvalue, κ2\kappa_{2} in Theorem 4, enables us to roughly estimate learning rates that make the gradient method converge to global minima (Theorem 7). We expect that it will help to alleviate the dependence of learning rates on heuristic settings.

2 Related works

Despite its importance in statistics and machine learning, study on the FIM for neural networks has been limited so far. This is because layer-by-layer nonlinear maps and huge parameter dimensions make it difficult to take analysis any further. Degeneracy of the eigenvalues of the FIM has been found in certain parameter regions . To understand the loss landscape, Pennington and Bahri has utilized random matrix theory and obtained the spectrum of FIM and Hessian under several assumptions, although the analysis is limited to special types of shallow networks. In contrast, this paper is the first attempt to apply the mean field approach, which overcomes the difficulties above and enables us to identify universal properties of the FIM in various types of DNNs.

LeCun et al. investigated the Hessian of the loss, which coincides with the FIM at zero training error, and empirically reported that very large eigenvalues exist, i.e., ”big killers”, which affects the optimization (discussed in Section 4.2). The eigenvalue distribution peaks around zero while its tail is very long; this behavior has been empirically known for decades , but its theoretical evidence and evaluation have remained unsolved as far as we know. Therefore, our theory provides novel theoretical evidence that this skewed eigenvalue distribution and its huge maximum appear universally in DNNs.

The theoretical tool we use here is known as the mean field theory of deep networks as briefly overviewed in Section 2.4. This method has been successful in analyzing neural networks with random weights under a large width limit and in explaining the performance of the models. In particular, it quantitatively coincides with experimental results very well and can predict appropriate initial values of parameters for avoiding the vanishing or explosive gradient problems . This analysis has been extended from fully connected deep networks to residual and convolutional networks . The evaluation of the FIM in this study is also expected to be extended to such cases.

Preliminaries

We focus on the Fisher information matrix (FIM) of neural network models, which previous works have developed and is commonly used . It is defined by

This study investigates the above empirical FIM for arbitrary TT. It converges to the expected FIM as TT\rightarrow\infty. Although the form of the FIM changes a bit in other statistical models (i.g., softmax outputs), these differences are basically limited to the multiplication of activations in the output layer . Our framework can be straightforwardly applied to such cases.

The FIM determines the asymptotic accuracy of the estimated parameters, as is known from a fundamental theorem of statistics, namely, the Cramér-Rao bound . Below, we summarize a more intuitive understanding of the FIM from geometric views.

Information geometric view. Let us define an infinitesimal squared distance dr2dr^{2}, which represents the Kullback-Leibler divergence between the statistical model p(x,y;θ)p(x,y;\theta) and p(x,y;θ+dθ)p(x,y;\theta+d\theta) against a perturbation dθd\theta. It is given by

Loss landscape view. The empirical FIM (2) determines the local landscape of the loss function around the global minimum. Suppose we have a squared loss function E(θ)=(1/2T)tTy(t)fθ(t)2E(\theta)=(1/2T)\sum_{t}^{T}||y(t)-f_{\theta}(t)||^{2}. The FIM is related to the Hessian of the loss function, H:=θθE(θ)H:=\nabla_{\theta}\nabla_{\theta}E(\theta), in the following way:

The Hessian coincides with the FIM when the parameter converges to the global minimum by learning, that is, the true parameter θ\theta^{*} from which the teacher signal y(t)y(t) is generated by y(t)=fθ(t)y(t)=f_{\theta^{*}}(t) or, more generally, with noise (i.e., y(t)=fθ(t)+εty(t)=f_{\theta^{*}}(t)+\varepsilon_{t}, where εt\varepsilon_{t} denotes zero-mean Gaussian noise) . In the literature on deep learning, its eigenvectors whose eigenvalues are close to zero locally compose flat minima, which leads to better generalization empirically . Modifying the loss function with the FIM has also succeeded in overcoming the catastrophic forgetting .

Note that the information geometric view tells us more than the loss landscape. While the Hessian (4) assumes the special teacher signal, the FIM works as the Riemannian metric to arbitrary teacher signals.

2 Network architecture

This study investigates a fully connected feedforward neural network. The network consists of one input layer with M0M_{0} units, L1L-1 hidden layers (L2L\geq 2) with MlM_{l} units per hidden layer (l=1,2,...,L1)(l=1,2,...,L-1), and one output layer with MLM_{L} units:

This study focuses on the case of linear outputs, that is, fθ,k(x)=hkL=ukLf_{\theta,k}(x)=h^{L}_{k}=u_{k}^{L}. We assume that the activation function ϕ(x)\phi(x) and its derivative ϕ(x):=dϕ(x)/dx\phi^{\prime}(x):=d\phi(x)/dx are square-integrable functions on a Gaussian measure. A wide class of activation functions, including the sigmoid-like and (leaky-) ReLU functions, satisfy these conditions. Different layers may have different activation functions. Regarding the network width, we set Ml=αlM  (lL1)M_{l}=\alpha_{l}M\ \ (l\leq L-1) and consider the limiting case of large MM with constant coefficients αl\alpha_{l}. This study mainly focuses on the case where the number of output units is given by a constant ML=CM_{L}=C. The higher-dimensional case of C=O(M)C=O(M) is argued in Section 4.3.

The FIM (2) of a deep network is computed by the chain rule in a manner similar to that of the backpropagation algorithm:

where δk,il:=fθ,k/uil\delta_{k,i}^{l}:=\partial f_{\theta,k}/\partial u_{i}^{l} for (k=1,...,Ck=1,...,C). To avoid the complicated notation, we omit the index of the output unit, i.e., δil=δk,il\delta_{i}^{l}=\delta_{k,i}^{l}, in the following.

3 Random connectivity

The parameter set θ={Wijl,bil}\theta=\{W_{ij}^{l},b_{i}^{l}\} is an ensemble generated by

and then fixed, where N(0,σ2)\mathcal{N}(0,\sigma^{2}) denotes a Gaussian distribution with zero mean and variance σ2\sigma^{2}, and we set σwl>0\sigma_{w^{l}}>0 and σbl>0\sigma_{b^{l}}>0. To avoid complicated notation, we set them uniformly as σwl2=σw2\sigma_{w^{l}}^{2}=\sigma_{w}^{2} and σbl2=σb2\sigma_{b^{l}}^{2}=\sigma_{b}^{2}, but they can easily be generalized. It is essential to normalize the variance of the weights by MM in order to normalize the output uilu_{i}^{l} to O(1)O(1). This setting is similar to how parameters are initialized in practice . We also assume that the input samples hi0(t)=xi(t)  (t=1,...,T)h_{i}^{0}(t)=x_{i}(t)\ \ (t=1,...,T) are generated in an i.i.d. manner from a standard Gaussian distribution: xi(t)N(0,1).x_{i}(t)\sim\mathcal{N}(0,1). We focus here on the Gaussian case for simplicity, although we can easily generalize it to other distributions with finite variances.

Let us remark that the above random connectivity is a common setting widely supposed in theories. Analyzing such a network can be regarded as the typical evaluation . It is also equal to analyzing the network randomly initialized . The random connectivity is often assumed in the analysis of optimization as a true parameter of the networks, that is, the global minimum of the parameters .

4 Mean-field approach

On neural networks with random connectivity, taking a large width limit, we can analyze the asymptotic behaviors of the networks. Recently, this asymptotic analysis is referred to as the mean field theory of deep networks, and we follow the previously reported notations and terminology .

First, let us introduce the following variables for feedforward signal propagations: q^l:=ihil(t)2/Ml\hat{q}^{l}:=\sum_{i}h^{l}_{i}(t)^{2}/M_{l} and q^stl:=ihil(s)hil(t)/Ml\hat{q}^{l}_{st}:=\sum_{i}h_{i}^{l}(s)h_{i}^{l}(t)/M_{l}. In the context of deep learning, these variables have been utilized to explain the depth to which signals can sufficiently propagate. The variable q^stl\hat{q}^{l}_{st} is the correlation between the activations for different input samples x(s)x(s) and x(t)x(t) in the ll-th layer. Under the large MM limit, these variables are given by integration over Gaussian distributions because the pre-activation uliu_{l}^{i} is a weighted sum of independent random parameters and the central limit theorem is applicable :

with q^0=1\hat{q}^{0}=1 and q^st0=0\hat{q}_{st}^{0}=0 (l=0,...,L1l=0,...,L-1). We can generalize the theory to unnormalized data with q^00\hat{q}^{0}\neq 0 and q^st00\hat{q}^{0}_{st}\neq 0, just by substituting them into the recurrence relations. The notation Du=duexp(u2/2)/2πDu=du\exp(-u^{2}/2)/\sqrt{2\pi} means integration over the standard Gaussian density. Here, the notation I[,]I_{\cdot}[\cdot,\cdot] represents the following integral: Iϕ[a,b]=Dz1Dz2ϕ(az1)ϕ(a(cz1+1c2z2))I_{\phi}[a,b]=\int Dz_{1}Dz_{2}\phi\left(\sqrt{a}z_{1}\right)\phi\left(\sqrt{a}(cz_{1}+\sqrt{1-c^{2}}z_{2})\right) with c=b/ac=b/a. The qstlq_{st}^{l} is linked to the compositional kernel and utilized as the kernel of the Gaussian process .

This assumption makes the dependence between ϕ(uil)\phi(u_{i}^{l}) (or ϕ(uil)\phi^{\prime}(u_{i}^{l})) and δjl+1\delta_{j}^{l+1}, which share the same parameter set, very weak, and one can regard it as independent. It enables us to apply the central limit theorem to the backpropagated chain (7). Thus, the previous studies derived the following recurrence relations (l=0,...,L1l=0,...,L-1):

Fundamental FIM statistics

Here, we report mathematical findings that the mean, variance, and maximum of eigenvalues of the FIM (2) are explicitly expressed by using macroscopic variables. Our theorems are universal for networks ranging in size from shallow (L=2L=2) to arbitrarily deep (L3L\geq 3) with various activation functions.

The FIM is a P×PP\times P matrix, where PP represents the total number of parameters. First, we compute the arithmetic mean of the FIM’s eigenvalues as mλ:=i=1Pλi/Pm_{\lambda}:=\sum_{i=1}^{P}\lambda_{i}/P. We find a hidden relation between the macroscopic variables and the statistics of FIM:

Suppose that Assumption 1 holds. In the limit of M1M\gg 1, the mean of the FIM’s eigenvalues is given by

Because the FIM is a positive semi-definite matrix and its eigenvalues are non-negative, this theorem means that most of the eigenvalues asymptotically approach zero when MM is large. Recall that the FIM determines the local geometry of the parameter space. The theorem suggests that the network output remains almost unchanged against a perturbation of the parameters in many dimensions. It also suggests that the shape of the loss landscape is locally flat in most dimensions.

Furthermore, by using Markov’s inequality, we can prove that the number of larger eigenvalues is limited, as follows:

Let us denote the number of eigenvalues satisfying λk\lambda\geq k by N(λk)N(\lambda\geq k) and suppose that Assumption 1 holds. For a constant k>0k>0, N(λk)min{ακ1CM/k,CT}N(\lambda\geq k)\leq\min\{\alpha\kappa_{1}CM/k,CT\} holds in the limit of M1M\gg 1.

The proof is shown in Supplementary Material A.2. When TT is sufficiently small, we have a trivial upper bound N(λk)CTN(\lambda\geq k)\leq CT and the number of non-zero eigenvalue is limited. The corollary clarifies that even when TT becomes large, the number of eigenvalues whose values are O(1)O(1) is O(M)O(M) at most, and still much smaller than the total number of parameters PP.

2 Variance of eigenvalues

Next, let us consider the second moment sλ:=i=1Pλi2/Ps_{\lambda}:=\sum_{i=1}^{P}\lambda_{i}^{2}/P. We now demonstrate that sλs_{\lambda} can be computed from the macroscopic variables:

Suppose that Assumption 1 holds. In the limit of M1M\gg 1, the second moment of the FIM’s eigenvalues is

The proof is shown in Supplementary Material A.3.

From Theorems 1 and 3, we can conclude that the variance of the eigenvalue distribution, sλmλ2s_{\lambda}-m_{\lambda}^{2}, is O(1)O(1). Because the mean mλm_{\lambda} is O(1/M)O(1/M) and most eigenvalues are close to zero, this result means that the edge of the eigenvalue distribution takes a huge value.

3 Maximum eigenvalue

As we have seen so far, the mean of the eigenvalues is O(1/M)O(1/M), and the variance is O(1)O(1). Therefore, we can expect that at least one of the eigenvalues must be huge. Actually, we can show that the maximum eigenvalue (that is, the spectral norm of the FIM) increases in the order of O(M)O(M) as follows.

Suppose that Assumption 1 holds. In the limit of M1M\gg 1, the maximum eigenvalue of the FIM is

The λmax\lambda_{max} is derived from the dual matrix FF^{*} (detailed in Supplemental Material A.4). If we take the limit TT\rightarrow\infty, we can characterize the quantity κ2\kappa_{2} by the maximum eigenvalue as λmax=ακ2M\lambda_{max}=\alpha\kappa_{2}M. Note that λmax\lambda_{max} is independent of CC. When C=O(M)C=O(M), it may depend on CC, as shown in Section 3.4.

This theorem suggests that the network output changes dramatically with a perturbation of the parameters in certain dimensions and that the local shape of the loss landscape is strongly distorted in that direction. Here, note that λmax\lambda_{max} is proportional to α\alpha, which is the summation over LL terms. This means that, when the network becomes deeper, the parameter space is more strongly distorted.

We confirmed the agreement between our theory and numerical experiments, as shown in Fig. 1. Three types of deep networks with parameters generated by random connectivity (8) were investigated: tanh, ReLU, and linear activations (L=3L=3, αl=C=1\alpha_{l}=C=1). The input samples were generated using i.i.d. Gaussian samples, and T=102T=10^{2}. When P>TP>T, we calculated the eigenvalues by using the dual matrix FF^{*} (defined in Supplementary Material A.3) because FF^{*} is much smaller and its eigenvalues are easy to compute. The theoretical values of mλm_{\lambda}, sλs_{\lambda} and λmax\lambda_{max} agreed very well with the experimental values in the large MM limit. We could predict mλm_{\lambda} even for small MM. In addition, In Supplementary Material C.1, we also show the results of experiments with fixed MM and changing TT. The theoretical values coincided with the experimental values very well for any TT as the theorems predict.

Connections to learning strategies

Here, we show some applications that demonstrate how our universal theory on the FIM can potentially enrich deep learning theories. It enables us to quantitatively measure the behaviors of learning strategies as follows.

Recently, Liang et al. proposed the Fisher-Rao norm for a capacity measure of generalization ability:

where θ\theta represents weight parameters. They reported that this norm has several desirable properties to explain the high generalization capability of DNNs. In deep linear networks, its generalization capacity (Rademacher complexity) is upper bounded by the norm. In deep ReLU networks, the Fisher-Rao norm serves as a lower bound of the capacities induced by other norms, such as the path norm and the spectral norm . The Fisher-Rao norm is also motivated by information geometry, and invariant under node-wise linear rescaling in ReLU networks. This is a desirable property to connect capacity measures with flatness induced by the rescaling .

Here, to obtain a typical evaluation of the norm, we define the average over possible parameters with fixed variances (σw2,σb2\sigma_{w}^{2},\sigma_{b}^{2}) by θ=iDθi()\langle\cdot\rangle_{\theta}=\int\prod_{i}D\theta_{i}(\cdot), which leads to the following theorem:

Suppose that Assumption 1 holds. In the limit of M1M\gg 1, the Fisher-Rao norm of DNNs satisfies

where αmin=miniαi\alpha_{min}=\min_{i}\alpha_{i}. Equality holds in a network with a uniform width Ml=MM_{l}=M, and then we have θFRθ=σw2(L1)Cκ1\langle||\theta||_{FR}\rangle_{\theta}=\sigma_{w}^{2}(L-1)C\kappa_{1}.

The proof is shown in Supplementary Material A.6. Although what we can evaluate is only the average of the norm, it can be quantified by κ1\kappa_{1}. This guarantees that the norm is independent of the network width in the limit of M1M\gg 1, which was empirically conjectured by .

Recently, Smith and Le argued that the Bayesian factor composed of the Hessian of the loss function, whose special case is the FIM, is related to the generalization. Similar analysis to the above theorem may enable us to quantitatively understand the relation between the statistics of the FIM and the indicators to measure the generalization ability.

2 Learning rate for convergence

Consider the steepest gradient descent method in a batch regime. Its update rule is given by

where η\eta is a constant learning rate. We have added a momentum term with a coefficient μ\mu because it is widely used in training deep networks. Assume that the squared loss function E(θ)E(\theta) of Eq. (4) has a global minimum θ\theta^{*} achieving the zero training error E(θ)=0E(\theta^{*})=0. Then, the FIM’s maximum eigenvalue is dominant over the convergence of learning as follows:

A learning rate satisfying η<2(1+μ)/λmax\eta<2(1+\mu)/\lambda_{max} is necessary for the steepest gradient method to converge to the global minimum θ\theta^{*}.

The proof is given by the expansion around the minimum, i.e., E(θ+dθ)=dθTFdθE(\theta^{*}+d\theta)=d\theta^{T}Fd\theta (detailed in Supplementary Material A.7). This lemma is a generalization of LeCun et al. , which proved the case of μ=0\mu=0. Let us refer to ηc:=2(1+μ)/λmax\eta_{c}:=2(1+\mu)/\lambda_{max} as the critical learning rate. When η>ηc\eta>\eta_{c}, the gradient method never converges to the global minimum. The previous work also claimed that η=ηc/2\eta=\eta_{c}/2 is the best choice for fastest convergence around the minimum. Although we focus on the batch regime, the eigenvalues also determine the bound of the gradient norms and the convergence of learning in the online regime .

Then, combining Lemma 6 with Theorem 4 leads to the following:

Suppose that Assumption 1 holds. Let a global minimum θ\theta^{*} be generated by Eq. (8) and satisfying E(θ)=0E(\theta^{*})=0. In the limit of M1M\gg 1, the gradient method never converges to θ\theta^{*} when

Theorem 7 quantitatively reveals that, the wider the network becomes, the smaller the learning rate we need to set. In addition, α\alpha is the sum over LL constant positive terms, so a deeper network requires a finer setting of the learning rate and it will make the optimization more difficult. In contrast, the expressive power of the network grows exponentially as the number of layers increases . We thus expect there to be a trade-off between trainability and expressive power.

To confirm the effectiveness of Theorem 7, we performed several experiments. As shown in Fig. 2, we exhaustively searched training losses while changing MM and η\eta, and found that the theoretical estimation coincides well with the experimental results. We trained deep networks (L=4L=4, αl=1\alpha_{l}=1, C=10C=10) and the loss function was given by the squared error.

The left column of Fig. 2 shows the results of training on artificial data. We generated training samples x(t)x(t) in the Gaussian manner (T=100T=100) and teacher signals y(t)y(t) by the teacher network with a true parameter set θ\theta^{*} satisfying Eq. (8). We used the gradient method (19) with μ=0.9\mu=0.9 and trained the DNNs for 100100 steps. The variances (σw2,σb2)(\sigma_{w}^{2},\sigma_{b}^{2}) of the initialization of the parameters were set to the same as the global minimum. We found that the losses of the experiments were clearly divided into two areas: one where the gradient exploded (gray area) and the other where it was converging (colored area). The red line is ηc\eta_{c} theoretically calculated using κ1\kappa_{1} and κ2\kappa_{2} on (σw2,σb2)(\sigma_{w}^{2},\sigma_{b}^{2}) of the initial parameters. Training on the regions above ηc\eta_{c} exploded, just as Theorem 7 predicts. The explosive region with η<ηc\eta<\eta_{c} got smaller in the limit of large MM.

We performed similar experiments on benchmark datasets and found that the theory can estimate the appropriate learning rates. The results on MNIST are shown in the right column of Fig. 2. As shown in Supplementary Material C.2, the results of training on CIFAR-10 were almost the same as those of MINIST. We used stochastic gradient descent (SGD) with a mini-batch size of 500500 and μ=0.9\mu=0.9, and trained the DNNs for 11 epoch. Each training sample was x(t)x(t) normalized to zero mean and variance 11 (T=50000T=50000). The initial values of (σw2,σb2)(\sigma_{w}^{2},\sigma_{b}^{2}) were set to the vicinity of the special parameter region, i.e., the critical line of the order-to-chaos transition, which the previous works recommended to use for achieving high expressive power and trainability. Note that the variances (σw2,σb2)(\sigma_{w}^{2},\sigma_{b}^{2}) may change from the initialization to the global minimum, and the conditions of the global minimum in Theorem 7 do not hold in general. Nevertheless, the learning rates estimated by Theorem 7 explained the experiments well. Therefore, the ideal conditions supposed in Theorem 7 seem to hold effectively. This may be explained by the conjecture that the change from the initialization to the global minima is small in the large limit .

Theoretical estimations of learning rates in deep networks have so far been limited; such gradients as AdaGrad and Adam also require heuristically determined hyper-parameters for learning rates. Extending our framework would be beneficial in guessing learning rates to prevent the gradient update from exploding.

3 Multi-label classification with high dimensionality

This study mainly focuses on the multi-dimensional output of C=O(1)C=O(1). This is because the number of labels is much smaller than the number of hidden units in most practice cases. However, since classification problems with far more labels are sometimes examined in the context of machine learning , it would be helpful to remark on the case of C=O(M)C=O(M) here. Denote the mean of the FIM’s eigenvalues in the case of C=O(M)C=O(M) as mλm_{\lambda}^{\prime} and so on. Straightforwardly, we can derive

The derivation is shown in Supplementary Material A.5. The mean of eigenvalues has the same form as Eq. (13) obtained in the case of C=O(1)C=O(1). The second moment and maximum eigenvalues can be evaluated by the form of inequalities. We found that the mean is of O(1)O(1) while the maximum eigenvalue is of O(M)O(M) at least and of O(M2)O(M^{2}) at most. Therefore, the eigenvalue distribution is more widely distributed than the case of C=O(1)C=O(1).

Conclusion and discussion

The present work elucidated the asymptotic statistics of the Fisher information matrix (FIM) common among deep networks with any number of layers and various activation functions. The statistics of FIM are characterized by the small mean of eigenvalues and the huge maximum eigenvalue, which are computed by the recurrence relations. This suggests that the parameter space determined by the FIM is locally flat in many directions while highly distorted in certain others. As examples of how one can connect the derived statistics to learning strategies, we suggest the Fisher-Rao norm and learning rates of steepest gradient descents.

We demonstrated that the experiments with the Gaussian prior on the parameters coincided well with the theory. Basically, the mean field theory is based on the central limit theorem with the parameters generated in an i.i.d. manner with finite variances. Therefore, one can expect that the good agreement with the theory is not limited to the experiments with the Gaussian prior. Further experiments will be helpful to clarify the applicable scope of the mean field approach.

The derived statistics are also of potential importance to other learning strategies, for instance, natural gradient methods. When the loss landscape is non-uniformly distorted, naive gradient methods are likely to diverge or become trapped in plateau regions, but the natural gradient, F1θE(θ)F^{-1}\nabla_{\theta}E(\theta), converges more efficiently . Because it normalizes the distortion of the loss landscape, the naive extension of Section 4.2 to the natural gradient leads to ηc=2(1+μ)\eta_{c}=2(1+\mu) and it seems to be much easier to choose the appropriately sized learning rate. However, we found that the FIM has many eigenvalues close to zero, and the inversion of it would make the gradient very unstable. In practice, several experiments showed that the choice of damping term ε\varepsilon, introduced in (F+εI)1θE(θ)(F+\varepsilon I)^{-1}\nabla_{\theta}E(\theta), is crucial to its performance in DNNs . The development of practical natural gradient methods will require modification such as damping.

It would also be interesting for our framework to quantitatively reveal the effects of normalization methods on the FIM. In particular, batch normalization may alleviate the larger eigenvalues because it empirically allows larger learning rates for convergence . It would also be fruitful to investigate the eigenvalues of the Hessian with a large error (4) and to theoretically quantify the negative eigenvalues that lead to the existence of saddle points and the loss landscapes without spurious local minima . The global structure of the parameter space should be also explored. We can hypothesize that the parameters are globally connected through the locally flat dimensions and compose manifolds of flat minima.

Our framework on FIMs is readily applicable to other architectures such as convolutional networks and residual networks by using the corresponding mean field theories . To this end, it may be helpful to remark that macroscopic variables in residual networks essentially diverge at the extreme depths . If one considers extremely deep residual networks, the statistics will require a careful examination of the order of the network width and the explosion of the macroscopic variables. We expect that further studies will establish a mathematical foundation of deep learning from the perspective of the large limit.

This work was partially supported by a Grant-in-Aid for Research Activity Start-up (17H07390) from the Japan Society for the Promotion of Science (JSPS).

References

Supplementary Materials

A.1 Theorem 1

To avoid complicating the notation, we first consider the case of the single output (C=1C=1). The general case is shown after. The network output is denoted by f(t)f(t) here. We denote the Fisher information matrix with full components as

The contribution of mλ(b)m_{\lambda}^{(b)} is negligible in the large MM limit as follows. The first term is

where αl\alpha_{l} comes from Ml=αlMM_{l}=\alpha_{l}M, and α\alpha comes from P=αM2P=\alpha M^{2}.

In contrast, the contributions of the bias entries are smaller than those of the weight entries in the limit of M1M\gg 1, as is easily confirmed:

mλ(W)m_{\lambda}^{(W)} is O(1/M)O(1/M) while mλ(b)m_{\lambda}^{(b)} is O(1/M2)O(1/M^{2}). Hence, the mean mλ(b)m_{\lambda}^{(b)} is negligible and we obtain mλ=κ1/Mm_{\lambda}=\kappa_{1}/M.

(ii) C>1𝐶1C>1 of O​(1)𝑂1O(1)

We can apply the above computation of C=1C=1 to each network output fk\nabla f_{k} (k=1,...,Ck=1,...,C):

Therefore, the mean of the eigenvalues becomes

A.2 Corollary 2

Because the FIM is a positive semi-definite matrix, its eigenvalues are non-negative. For a constant k>0k>0, we obtain

This is known as Markov’s inequality. When M1M\gg 1, combining this with Theorem 1 immediately yields

Because CTCT is also a trivial upper bound of N(λk)N(\lambda\geq k), we obtain Corollary 2. \blacksquare

A.3 Theorem 3

Here, let us express the FIM as F=θfθfT/TF=\nabla_{\theta}f\nabla_{\theta}f^{T}/T, where θf\nabla_{\theta}f is a P×TP\times T matrix whose columns are the gradients on each input sample, i.e., θf(t)\nabla_{\theta}f(t) (t=1,...,T)(t=1,...,T). We also introduce a dual matrix of FF, that is, FF^{*}:

Note that FF is a P×PP\times P matrix while FF^{*} is a T×TT\times T matrix. We can easily confirm that these FF and FF^{*} have the same non-zero eigenvalues.

Thus, under the limit of M1M\gg 1, the dual matrix is asymptotically given by

Neglecting the lower order term, we obtain

(ii) C>1𝐶1C>1 of O​(1)𝑂1O(1)

Here, we introduce the following dual matrix FF^{*}:

where θfk\nabla_{\theta}f_{k} is a P×TP\times T matrix whose columns are the gradients on each input sample, i.e., θfk(t)\nabla_{\theta}f_{k}(t) (t=1,...,T)(t=1,...,T), and BB is a P×CTP\times CT matrix. The FIM is represented by F=BBT/TF=BB^{T}/T. FF^{*} is a CT×CTCT\times CT matrix and consists of T×TT\times T block matrices,

The diagonal block F(k,k)F^{*}(k,k) is evaluated in the same way as the case of C=1C=1. It becomes αMK/T\alpha MK/T as shown in Eq. (A.23). The non-diagonal block F(k,k)F^{*}(k,k^{\prime}) has the following stst-th entries:

where δk,k\delta_{k,k^{\prime}} is the Kronecker delta.

where the first term comes from the diagonal blocks of O(M)O(M) and the second one is their lower order term. The third term comes from the non-diagonal blocks of O(M)O(\sqrt{M}). As one can see from here, when C=O(M)C=O(M), the thrid term becomes non-negligible. This case is examined in Section 4.3. \blacksquare

A.4 Theorem 4

Because FF and FF^{*} have the same non-zero eigenvalues, what we should derive here is the maximum eigenvalue of FF^{*}. As shown in Eq. (A.23), the leading term of FF^{*} asymptotically becomes αMK/T\alpha MK/T in the limit of M1M\gg 1. The eigenvalues of αMK/T\alpha MK/T are explicitly obtained as follows: λmax=α(T1Tκ2+1Tκ1)M\lambda_{max}=\alpha\left(\frac{T-1}{T}\kappa_{2}+\frac{1}{T}\kappa_{1}\right)M for an eigenvector e=(1,...,1)e=(1,...,1), and λi=α(κ1κ2)M/T\lambda_{i}=\alpha(\kappa_{1}-\kappa_{2})M/T for eigenvectors e1eie_{1}-e_{i} (i=2,...,Ti=2,...,T) where eie_{i} denotes a unit vector whose entries are 11 for the ii-th entry and otherwise. Thus, we obtain λmax=α(T1Tκ2+1Tκ1)M\lambda_{max}=\alpha\left(\frac{T-1}{T}\kappa_{2}+\frac{1}{T}\kappa_{1}\right)M.

(ii) C>1𝐶1C>1 of O​(1)𝑂1O(1)

Let us denote FF^{*} shown in Eq. (A.31) by F:=Fˉ+RF^{*}:=\bar{F}^{*}+R. Fˉ\bar{F}^{*} is the leading term of FF^{*} and given by a CT×CTCT\times CT block diagonal matrix whose diagonal blocks are given by αMK/T\alpha MK/T. RR denotes the residual term of O(M)/TO(\sqrt{M})/T. In general, the maximum eigenvalue is denoted by the spectral norm 2||\cdot||_{2}, that is, λmax=F2\lambda_{max}=||F^{*}||_{2}. Using the triangle inequality, we have

We can obtain Fˉ2=α(T1Tκ2+1Tκ1)M||\bar{F}^{*}||_{2}=\alpha\left(\frac{T-1}{T}\kappa_{2}+\frac{1}{T}\kappa_{1}\right)M because the maximum eigenvalues of the diagonal blocks are the same as the case of C=1C=1. Regarding R2||R||_{2}, this is bounded by R2RF=C2st(O(M)/T)2=O(CM)||R||_{2}\leq||R||_{F}=\sqrt{C^{2}\sum_{st}(O(\sqrt{M})/T)^{2}}=O(C\sqrt{M}). Therefore, when C=O(1)C=O(1), we can neglect R2||R||_{2} of O(M)O(\sqrt{M}) compared to Fˉ2||\bar{F}^{*}||_{2} of O(M)O(M).

On the other hand, we can also derive the lower bound of λmax\lambda_{max} as follows. In general, we have

where v1v_{1} is a CTCT-dimensional vector whose first TT entries are 1/T1/\sqrt{T} and the others are , that is, v1=(1,...,1,0,...,0)/T.v_{1}=(1,...,1,0,...,0)/\sqrt{T}. We can compute this lower bound by taking the sum over the entries of F(1,1)F^{*}(1,1), which is equal to Eq. (A.23):

Finally, we find that the upper bound (A.34) and lower bound (A.37) asymptotically take the same value of O(M)O(M), that is, λmax=(T1Tκ2+1Tκ1)M\lambda_{max}=\left(\frac{T-1}{T}\kappa_{2}+\frac{1}{T}\kappa_{1}\right)M.

A.5 Case of C=O​(M)𝐶𝑂𝑀C=O(M)

The mean of eigenvalues mλm_{\lambda}^{\prime} is derived in the same way as shown in Section A.1 (ii), that is, mλ=Cκ1/Mm^{\prime}_{\lambda}=C\kappa_{1}/M.

Regarding the second moment sλs_{\lambda}^{\prime}, the lower order term becomes non-negligible as remarked in Eq. (A.33). We evaluate this sλs_{\lambda}^{\prime} by using inequalities as follows:

As shown in Section A.3, for any kk, we obtain θfkT(s)θfk(t)F2/P=α(T1Tκ22+1Tκ12)||\nabla_{\theta}f_{k}^{T}(s)\nabla_{\theta}f_{k}(t)||_{F}^{2}/P=\alpha\left(\frac{T-1}{T}\kappa_{2}^{2}+\frac{1}{T}\kappa_{1}^{2}\right) in the limit of M1M\gg 1. Thus, the lower bound becomes the same form as sλs_{\lambda}, That is, sλ=Cα(T1Tκ22+1Tκ12)s_{\lambda}=C\alpha(\frac{T-1}{T}\kappa_{2}^{2}+\frac{1}{T}\kappa_{1}^{2}) . In contrast, the upper bound is given by

where FkF_{k} denotes the FIM of the kk-th output, i.e., Fk:=tθfk(t)θfk(t)T/TF_{k}:=\sum_{t}\nabla_{\theta}f_{k}(t)\nabla_{\theta}f_{k}(t)^{T}/T. Therefore, the upper bound is reduced to the summation over sλs_{\lambda} of C=1C=1. In the limit of M1M\gg 1, we obtain sλC2FkF2/P=C2α(T1Tκ22+1Tκ12)=Csλs_{\lambda}^{\prime}\leq C^{2}||F_{k}||_{F}^{2}/P=C^{2}\alpha\left(\frac{T-1}{T}\kappa_{2}^{2}+\frac{1}{T}\kappa_{1}^{2}\right)=Cs_{\lambda}.

Next, we show inequalities for λmax\lambda_{max}. We have already derived the lower bound (A.37) and this bound holds in the case of C=O(M)C=O(M) as well. In contrast, the upper bound (A.34) may become loose when CC is larger than O(1)O(1) because of the residual term R2||R||_{2}. Although it is hard to explicitly obtain the value of R2||R||_{2}, the following upper bound holds and is easy to compute by using sλs_{\lambda} of Eq. (14). Because the FIM is a positive semi-definite matrix, λi0\lambda_{i}\geq 0 holds by definition. Then, we have λmaxiλi2\lambda_{max}\leq\sqrt{\sum_{i}\lambda_{i}^{2}}. Combining this with sλ=iλi2/Ps_{\lambda}^{\prime}=\sum_{i}\lambda_{i}^{2}/P, we have λmaxαsλMαCsλM\lambda_{max}\leq\sqrt{\alpha s_{\lambda}^{\prime}}M\leq\sqrt{\alpha Cs_{\lambda}}M.

A.6 Theorem 5

where F(l,ij),(l,ab)F_{(l,ij),(l^{\prime},ab)} represents an entry of the FIM, that is, kCtWijlfk(t)Wablfk(t)/T\sum_{k}^{C}\sum_{t}\nabla_{W^{l}_{ij}}f_{k}(t)\nabla_{W^{l^{\prime}}_{ab}}f_{k}(t)/T. Because F(l,ij),(l,ab)F_{(l,ij),(l^{\prime},ab)} includes the random variables WijlW^{l}_{ij} and WablW^{l^{\prime}}_{ab}, we consider the following expansion. Note that WijlW^{l}_{ij} and WablW^{l^{\prime}}_{ab} are infinitesimals generated by Eq. (8)(8). Performing a Taylor expansion around Wijl=Wabl=0W^{l}_{ij}=W^{l^{\prime}}_{ab}=0, we obtain

where θ{\theta}^{*} is the parameter set {Wijl,bil}\{W^{l}_{ij},b^{l}_{i}\} with Wijl=Wabl=0W^{l}_{ij}=W^{l^{\prime}}_{ab}=0. By substituting the above expansion into the Fisher-Rao norm and taking the average θ\langle\cdot\rangle_{\theta}, we obtain the following leading term:

For, (l,ij)(l,ab)(l,ij)\neq(l^{\prime},ab), the last line becomes zero because of WijlWabl{Wijl,Wabl}=WijlWijlWablWabl=0\langle W^{l}_{ij}W^{l^{\prime}}_{ab}\rangle_{\{W^{l}_{ij},W^{l^{\prime}}_{ab}\}}=\langle W^{l}_{ij}\rangle_{W^{l}_{ij}}\langle W^{l^{\prime}}_{ab}\rangle_{W^{l^{\prime}}_{ab}}=0. For (l,ij)=(l,ab)(l,ij)=(l^{\prime},ab), we have (Wijl)2{Wijl}=σw2/Ml1\langle(W^{l}_{ij})^{2}\rangle_{\{W^{l}_{ij}\}}=\sigma^{2}_{w}/M_{l-1}. After all, in the limit of M1M\gg 1, we obtain

A.7 Lemma 6

Suppose a perturbation around the global minimum: θt=θ+Δt\theta_{t}=\theta^{*}+\Delta_{t}. Then, the gradient update becomes

where we have used E(θ)=0E(\theta^{*})=0 and E(θ)/θ=0\partial E(\theta^{*})/\partial\theta=0.

Consider a coordinate transformation from Δt\Delta_{t} to Δˉt\bar{\Delta}_{t} that diagonalizes FF. It does not change the stability of the gradients. Accordingly, we can update the ii-th component as follows:

Solving its characteristic equation, we obtain the general solution,

where AA and BB are constants. This recurrence relation converges if and only if ηλi<2(1+μ)\eta\lambda_{i}<2(1+\mu) for all ii. Therefore, η<2(1+μ)/λmax\eta<2(1+\mu)/\lambda_{max} is necessary for the steepest gradient to converge to θ\theta^{*}. \blacksquare

Appendix B Analytical recurrence relations

Consider the following error function as an activation function ϕ(x)\phi(x):

Hence, the recurrence relations for the feedforward signals (9) have the following analytical forms:

Because the derivative of the error function is Gaussian, we can also easily integrate ϕ(x)\phi^{\prime}(x) over the Gaussian distribution and obtain the following analytical representations of the recurrence relations (11):

To compute the recurrence relations for the feedforward correlations (10), note that we can generally transform Iϕ[a,b]I_{\phi}[a,b] into

This is the analytical form of the recurrence relation for q^stl\hat{q}^{l}_{st}.

Finally, because the derivative of the error function is Gaussian, we can also easily obtain

B.2 ReLU networks

We define a ReLU activation as ϕ(x)=0  (x<0),  x  (0x)\phi(x)=0\ \ (x<0),\ \ x\ \ (0\leq x). For a network with this ReLU activation function, the recurrence relations for the macroscopic variables require no numerical integrations.

B.3 Linear networks

We define a linear activation as ϕ(x)=x\phi(x)=x. For a network with this linear activation function, the recurrence relations for the macroscopic variables do not require numerical integrations.

Appendix C Additional Experiments

C.2 Training on CIFAR-10