On the Variance of the Fisher Information for Deep Learning
Year: 2021

Abstract:
In the realm of deep learning, the Fisher information matrix (FIM) gives novel insights and useful tools to characterize the loss landscape, perform second-order optimization, and build geometric learning theories. The exact FIM is either unavailable in closed form or too expensive to compute. In practice, it is almost always estimated based on empirical samples. We investigate two such estimators based on two equivalent representations of the FIM — both unbiased and consistent. Their estimation quality is naturally gauged by their variance given in closed form. We analyze how the parametric structure of a deep neural network can affect the variance. The meaning of this variance measure and its upper bounds are then discussed in the context of deep learning.

heading: 1 Introduction
text: The Fisher information is one of the most fundamental concepts in statistical machine learning. Intuitively, it measures the amount of information carried by a single random observation when the underlying model varies along certain directions in the parameter space: if such a variation does not change the underlying model, then a corresponding observation contains zero (Fisher) information and is non-informative regarding the varied parameter. Parameter estimation is impossible in this case. Otherwise, if the variation significantly changes the model and has large information, then an observation is informative and the parameter estimation can be more efficient as compared to parameters with small Fisher information. In machine learning, this basic concept is useful for defining intrinsic structures of the parameter space, measuring model complexity, and performing gradient-based optimization.
Given a statistical model that is specified by a parametric form p(z |θ) and a continuous domain θ ∈ M, the Fisher information matrix (FIM) is a 2D tensor varying with θ ∈ M, given by
I(θ) = Ep(z | θ) ( ∂ℓ
∂θ
∂ℓ
∂θ⊤
) , (1)
where Ep(z | θ)(·), or simply Ep(·) if the model p is clear from the context, denotes the expectation w.r.t. p(z |θ), and ℓ := log p(z |θ) is the log-likelihood function. All vectors are column vectors throughout this paper. Under weak conditions (see Lemma 5.3 in Lehmann and Casella [16] for the univariate case), the FIM has the equivalent expression I(θ) = Ep(z | θ) ( −∂2ℓ/∂θ∂θ⊤ ) . Given N i.i.d. observations z1, . . . , zN , these two equivalent expressions of the FIM lead to two different estimators
Î1(θ) = 1
N N∑ i=1 ( ∂ℓi ∂θ ∂ℓi ∂θ⊤ ) and Î2(θ) = 1 N N∑ i=1 ( − ∂ 2ℓi ∂θ∂θ⊤ ) , (2)
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
where ℓi := log p(zi |θ) is the log-likelihood of the i’th observation zi. The notations Î1(θ) and Î2(θ) are abused for simplicity as they depend on both θ and the random observations zi. These estimators are universal and independent to the parametric form p(z |θ). They are expressed in terms of the 1st- or 2nd-order derivatives of the log-likelihood. Usually, we already have these derivatives to perform gradient-based learning. Therefore, we can save computational cost and reuse these derivatives to estimate the Fisher information, which in turn can be useful, e.g., to perform natural gradient optimization [1, 25]. Estimating the FIM is especially meaningful for deep learning, where the computational overhead of the exact FIM can be significant.
It is straightforward from the law of large numbers and the central limit theorem that both estimators in Eq. (2) are unbiased and consistent. This is formally stated as follows.
Proposition 1.
Ep(z | θ) ( Î1(θ) ) = Ep(z | θ) ( Î2(θ) ) = I(θ).
∀ϵ > 0, lim N→∞
Prob (∥∥∥Î1(θ)− I(θ)∥∥∥ F + ∥∥∥Î2(θ)− I(θ)∥∥∥ F > ϵ ) = 0,
where Prob(·) denotes the probability of the parameter statement being true and ∥ · ∥F is the Frobenius norm of a tensor (with ∥ · ∥2 as the regular vector L2-norm)
The Fisher information can be zero for non-regular models or infinite [7]. However, these properties may not be preserved by the empirical estimators.
How far can Î1(θ) and Î2(θ) deviate from the “true FIM” I(θ), and how fast can they converge to I(θ) as the number of observations increases? To answer these questions, it is natural to think of the variance of Î1(θ) and Î2(θ). For example, an estimator with a large variance means the estimation does not accurately reflect I(θ); and any procedure depending on the FIM consequently suffers from the estimation error. Through studying the variance, we can control the estimation quality and reliably perform subsequent measurements or algorithms based on the FIM.
Towards this direction, we made the following contributions that will unfold in the following Sections 2 to 4:
• We review and rediscover two equivalent expression of the FIM in the context of deep feed-forward networks (Section 2);
• We give in closed form the variance (extending to meaningful upper bounds) and discuss the convergence rate of the estimators Î1(θ) and Î2(θ) (Section 3);
• We analyze how the 1st- and 2nd-order derivatives of the neural network can affect the estimation of the FIM (Section 4).
We discuss related work in Section 5 and conclude in Section 6.

heading: 2 Feed-forward Networks with Exponential Family Output
text: This section realizes the concept of Fisher information in a feed-forward network with exponentially family output and explains why its estimators are useful in theory and practice.
Consider supervised learning with a neural network. The underlying statistical model is p(z |θ) = p(x)p(y |x,θ), where z = (x,y), the random variable x represents features, and y is the target variable. The marginal distribution p(x) is parameter-free, usually fixed as the empirical distribution p(x) = 1M ∑M i=1 δ(x−xi) w.r.t. a set of observations {xi}Mi=1, where δ(·) is the Dirac delta. In this paper, we consider w.l.o.g. M = 1 as the FIM w.r.t. observations {xi}Mi=1 is simply the average over FIMs of each individual observation. All results generalize to multiple observations by taking the empirical average. The predictor p(y |x,θ) is a neural network with parameters θ = {Wl−1}Ll=1
and exponential family output units, given by p(y |x) = exp ( t⊤(y)hL − F (hL) ) ,
hL = WL−1h̄L−1,
hl = σ(Wl−1h̄l−1), (l = 1, . . . , L− 1) h̄l = (h ⊤ l , 1) ⊤,
h0 = x, (3) where t(y) is the sufficient statistics of the prediction model, F (h) = log ∫ exp(t⊤(y)h)dy is the log-partition function, and σ : ℜ → ℜ is an element-wise non-linear activation function. Moreover, Wl is a nl+1 × (nl + 1) matrix, representing the neural network parameters (weights and biases) in the l’th layer, where nl := dim(hl) denotes the size of layer l. We use W−l for the nl+1 × nl weight matrix without the bias terms, obtained by removing the last column of Wl. hl is a learned representation of the input x. All intermediate variables hl are extended to include a constant scalar 1 in h̄l, so that a linear layer can simply be expressed as Wlh̄l. The last layer’s output hL with dimensionality nL specifies the natural parameter of the exponential family.
We need the following Lemma which gives the FIM w.r.t. hL, which is a nL × nL matrix in simple closed form for commonly used probability distributions. Lemma 2. For the neural network model specified in Eq. (3),
I(hL) = Cov(t(y)) = ∂η
∂hL ,
where Cov(·) denotes the covariance matrix w.r.t. p(y |x,θ), η := η(hL) := ∂F/∂hL is the expectation parameters, and the vector-vector-derivative ∂η/∂hL denotes the Jacobian matrix of the mapping hL → η.
The derivatives of the log-likelihood ℓ(θ) := log p(x,y |θ) characterize its landscape and are essential to compute the FIM. By Eq. (3), the score function (gradient of ℓ) is
∂ℓ ∂θ = ( ∂hL ∂θ )⊤ (t(y)− η(hL)) = ∂haL ∂θ (ta − ηa). (4)
In this paper, we mix the usual Σ-notation of summation with the Einstein notation: in the same term, an index appearing in both upper- and lower-positions indicates a sum over this index. For example, tah a = ∑
a taha. Hence, in our equations, upper- and lower-indexes have the same meaning: both ha and ha mean the a’th element of h. For convenience and consistency, we take quantities w.r.t. θ as upper indexed and other quantities as lower indexed, i.e., Iij(θ) versus Iij(hL). This mixed representation of sums helps to simplify our expressions without causing confusion. From Eq. (4) and Lemma 2, the Hessian of ℓ is given by
∂2ℓ
∂θ∂θ⊤ = (ta − ηa) ∂2haL ∂θ∂θ⊤ − ∂h a L ∂θ ∂ηa ∂θ⊤ = (ta − ηa) ∂2haL ∂θ∂θ⊤ − ∂h a L ∂θ Iab(hL) ∂hbL ∂θ⊤ . (5)
Similar to the case of a general statistical model, the FIM is equivalent to the expectation of the Hessian of −ℓ as long as the activation function is smooth enough. Theorem 3. Consider the neural network model in Eq. (3). For any activation function σ ∈ C2(ℜ) (both σ′(z) and σ′′(z) exist and are continuous), we have I(θ) = Ep ( − ∂ 2ℓ ∂θ∂θ⊤ ) .
Remark 3.1. ReLU networks do not have this equivalent expression as ReLU(z) is not differentiable at z = 0.
Through the definition of the FIM, or alternatively its equivalent formula in Theorem 3, we arrive at the same expression
I(θ) = ( ∂hL ∂θ )⊤ I(hL) ∂hL ∂θ⊤ = ∂haL ∂θ Iab(hL) ∂hbL ∂θ⊤ . (6)
Equation (6) takes the form of a generalized Gauss-Newton matrix [18]. This general expression of the FIM has been known in the literature [24, 25]. Under weak conditions, I(θ) is a pullback
metric [30] of I(hL) in Lemma 2 associated with the mapping θ → hL. To compute I(θ) in closed form, one need to first compute the Jacobian matrix of size nL × dim(θ) then perform the matrix multiplication in Eq. (6). The naive algorithm to evaluate Eq. (6) has a computational complexity of O(n2L dim(θ)+nL dim
2(θ)), where the term O(nL dim2(θ)) is dominant as dim(θ) ≫ nL in deep architectures. Once the parameter θ is updated, the FIM has to be recomputed. This is infeasible in practice for large networks where dim(θ) can be millions or billions.
The two estimators Î1(θ) and Î2(θ) in Eq. (2) provide a computationally inexpensive way to estimate the FIM. Given θ and x, one can draw i.i.d. samples y1, . . . ,yN ∼ p(y |x,θ). Both ∂ℓi/∂θ and ∂2ℓi/∂θ∂θ⊤ can be evaluated directly through auto-differentiation (AD) that is highly optimized for modern GPUs. For Î1(θ), we already have ∂ℓi/∂θ to perform gradient descent. For Î2(θ), efficient methods to compute the Hessian are implemented in AD frameworks such as PyTorch [26]. Using these derivatives, the computational cost only scales with the number N of samples but does not scale with nL.
We rarely need the full FIM of size dim(θ)×dim(θ). Most of the time, only its diagonal blocks are needed, where each block corresponds to a subset of parameters, e.g., the neural network weights of a particular layer. Therefore the computation of both estimators can be further reduced.
If p(y |x,θ) has the parametric form in Eq. (3), from Eqs. (4) and (5), the FIM estimators become
Î1(θ) = ∂haL ∂θ · 1 N N∑ i=1 (ta(yi)− ηa)(tb(yi)− ηb) · ∂hbL ∂θ⊤ , (7)
Î2(θ) = ( ηa − 1
N N∑ i=1 ta(yi) ) ∂2haL ∂θ∂θ⊤ + ∂haL ∂θ Iab(hL) ∂hbL ∂θ⊤ . (8)
Recall that the notation of Î1(θ) and Î2(θ) is abused as they depend on x and y1 · · ·yN . Notably, in Eq. (7), Î1(θ) is expressed in terms of the Jacobian matrix of the mapping θ → hL and the empirical variance of the minimal sufficient statistic t(yi) of the output exponential family. In Eq. (8), Î2(θ) depends on both the Jacobian and Hessian of θ → hL and the empirical average of t(yi). The second term on the right-hand side (RHS) of Eq. (8) is exactly the FIM, and therefore the first term serves as a bias term. Eqs. (7) and (8) are only for the case with exponential family output. If the output units belong to non-exponential families, e.g., a statistical mixture model, one falls back to the general formulae, i.e., Eq. (2) for the FIM.
As an application of the Fisher information, the Cramér-Rao lower bound (CRLB) states that any unbiased estimator θ̂ of the parameters θ satisfies Cov(θ̂) ⪰ [I(θ)]−1. For example, in Lemma 2, the FIM is w.r.t. the output of the neural network. As such, I(hL) can be used to study the estimation covariance of hL based on random samples y1 · · ·yN drawn from p(y |x,θ). Similarly for Eq. (6), we can consider unbiased estimators of the weights of the neural network. In any case, to apply the CRLB, one needs an accurate estimation of I(θ). If the scale of I(θ) is relatively small when compared to its covariance, its estimation Î(θ) is more likely to be a small positive value (or even worse, zero or negative). The empirical computation of the CRLB is not meaningful in this case.

heading: 3 The Variance of the FIM Estimators
text: Based on the deep learning architecture specified in Eq. (3), we measure the quality of the two estimators Î1(θ) and Î2(θ) given by their variances. Given the same sample size N , a smaller variance is preferred as the estimator is more accurate and likely to be closer to the true FIM I(θ). We study how the structure of the exponential family has an impact on the variance.

heading: 3.1 Variance in closed form
text: We first consider Î1(θ) and Î2(θ) in Eq. (2) as real matrices of dimension dim(θ)× dim(θ). As Î1(θ) is a square matrix, the corresponding covariance is a 4D tensor [ Cov ( Î1(θ) )]ijkl of dimen-
sion dim(θ)× dim(θ)× dim(θ)× dim(θ), representing the covariance between the two elements
Îij1 (θ) and Îkl1 (θ). The element-wise variance of Î1(θ) is a matrix with the same size of Î1(θ), which we denote as Var(Î1(θ)). Thus,
Var(Î1(θ))ij = [ Cov ( Î1(θ) )]ijij . (9)
Similarly, the covariance and element-wise variance of Î2(θ) are denoted as [ Cov ( Î2(θ) )]ijkl and
Var(Î2(θ))ij , respectively. As the samples y1, . . . ,yN are i.i.d., we have
Cov(Î1(θ)) = 1
N Cov
( ∂ℓ
∂θ
∂ℓ
∂θ⊤
) and Cov(Î2(θ)) = 1
N Cov
( ∂2ℓ
∂θ∂θ⊤
) . (10)
Both Cov(Î1(θ)) and Cov(Î2(θ)) have an order of O(1/N). For the neural network model in Eq. (3), we further have those covariance tensors in closed form. Theorem 4.[ Cov ( Î1(θ) )]ijkl = 1
N · Cov
( ∂ℓ
∂θi
∂ℓ ∂θj , ∂ℓ ∂θk ∂ℓ ∂θl ) = 1
N · ∂ihaL(x)∂jhbL(x)∂khcL(x)∂lhdL(x) · (Kabcd(t)− Iab(hL) · Icd(hL)) ,
where the 4D tensor
Kabcd(t) := E [(ta − ηa(hL(x)))(tb − ηb(hL(x)))(tc − ηc(hL(x)))(td − ηd(hL(x)))] is the 4th (unscaled) central moment 1 of t(y) and ∂ihL(x) := ∂hL(x)/∂θi. Remark 4.1. The 4D tensor (Kabcd(t)− Iab(hL) · Icd(hL)) is the covariance of the random matrix
∂ℓ
∂hL
∂ℓ
∂h⊤L = (t(y)− η)(t(y)− η)⊤,
where y ∼ p(y |hL). This random matrix is an estimator of I(hL), i.e. the FIM w.r.t. the natural parameters hL. Theorem 4 describes how the covariance tensor adapts w.r.t. the coordinate transformation hL → θ.
Notably, as t(y) is the sufficient statistics of an exponential family, the derivatives of the logpartition function F (h) w.r.t. the natural parameters h are equivalent to the cumulants of t(y). The cumulants correspond to the coefficients of the Taylor expansion of the logarithm of the moment generating function [20]. Importantly, the cumulants of order 3 and below are equivalent to the central moments (see e.g. Lemma 2). However, this is not the case for the 4th central moment which must be expressed as a combination of the 2nd and 4th cumulants, as stated in the following Lemma. Lemma 5.
Kabcd(t) = κabcd + Iab(hL) · Icd(hL) + Iac(hL) · Ibd(hL) + Iad(hL) · Ibc(hL), where
κabcd := ∂4F (h)
∂ha∂hb∂hc∂hd ∣∣∣∣ h=hL(x) .
Remark 5.1. In the 1D case, the 4th central moment simplifies to K(t) = F ′′′′(hL)+ 3(F ′′(hL))2.
For the second estimator Î2(θ), the covariance only depends on the 2nd central moment of t(y). Theorem 6.[
Cov ( Î2(θ) )]ijkl = 1
N · Cov
( − ∂ 2ℓ
∂θi∂θj ,− ∂
2ℓ
∂θk∂θl
) = 1
N · ∂2ijhαL(x)∂2klh β L(x)Iαβ(hL),
where ∂2ijhL(x) := ∂ 2hL(x)/∂θi∂θj 2. 1The kurtosis of a random variable is defined by its 4th standardized (both centered and normalized) moment. Here, K(·) denotes the 4th central moment but not the kurtosis. 2In this paper, the derivatives are by default taken w.r.t. θ. Therefore, ∂i := ∂∂θi and ∂ 2 ij := ∂2 ∂θi∂θj .
Remark 6.1. By Lemma 2, the matrix Iαβ(hL) is the covariance of the sufficient statistic t(y). Hence, the covariance of Î2(θ) scales with the covariance of t(y). If t(y) tends to be deterministic, then the covariance of Î2(θ) shrinks towards 0 and its estimation of the FIM becomes accurate.
The covariance in Theorems 4 and 6 has two different components: À the derivatives of the deep neural network; and Á the central (unscaled) moments of t(y). The 4D tensor Kabcd(t) and the 2D FIM Iαβ(hL) correspond to the 4th and 2nd central moments of t(y), respectively. Intuitively, the larger the scale of the Jacobian or the Hessian of the neural network mapping θ → hL and/or the larger the central moments of the exponential family, the lower the accuracy when estimating the FIM.

heading: 3.2 Variance Bounds
text: We aim to derive meaningful upper bounds of the covariances presented in Theorems 4 and 6. Using the Cauchy-Schwarz inequality, we can “decouple” the derivatives of the neural network mapping and the central moments of the exponential family into different terms. This provides various bounds on the scale of covariance quantities.
Theorem 7. ∥∥∥Cov (Î1(θ))∥∥∥ F ≤ 1 N · ∥∥∥∥∂hL∂θ ∥∥∥∥4 F · ∥K(t)− I(hL)⊗ I(hL)∥F ,
where ⊗ is the tensor-product: (I(hL)⊗ I(hL))abcd := Iab(hL) · Icd(hL). The scale ∥∥∥Cov (Î1(θ))∥∥∥
F measures how much the estimator Î1(θ) deviates from I(θ). Theo-
rem 7 says that this deviation is bounded by the scale of the Jacobian matrix ∂hL/∂θ as well as the scale of (K(t) − I(hL) ⊗ I(hL)). Recall from Remark 4.1 the latter measures the variance when estimating the FIM I(hL) of the exponential family. Theorem 7 allows us to study these two different factors separately. Similarly, we have an upper bound on the scale of the covariance of Î2(θ). Theorem 8. ∥∥∥Cov (Î2(θ))∥∥∥
F ≤ 1 N · ∥∥∥∥∂2hL(x)∂θ∂θ⊤ ∥∥∥∥2 F · ∥I(hL)∥F .
On the RHS, the Hessian ∂2hL(x)/∂θ∂θ⊤ is a 3D tensor of shape nL × dim(θ)× dim(θ). Therefore, the variance of Î2(θ) is bounded by the scale of the Hessian, as well as the scale of the FIM I(hL) of the output exponential family. We consider an upper bound to further simplify related terms in Theorems 7 and 8.
Lemma 9.
∥K(t)− I(hL)⊗ I(hL)∥F ≤ √ 2 ( nL∑ a=1 (√ Kaaaa(t) + Iaa(hL) ))2 ,
∥I(hL)∥F ≤ nL∑ a=1 Iaa(hL).
Remark 9.1. Using Lemma 9, it is straightforward to bound the scale of the covariance tensors with the size of the Jacobian/Hessian, as well as the central moments Kaaaa(t) and Iaa(hL). These bounds are meaningful but omitted for brevity.
Remark 9.2. By Lemma 9, ∥K(t)− I(hL)⊗ I(hL)∥F is in the order of O(n2L) and ∥I(hL)∥F is in the order of O(nL).
The scale of the tensors K(t)− I(hL)⊗ I(hL) and I(hL) is bounded by the diagonal elements of K(t) and I(hL), or the element-wise central moments of t(y). Understanding the scale of these 1D central moments helps to understand the scale of the moment terms in our key statements.
Table 1 presents some 1D exponential families and their cumulants. Figure 1 displays K(t)−Var2(t) and Var(t) against the mean of these distributions. Based on Fig. 1a, if the neural network has Bernoulli output units, then the scale of Kaaaa(t)− (Iaa(hL))2 is smaller than Iaa(hL) regardless of hL. Notably, when p = 0.5, the variance of the first estimator Î1(θ) is 0 — regardless of hL. For normal distribution output units (corresponding to the mean squared error loss) in Fig. 1b, both central moment quantities are constant. For Poisson output units in Fig. 1c, Iaa(hL) increases linearly with the average number of events λ, while Kaaaa(t)− (Iaa(hL))2 increases quadratically. Thus, the upper bound of ∥Cov(Î1(θ))∥F increases faster than the upper bound of ∥Cov(Î2(θ))∥F as hL enlarges. In this case, one may prefer Î2(θ) rather than Î1(θ) and/or control the scale of hL. In general, hL is desired to be in certain regions in the parameter space of the exponential family to control the estimation variance of the FIM. Techniques to achieve this include regularization on the scale of hL; temperature scaling [12]; or normalization layers [4, 28]. Of course, they could inversely increase the scale of the derivatives of the neural network, which can be controlled by imposing additional constraints, i.e., Lipschitz requirements.

heading: 3.3 Positive Definiteness
text: By definition, the FIM of any statistical model is positive semidefinite (p.s.d.). The first estimator Î1(θ) is naturally on the p.s.d. manifold (space of p.s.d. matrices). On the other hand, Î2(θ) can “fall off” the p.s.d. manifold. It is important to examine the likelihood for Î2(θ) having a negative spectrum and the corresponding scale, so that any algorithm (e.g. natural gradient) relying of the FIM being p.s.d. can be adapted.
Eq. (8) can be re-expressed as the sum of a p.s.d. matrix and a linear combination of nL symmetric matrices. We provide the likelihood for Î2(θ) staying on the p.s.d. manifold given conditions on the spectrum of the Hessian. Theorem 10. Let λmin(·), λmax(·), and ρ(·) denote the smallest eigenvalue, the largest eigenvalue, and the spectral radius (largest absolute value of the spectrum), respectively. Let ρ := (ρ(∂2h1L), · · · ρ(∂2h nL L )). If λmin(I(θ)) > 0, then with probability at least
1− nL · ∥ρ∥ 2 2 · λmax(I(hL))
N · λ2min(I(θ)) ,
the estimator Î2(θ) with N samples is a p.s.d. matrix.
The bound becomes uninformative as the output layer size nL increases, as the spectrum of the Hessian of hL scales up, or as the spectrum of the FIM I(hL) enlarges. On the other hand, as the
minimal eigenvalue of the FIM I(θ) increases, Theorem 10 can give meaningful lower bounds. In particular, with sample rate O(N−1), estimator Î2(θ) will be a p.s.d. matrix. In practice for overparametrized networks, λmin(I(θ)) is close to or equals 0 and Theorem 10 is not meaningful. In any case, we need to consider the scale of the negative spectrum of Î2(θ). Theorem 11.
λmin ( Î2(θ) ) ≥ −ρ(∂2haL(x)) ∣∣∣∣∣ηa − 1N N∑ i=1 ta(yi) ∣∣∣∣∣ . Theorem 11 guarantees that in the worst case, the scale of the negative spectrum of Î2(θ) is controlled. By Lemma 2, Var(ηa − 1N ∑N i=1 ta(yi)) = 1 N I
aa(hL). Therefore, as N increases or Iaa(hL) decreases, the negative spectrum of Î2(θ) will shrink. Further analysis on the spectrum of Î1(θ) and Î2(θ) can utilize the geometric structure of the p.s.d. manifold. This is left for future work.

heading: 3.4 Convergence Rate
text: The rate of convergence for each of the estimators is of particular interest when considering their practical viability. Through a generalized Chebyshev inequality [8], we can get a simple Frobenius norm convergence rate. Lemma 12. Let 0 < ε < 1. Then∥∥∥Î1(θ)− I(θ)∥∥∥
F ≤ 1√ εN · √√√√dim(θ)∑ i,j=1 Var ( ∂ℓ ∂θi ∂ℓ ∂θj ) holds with probability at least 1− ε; and
∥∥∥Î2(θ)− I(θ)∥∥∥ F ≤ 1√ εN · √√√√dim(θ)∑ i,j=1 Var ( − ∂ 2ℓ ∂θi∂θj ) hold with probability at least 1− ε.
Each of these convergence rates only depends on the element-wise variance of the estimator terms in Eq. (9). Moreover, each of the estimators has a convergence rate of O(N−1/2). The rate’s constants are determined by the variance of the estimators given by Theorems 4 and 6, which are influenced by the derivatives of the neural network and the moments of the output exponential family.

heading: 4 Effect of Neural Network Derivatives
text: The derivatives of the deep learning network can affect the estimation variance of the FIM. By Theorem 4, the variance of the first estimator Î1(θ) scales with the Jacobian of the neural network mapping θ → hL(x). By Theorem 6, the variance of Î2(θ) scales with the Hessian of θ → hL(x). The larger the scale of the Jacobian or the Hessian, the larger the estimation variance. In this section, we examine these derivatives in more detail.
We give the closed form gradient of the log-likelihood ℓ and the last layer’s output hL w.r.t. the neural network parameters. Lemma 13.
∂ℓ
∂Wl = Dl
∂ℓ
∂hl+1 h̄⊤l ,
∂ℓ
∂hl = B⊤l (t(y)− η(hL)) , ∂haL ∂Wl = DlB ⊤ l+1eah̄ ⊤ l ,
where ea is the ath standard basis vector, Bl and Dl are recursively defined by
BL = I, Bl = Bl+1DlW − l , DL−1 = I, Dl = diag ( σ′(Wlh̄l) ) ,
I is the identity matrix, and diag(·) means a diagonal matrix with given diagonal entries.
By Lemma 13, we can estimate the FIM w.r.t. the hidden representations hl through
Î1(hl) = 1
N N∑ i=1 ∂ℓi ∂hl ∂ℓi ∂h⊤l = B⊤l
( 1
N N∑ i=1
(t(yi)− η(hL)) (t(yi)− η(hL))⊤ ) Bl. (11)
As Bl is recursively evaluated from the last layer to previous layers, the FIM can also be recursively estimated based on Î1(θ). It is similar to back-propagation, except that the FIMs are backpropagated instead of gradients of the network. This is similar to the backpropagated metric [23].
To investigate how the first estimator Î1(θ) is affected by the loss landscape, we bound the Frobenius norm of the parameter-output Jacobian ∂hL/∂θ.
Lemma 14. If the activation function has bounded gradient and ∀z ∈ ℜ, |σ′(z)| ≤ 1, then∥∥∥∥ ∂hL∂Wl ∥∥∥∥ F = ∥Bl+1Dl∥F · ∥h̄l∥2 ≤ L−1∏ i=l+1 ∥W−i ∥F · ∥h̄l∥2, (12) where ∂hL/∂Wl = [ ∂h1L/∂Wl, · · · , ∂h nL L /∂Wl ] is the derivative of a vector w.r.t. a matrix that is a 3D tensor.
Given Lemma 14, we see that the gradient ∂hL/∂Wl scales with both the neural network weights Wi and the gradient of the activation function Dl. Common activation functions have both bounded outputs and 1st-order derivatives; or at least are locally Lipschitz, i.e., sigmoid and ReLU activation functions. During training, regularizing the scale of the neural network weights is a sufficient condition for bounding the variance of Î1(θ). An alternative bound can be established which depends on the maximum singular values of the weight matrices.
Lemma 15. Suppose that the activation function has bounded gradient ∀z ∈ ℜ, |σ′(z)| ≤ 1. Then∥∥∥∥ ∂hL∂Wl ∥∥∥∥ 2σ ≤ ( L−1∏ i=l+1 smax(W − i ) ) · ∥h̄l∥2, (13)
where smax(·) denotes the maximum singular value and ∥T ∥2σ denotes the tensor spectral norm for a 3D tensor T , defined by
∥T ∥2σ = max {⟨T ,α⊗ β ⊗ γ⟩ : ∥α∥2 = ∥β∥2 = ∥γ∥2 = 1} .
Therefore, regularizing smax(W−i ), or the spectral norm of the weight matrices, also helps to improve the estimation accuracy of the FIM.
We further reveal the relationship between the loss landscape and the FIM estimators. For a given target ỹ, the log-likelihood is denoted as l̃ := log p(ỹ |x,θ). Furthermore, let us define ∆Î1(θ) := (∂l̃/∂θ)(∂l̃/∂θ⊤)− Î1(θ) and ∆Î2(θ) := −∂2 l̃/∂θ∂θ⊤ − Î2(θ). By Eqs. (4) and (5),
∆Î1(θ) = ∂haL ∂θ
[ (ta(ỹ)− ηa)(tb(ỹ)− ηb)− 1
N N∑ i=1 (ta(yi)− ηa)(tb(yi)− ηb) ] ∂hbL ∂θ⊤ ,
∆Î2(θ) =
[ 1
N N∑ i=1 ta(yi)− ta(ỹ) ] ∂2haL ∂θ∂θ⊤ .
Hence, the difference between Î1(θ) (resp. Î2(θ)) and the squared gradient (resp. Hessian) of the loss −ℓ̃ depends on how yi differs from ỹ. If the network θ is trained, then the random samples yi ∼ p(y |x,θ) are close to the given target ỹ. In this case, Î1(θ) corresponds to the squared gradient, and Î2(θ) corresponds to the Hessian. This is not true for untrained neural networks with random weights.

heading: 5 Related Work
text: The two estimators Î1(θ) and Î2(θ) are not new as one usually utilizes one of them to compute the FIM. Guo and Spall [11] analyzed their accuracy for univariate symmetric density functions based on the central limit theorem. Fisher information estimation is also examined in latent variable models [9]. Under the same topic, our work analyzes the factors affecting the variance of Î1(θ) and Î2(θ) for deep neural networks. A large body of work tries to approximate the FIM or define similar curvature tensors for performing natural gradient descent [1, 18, 19, 14, 23, 31]. If the loss is an empirical expectation of − log p(y |x,θ), its Hessian w.r.t. hL is exactly I(hL). Then, the FIM in Eq. (6) is in the form of a Generalized Gauss-Newton matrix (GNN) [29, 18]. Ollivier [23] provided algorithm procedures to compute the unit-wise FIM and discusses Monte Carlo natural gradient. The FIM can be computed locally [31, 3] based on a joint distribution representation of the neural network. Efficient computational methods are developed to evaluate the FIM inverse [24].
The estimator Î1(θ) is not the “empirical Fisher” (see e.g. [18, Section 11]) as yi is randomly sampled from p(y |x,θ), making Î1(θ) an unbiased estimator. The difference between the empirical Fisher and the FIM is clarified [15]. Similarly, the estimator Î2(θ) is not the Hessian of the loss, as yi is randomly sampled rather than fixed to the given target.
Recently, the structure of the FIM (or its partial approximations) are examined in deep learning. The FIM of randomized networks is analyzed [2], where the weights of the neural network are assumed to be random. Often the analysis of randomized networks uses spectral analysis and random matrix theory [27]. An insight from this body of work is that most of the eigenvalues of the FIM are close to 0; while the high end of spectrum has large values [13]. In this paper, the estimators Î1(θ) and Î2(θ) are random matrices due to the sampling of yi ∼ p(y |x) (the weights are considered fixed). In information geometry [1, 22], the FIM serves as a Riemannian metric in the space of probability distributions. The FIM is a covariant tensor and is invariant to diffeomorphism on the sample space [22]. Higher order tensors are used to describe the intrinsic structure in this space. For example, the Riemannian curvature is a 4D tensor, while the Ricci curvature is 2D. The third cumulants of the sufficient statistics give an affine connection (belonging to the α-connections or the AmariČensov tensor) of the exponential family. The FIM is generalized to a one-parameter family [21].
In statistics, our estimator Î2(θ) is Efron and Hinkley [10]’s “observed Fisher information”, which is usually evaluated at the maximum likelihood estimation θ̂. Higher order moments of the maximum likelihood estimator (MLE) were discussed (see e.g. Bowman and Shenton [5]). These moments are associated with parameter estimators and differ from the concept of the FIM estimators. Similar concepts are examined in higher-order asymptotic theory [1, Chapter 7].

heading: 6 Conclusion
text: The FIM is a covariant p.s.d. tensor revealing the intrinsic geometric structure of the parameter manifold. It yields useful practical methods such as the natural gradient. In practice, the true FIM I(θ) is usually expensive or impossible to obtain. Estimators of the FIM based on empirical samples is used in the deep learning practice. We analyzed two different estimators Î1(θ) and Î2(θ) of the FIM of a deep neural network. These estimators are convenient to compute using auto-differentiation frameworks but randomly deviates from I(θ). Our central results, Theorems 4 and 6, present the variance of Î1(θ) and Î2(θ) in closed form, which is further extended to upper bounds in simpler forms. Two factors affecting the estimation variance are À the derivatives of neural network output hL w.r.t. the weight parameters θ; and Á the property of hL as an exponential family distribution. A large scale of the 1st- and/or 2nd-order derivatives leads to a large variance when estimating the FIM. Our analytical results can be useful to measure the quality of the estimated FIM and could lead to variance reduction techniques.
Acknowledgments and Disclosure of Funding
We thank the anonymous NeurIPS reviewers for their constructive comments. We thank Frank Nielsen for the insightful feedback. We thank James C. Spall for pointing us to related work.

Review summaries:
Update: I thank the authors for their response, and I will keep the score as is. It would be nice to have a simple empirical analysis in the final supplemental material.   The paper discusses the covariance of two FIM estimators, which is useful in deep learning where the exact FIM over the full dataset can be hard to compute. One is based on the variance of the score (type 1) and the other is based on the Hessian of the relative entropy (type 2). Comparing these covariances can depend on the exponential family used in the model, type 1 has a higher variance for Normal and Poisson whereas type 2 has a higher variance for Bernoulli. Convergence rate results are given which is $O(N^{-0.5))$ with constant determined by the variances. These results provide some theoretical justification as to which estimator of FIM to use when using it to perform natural gradient in deep learning. 

This paper investigates the Fisher information matrix (FIM) theoretically. Two equivalent expressions of FIM and their estimators are discussed (rediscovered) for deep neural network with exponential family output units. The closed-form variance (and simpler-form upper bounds) for the two estimators is proposed; the convergence rate of the positiveness definiteness are also analyzed. Finally, the variance is connected with the structure of the corresponding deep neural networks.