Trade-Offs of Diagonal Fisher Information Matrix Estimators

Authors: alexander  soen
ke  sun
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-13 13:16:07

Abstract: The Fisher information matrix can be used to characterize the local geometry of
the parameter space of neural networks. It elucidates insightful theories and
useful tools to understand and optimize neural networks. Given its high
computational cost, practitioners often use random estimators and evaluate only
the diagonal entries. We examine two popular estimators whose accuracy and sample
complexity depend on their associated variances. We derive bounds of the
variances and instantiate them in neural networks for regression and
classification. We navigate trade-offs for both estimators based on analytical
and numerical studies. We find that the variance quantities depend on the
non-linearity w.r.t. different parameter groups and should not be neglected when
estimating the Fisher information.

Main paper: 
 settings: 
 In the parameter space of neural networks (NNs), i.e. the neuromanifold [1] , the network weights and biases play the role of a coordinate system and the local metric tensor can be described by the Fisher Information Matrix (FIM). As a result, empirical estimation of the FIM helps reveal the geometry of the loss landscape and the intrinsic structure of the neuromanifold. Utilizing these insights has lead to efficient optimization algorithms, e.g., the natural gradient [1] and Adam [16] .A NN with inputs x and stochastic outputs y can be specified by a conditional p.d.f. p(y | x; θ), where θ is the NN's weights and biases. This paper considers the general parametric formEQUATIONwhere h θ : ℜ I → ℜ T maps I-dimensional inputs x to T -dimensional exponential family parameters, t(y) is a vector of sufficient statistics, π(y) is a base measure, and F (•) is the log-partition function (normalizing the exponential). For example, if y denotes class labels and t(y) maps to its corresponding one-hot vectors, then Eq. ( 1) is associated with a multi-class classification network.Assuming that the marginal distribution q(x) is parameter-free, we define parametric joint distributions p(x, y; θ) = q(x)p(y | x; θ). The (joint) FIM is defined as I(θ) .is the 'conditional FIM'. The second equality (*) holds if h θ 's activation functions are in C 2 (ℜ) (i.e., h θ is a sufficiently smooth NN). I(θ | x) does not have this equivalent expression (*) for NNs with ReLU activation functions [37] . Both I(θ) and I(θ | x) define dim(θ) × dim(θ) positive semi-definite (PSD) matrices. The distinction in notation is to emphasize that the joint FIM I(θ) (depending only on θ) is simply the average over individual conditional FIMs I(θ | x) (depending on both θ and x).38th Conference on Neural Information Processing Systems (NeurIPS 2024). In practice, the FIM is typically computationally expensive and needs to be estimated. Given q(x) and a NN with weights and biases θ parameterizing p(y | x; θ), as per Eq. ( 1), we consider two commonly used estimators of the FIM [11, 37] given byEQUATIONwhere p(y k | x k ) .= p(y k | x k ; θ) and (x 1 , y 1 ), . . . , (x N , y N ) are i.i.d. sampled from p(x, y; θ). A conditional variant of the estimators, denoted as Î1 (θ | x) and Î2 (θ | x), can be defined by fixing x = x 1 = • • • = x N and sampling y 1 , . . . , y N independently from p(y | x) in Eq. ( 3) -details omitted for brevity.Both estimators, Î1 (θ) and Î2 (θ), are random matrices with the same shape as I(θ). By Eq. ( 2), they are unbiased -for Î2 (θ), this only holds if activations functions are in C 2 (ℜ). Following Eq. ( 1)'s setting, the estimation variances of Î1 (θ) and Î2 (θ) can be expressed in closed form and upper bounded [37] . This provides an important, yet not widely discussed, tool for quantifying the estimators' accuracy [11] and hence insights for where / when different estimators should be used. Despite this, for deep NNs, neither these variances nor their bounds can be computed efficiently due to the huge dimensionality of θ.This work focuses on estimating the diagonal entries of the FIM and their associated variances. Our results -including estimators of the FIM, their variances, and their variance bounds -can be implemented through automatic differentiation. These computational tools empower us to practically explore the trade-offs between the two estimators. For example, Fig. 1 shows natural gradient descent [1] for generalized linear models on a toy dataset, where Î2 (θ) is preferable (especially for regression) and Î1 (θ) suffers from high variance and an unstable learning curve. Our analytical results reveal how moments of the output exponential family and gradients of the NN in Eq. (1) affects the FIM estimators. We discover a general decomposition of the estimators' variances corresponding to the samples of x and y. We investigate different scenarios where each FIM estimator is the preferred one and then connect our analysis to the empirical FIM. 
related work: 
 Prior efforts aim to analyze the structure of the FIM of NNs with random weights [34, 14, 15, 3, 31] . This body of work hinges on utilizing tools from random matrix theory and spectral analysis, characterizing the behavior and statistics of the FIM. One insight is that randomly weighted NNs have FIMs with a majority of eigenvalues close to zero; with the other eigenvalues taking large values [14, 15] . In our work, the randomness stems from sampling from data distributions p(x, y) -which follows the principle of Monte Carlo (MC) information geometry [29] that approximates information geometric quantities via MC estimation. We examine a different subject on how the distribution of the FIM on a matrix manifold is affected by finite sampling of the data distribution.In the literature of NN optimization, a main focus is on deriving a computationally friendly proxy for the FIM. One can consider the unit-wise FIM [30, 20, 39, 3] (also known as quasi-diagonal FIM [30] ), where a block-diagonal approximation of the FIM is taken to capture intra-neuron curvature information. Or one can consider the block-diagonal layer-wise FIM where each block corresponds to parameters within a layer [19, 27, 32, 26, 12, 35, 13] . NN optimizers can approximate the inverse FIM [36] or approximate the product of the inverse FIM and the gradient vector [35] .Much less attention is paid to how related approximations deviate from the true FIM [11, 37] or how optimization is affected by such deviation [41] . For the univariate case, one can study the asymptotic variance of the Fisher information [11] with the central limit theorem. In deep NNs, the estimation variance of the FIM can be derived in closed form and bounded [37] . However, our former analysis [37] has two limitations: (1) the variance tensors are 4D and can not be easily computed;(2) only the norm of these tensors are bounded, and it is not clear how the variance is distributed among individual parameters. The current work tackles these limitations by focusing on the diagonal elements of the FIM. Our results can be computed numerically at a reasonable cost in typical learning settings. We provide novel bounds so that one can quantify the accuracy of the FIM computation w.r.t. individual parameters or subgroup of parameters.Issues of utilizing the empirical FIM to approximate the FIM have been highlighted [32, 25] . For example, estimators of the FIM do not in general capture any second-order information about the log-likelihood [18] . The empirical FIM is a biased estimator and can be connected with our unbiased estimators via a generalized definition of the Fisher matrix in Section 6.Alternative to the FIM, the Generalized Gauss-Newton (GGN) matrix -a Hessian approximatorwas originally motivated through the squared loss for non-linear models [25] . The GGN is equivalent to the FIM when a loss function is taken to be the empirical expectation of the negative log-likelihood of Eq. ( 1) [12, 32, 25] .Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: Limitations of our results are expressed through out the paper. See for instance Section 3 which discusses computation implications.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
variance of diagonal fim estimators: 
 In our notations, all vectors such as x, y, and θ are column vectors. We use k to index random samples x and y and use i and j to index the NN weights and biases θ. We shorthand h = p(y, x; θ) whenever the parameters θ is clear from context. To be consistent, we use ' | x conditioning' to distinguish between jointly calculated values versus conditioned values with fixed x. By default, the derivatives are w.r.t. θ. For example,∂ i h . = ∂h/∂θ i and ∂ 2 i h . = ∂ 2 h/∂θ 2 i .We adopt Einstein notation to express tensor summations, so that an index appearing as both a subscript and a superscript in the same term indicates a summation. For example, x a y a denotes a x a y a . For clarity, we mix standard Σ-sum and Einstein notation. We denote the variance and covariance of random variables by Var(•) and Cov(•), respectively.Based on the parametric form of the model in Eq. ( 1), the diagonal entries of the FIM estimators in Eq. ( 3) can be written as1 :Î1 (θ i ) . = Î1 (θ) ii = 1 N N k=1 ∂F (h(x k )) ∂θ i - ∂h a (x k ) ∂θ i • t a (y k ) 2 ; Î2 (θ i ) . = Î2 (θ) ii = 1 N N k=1 ∂ 2 F (h(x k )) ∂ 2 θ i - ∂ 2 h a (x k ) ∂ 2 θ i • t a (y k ) .Correspondingly, the i'th diagonal entry of the FIM I(θ), which is the expected value of Î1 (θ i ) and Î2 (θ i ), is denoted as I(θ i ). Notation is abused in I(θ i ), Î1 (θ i ), and Î2 (θ i ) as they depend on the whole θ vector rather than solely on θ i . Clearly Î1 (θ i ) ≥ 0, while there is no guarantee for Î2 (θ i ) which can be negative. Our results will be expressed in terms of the (central) moments of t(y):η a (x) . = E p(y | x) [t a (y)]; I(h | x) . = E p(y | x) [(t(y) -η(x))(t(y) -η(x)) ⊤ ]; K p (t | x) . = E p(y | x) [(t(y) -η(x)) ⊗ (t(y) -η(x)) ⊗ (t(y) -η(x)) ⊗ (t(y) -η(x))] ,where "⊗" denotes the tensor product. We denote the covariance of t w.r. Regression (Iso.) Gaussian ℜ T y 1 † 3 † Classification Categorical [C] ⊂ ℜ ( y = 0 , . . . , y = C ) min {σmax(x), 2 • min {σmax(x), 1 -∥σ(x)∥ 2 2 1 -∥σ(x)∥ 2 2In practice, both estimators can be computed via automatic differentiation [33, 6] . In terms of complexity, by restricting to just the diagonal elements I(θ i ), we need to calculate O(dim(θ)) elements (originally O(dim(θ) × dim(θ)) for the full FIM). Although the log-partition function for general exponential family distributions can be complicated, for the ones used in NNs (determined by the loss functions used in optimization) [37] the log-partition function F is usually in closed-form; and thus the cumulants η(x) and I(h | x) can be calculated efficiently.Indeed, the primary cost of the estimators comes from evaluating the gradient information of the NN, given by ∂h(x)/∂θ i and ∂ 2 h(x)/∂ 2 θ i . The former can be calculated easily. The latter is costly even when restricted to the diagonal elements of the FIM. With the Hessian's quadratic complexity, in practice approximations are used to reduce the computational overhead [4, 45, 46, 8] . In this case, additional error and (potentially) variance may be introduced as a result of the Hessian approximation. Note, the computational cost of the Hessian can still be manageable for the last few layers close to the output. By the chain rule, we only require a sub-computational graph from the output layer to a certain layer to compute the Hessian of that layer. Despite this, there is still a memory cost that scales quadratically with the number of parameters for non-linear activation functions [6] .The high cost of Hessian computation does not justify refraining from using Î2 . Depending on the setting (chosen loss function), an estimator's variance can outweigh the benefits of lower computational costs [37] . This is especially true when the FIM is used in an offline setting -where the Hessian's cost can be tolerated -to study, e.g., the singular structure of the neuromanifold [2, 40] , the curvature of the loss [7] , to quantify model sensitivity [28] , and to evaluate the quality of the local optimum [14, 15] , etc.To study the quality of Î1 (θ) and Î2 (θ), it is natural to examine the variance of the estimators [37] :V j (θ i | x) . = Var( Îj (θ i | x)), where Îj (θ i | x) . = Îj (θ | x) ii (j ∈ {1, 2}) is the i'th diagonal element of Îj (θ | x). Similar to Î1 (θ i ) and Î2 (θ i ), V j (θ i | x)and Îj (θ i | x) depend on the vector θ and are abuses of notation. An estimator with a smaller variance indicates that it is more accurate and more likely to be close to the true FIM. Based on the variance, one can derive sample complexity bounds of the diagonal FIM via Chebyshev's inequality, see for instance [37, Section 3.4] .By its definition, V j (θ i | x) has a simple closed form, which was proved in [37] and is restated below. Lemma 3.1. ∀x ∈ ℜ I , ∀i = 1, . . . , dim(θ),EQUATIONEQUATIONEQUATION) Given a fixed x ∈ ℜ I , both V 1 (θ i | x) and V 2 (θ i | x) have an order of O(1/N ), with N denoting the number of samples of y k . They further depend on two factors: ① the derivatives of the parameteroutput mapping θ → h stored in a T × dim(θ) matrix, either ∂ i h a (x) or ∂ 2 i h a (x), where the latter can be expensive to calculate; and ② the central moments of t(y), whose computation only scales with T (the number of output units) and is independent to dim(θ).From an information geometry [1] perspective, I(θ), V 1 (θ), and V 2 (θ) are all pullback tensors of different orders. For example, I(θ) is the pullback tensor of I(h) and the singular semi-Riemannian metric [40] . They induce the geometric structures of the neuromanifold (parameterized by θ) based on the corresponding low dimensional structures of the exponential family (parameterized by h).Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?• Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
practical variance estimation: 
 To further understand the dependencies of the derivative and central moment terms, the FIM I(θ i | x) and variances of estimators Îj (θ i | x) can be bounded to strengthen intuition and to provide a computationally convenient proxy of the interested quantities. Theorem 4.1. ∀x ∈ ℜ I ,EQUATIONEQUATIONEQUATIONwhere M = K p (t | x) -I(h | x) ⊗ I(h |EQUATIONTo help ground Theorem 4.1, we summarize different sufficient statistics quantities for common learning settings in Table 1 -with further learning setting implications presented in Section 5. Note that Eqs. ( 8) and ( 9) (and many subsequent results) can be further generalized for off-diagonal elements. See Appendix C for details. Compared to prior work [37] , Theorem 4.1 provides bounds for individual elements of the variance tensors, where the NN weights (the derivatives) and sufficient statistics (the eigenvalues) are neatly disentangled into a product. From a technical point of view, this comes from a difference in proof technique: we utilize variational definitions and computations of eigenvalues to establish bounds whereas [37] primarily applies Hölder's inequality.We stress that λmin (T ) and λmax (T ) in Eq. ( 10) correspond to tensor eigenvalues iff T is a supersymmetric tensor [23] (a.k.a. totally symmetric tensor), i.e., indices are permutation invariant. In this case, Eq. ( 10) is exactly the maximum and minimum Z-eigenvalues. These variational forms mirror the Courant-Fischer min-max theorem for symmetric matrices [42] . In the case of Eq. ( 8), with M = K p (t | x) -I(h | x) ⊗ I(h | x), the tensor is not a supersymmetric tensor in general. Despite this, we note that the lower bound of Eq. ( 8) is non-trivial. A weaker bound than Eq. ( 8) can be established based on the Z-eigenvalue of the supersymmetric tensorEQUATIONThe tensor eigenvalue is typically expensive to calculate. However in our case, the eigenvalues λmin (K p (t | x)) and λmax (K p (t | x)) on the RHS of Eqs. ( 11) and ( 12) can be calculated via [17] (11) .Examining Theorem 4.1 reveals several trade-offs. An immediate observation is that the first order gradients of h(x) correspond to the robustness of h to parameter misspecification (w.r.t. an input x). As such, from the bounds in Eqs. (7) and (8) , the scale of I(θ i | x) and V 1 (θ i | x) will be large when small shifts in parameter space yield large changes in the output h(x). Another observation is how the spectrum of I(h | x) affects the scale of I(θ i | x) and the estimator variances. In particular, when λ min (I(h | x)) increases, the scale of V 1 (θ i | x) decreases but the scale of I(θ i | x) and V 2 (θ i | x) increases. When λ max (I(h | x)) decreases, then the opposite scaling occurs. With these two observations, there is a tension in how the scale ofI(θ i | x) follows the different variances V 1 (θ i | x) and V 2 (θ i | x). The element-wise FIM I(θ i | x) follows V 1 (θ i | x) interms of the scale of NN derivatives ∥∂ i h(x)∥ 2 ; at the same time, I(θ i | x) follows V 2 (θ i | x) in terms of the spectrum of sufficient statistics moment I(h | x). Remark 4.4. Typically, h is the linear output units: h(x) = W -1 h -1 (x), where W -1 is the weights of the last layer, and h -1 (x) is the second last layer's output. We haveV 2 (θ i | x) = 0 ≤ V 1 (θ i | x) for any θ i in W -1 . A smaller variance V 2 (θ i | x) is guaranteed for the last layer regardless of the choice of the exponential family in Eq. ( 1). Remark 4.5. h(x) = w j -1 ϕ(h ⊤ -2 (x)w j -2 + C -2 ) + c -1 defines the NN mapping w.r.t. the j'th neuron in the second last layer, where w j -2 and w j -1 are incoming and outgoing links of the interested neuron, respectively; h -2 (x) is the output of the third last layer; and ϕ is the activation function. The 'constants' C -2 and c -1 denote an aggregation of all terms which are independent of w j -2 and w j -1 in their respective layers. The Hessian ofh k (x) w.r.t. w j -2 is ∂ 2 h k (x) = (w j -1 ) k • ϕ ′′ (h ⊤ -2 (x)w j -2 + C -2 ) • (h -2 (x)h ⊤ -2 (x)). By Theorem 4.1, V 2 (θ i | x) can be arbitrarily small depending on ϕ ′′ (h ⊤ -2 (x)w j -2 + C -2 ). For example, if ϕ(t) = 1/(1 + exp(-t)), then ϕ ′′ (t) = ϕ(t)(1 -ϕ(t))(1 -2ϕ(t)).In this case, for a neuron in the second last layer, a sufficient condition for V 2 (θ i | x) = 0 (and having Î2 favored against Î1 ) is h ⊤ -2 (x)w j -2 + C -2 = 0 for the neuron's pre-activation. When the pre-activation value is saturated (-∞ or ∞), we also have thatV 1 (θ i | x) = V 2 (θ i | x) = 0. Alternatively, suppose that ϕ(t) = SoftPlus(t) . = log(1 + exp(t)), a continuous relaxation of ReLU, then ϕ ′′ (t) = ϕ ′ (t)(1 -ϕ ′ (t)) where ϕ ′ (t) = 1/(1 + exp(-t)). Then a sufficient condition for V 2 (θ i | x) = 0 with V 1 (θ i | x) ̸ = 0 for a neuron in the second last layer is h ⊤ -2 (x)w j -2 + C -2 → +∞.These observations are further clarified by looking at related quantities over multiple parameters. So far we have only examined the variance of the FIM element-wise w.r.t. parameters θ i . To study all parameters θ jointly, we consider the trace variances of the FIM estimators: for any j ∈ {1, 2}, V j (θ | x) denotes the trace of the covariance matrix of diag( Îj (θ | x)), where diag(•) extracts a matrix's diagonal elements into a column vector. We present upper bounds of these joint quantities. Corollary 4.6. For any x ∈ ℜ I ,EQUATIONEQUATIONwhere dHes(h| x) . = (diag(Hes(h 1 | x)), . . . , diag(Hes(h T | x))) and ∥ • ∥ F is the Frobenius norm.This upper bound comes from integrating the parameter-wise variances in Theorem 4.1 and incorporating a trace variance bound which utilizes the full spectrum of the NN derivatives and sufficient statistics quantities. This is fully depicted in Theorem D.1. Lower bounds can also be derived in terms of singular values (deferred to the Appendix). Note the upper bounds in Corollary 4.6 can be improved by expressing the min function's first term with singular value quantities.Having the min function in Corollary 4.6 is helpful as it shows a trade-off between two upper bounds: the scale of NN derivatives ∂h(x) and ∂ 2 h(x) versus the spectrum of the sufficient statistic terms.In the case of Eqs. ( 13) and ( 15), the trace of I(h | x) is exactly the sum of all eigenvalues, including λ max (I(h | x)). This can be helpful when the scale of the NN derivatives are not bounded by a small value. It should be noted that, by the chain rule, these NN derivatives scale with the overall sharpness / flatness [22] of the landscape of the loss, i.e., the log-likelihood of Eq. ( 1). For NNs with large derivatives, the first term of the min could yield tight bounds of the variance, and one can therefore avoid dealing with the quadratic scaling of ∥∂h(x)∥ in the second term. On the other hand, if the sharpness of the NN h can be controlled, e.g. via sharpness aware minimization [10] , then one can benefit from the second term of the min and avoid computing the full spectrum of I(h | x) in the first term.Joint FIM Estimators In the above, we considered the variance of conditional FIMs, which can scale differently depending on the input x. Prior work's analysis was limited to that of conditional FIMs (and their estimators) [37] . Nevertheless, the 'joint FIM' estimators Îj (θ i ) depend on sampling of x w.r.t. the data distribution q(x). The bounds in Eq. ( 7) can be extended to the joint FIM I(θ i ) .= Eq(x) I(θ i | x) by simply taking an expectation Eq(x) over the bounds. To analyze the variances of the joint FIM estimators V j (θ i ), we present the following theorem which connects the prior results established for V j (θ i | x), e.g. Theorem 4.1, via the law of total variance. Theorem 4.7. For any j ∈ {1, 2}, given N x samples of x ∼ q(x) and N samples ofEQUATIONwhereVar (I(θ i | x)) is the variance of I(θ i | x) w.r.t. q(x).The dependence on N , the number of samples of y for each fixed x, is hidden in V j (θ i | x). When N = 1, the hierarchical sampling described in the Theorem corresponds to an i.i.d. sampling of the joint distribution p(x, y).The variance incurred when estimating the FIM has two components. The first term on the RHS of Eq. ( 16) characterizes the randomness of the FIM w.r.t. q(x), i.e., the input randomness. It vanishes when the FIM is estimated by taking the expectation w.r.t. q(x), or the number of samples x is large enough. The second term (although also depending on N x ) comes from the sampling of y | x according to p(y | x), i.e., the output randomness, which scales with the central moments of t(y).If the NN is trained so that p(y | x) tends to be deterministic, this term will disappear leaving the first term to dominate. Eq. ( 16) can be further generalized using the law of total covariance to extend prior work considering conditional FIM covariances [37] to joint FIM covariances. Theorem 4.7 connects the variance of assuming a fixed input x with multiple samples y k with the variance of pairs of samples (x k , y k ). The V j (θ i | x) bounds in this section can thus be applied to the corresponding joint variance V j (θ i ) by using this theorem. This is straightforward and omitted.The first variance term in Theorem 4.7 is difficult to compute in practice: it relies on how the closed-form FIM varies w.r.t. q(x). As such, it is useful to bound the first term into computable quantities. Lemma 4.8.Var (I(θ i | x)) ≤ Eq(x) ∥∂ i h(x)∥ 4 2 • λ 2 max (I(h | x)) .This upper bound is very similar to the 4th central moment term λmax (K (I(h | x)) ̸ = λmax (K p (t | x)). This is especially true for the classification and regression problems explored in this paper (see Table 1 ). However, the maximum eigenvalues can be related for exponential families with bounded sufficient statistic via Proposition 4.3, making both bounds depend only on λ max (I(h | x)).The total number of samples of (x k , y k ) is N x • N . In terms of sample complexity of the joint variance, using Theorem 4.7 and Lemma 4.8, the bound's rate is given by O(1/N x + 1/(N x • N )). 
case studies: 
 To make our theoretic results more concrete, we consider regression and classification settings, which correspond to specifying the exponential family in Eq. ( 1) to an isotropic Gaussian distribution and a categorical distribution, respectively. We include an empirical analysis of NNs trained on MNIST. Notably, our analysis considers general multi-dimensional NN output. This extends the case studies of [37] which was limited to 1D distributions due to the limitations of their bounds (and their associated computational costs of dealing with a 4D tensor of the full covariance).To characterize regression, we consider Gaussian distributions. As per Eq. ( 1), we have t(y) = y ∈ ℜ T and base measure π(y) ∝ exp(-1 2 y ⊤ y).This corresponds to the case where h(x) is trained via the squared loss. In this case,I(h | x) = I and K p abcd (t | x) = I ab (h | x) • I cd (h | x) + I ac (h | x) • I bd (h | x) + I ad (h | x) • I bc (h | x).The derivatives of the log-partition function F (h) yields these central moments [37, Lemma 5] . To apply Theorem 4.1, we examine the extreme eigenvalues of I(h | x) and K p (t | x) -I(h | x) ⊗ I(h | x). Proposition 5.1. Suppose that Eq. ( 1) is an isotropic Gaussian distribution. Then:λ min (I(h | x)) = λ max (I(h | x)) = 1; λmin (K p (t | x) -I(h | x) ⊗ I(h | x)) = λmax (K p (t | x) -I(h | x) ⊗ I(h | x)) = 2.Hence, for regression the eigenvalues of sufficient statistics quantities in our bounds are equal. As such, in this case, the bound forI(θ i | x), V 1 (θ i | x), and V 2 (θ i | x) in Theorem 4.1 are all equalities.As I(θ i | x) can be written exactly in terms of the gradients of h, in practice one does not need to utilize a random estimator when computing the conditional FIM, which simplifies to a Gauss-Newton matrix for the squared loss [25] . When computing the FIM over a random sample x, a variance still appears due to Theorem 4.7. By Lemma 4.8 and Proposition 5.1, the variance over a joint distribution is bounded by a function of the derivative ∂ i h(x) over the marginal input distributions:V(θ i ) ≤ Eq(x)[∥∂ih(x)∥ 4 2 ] ≤ max x∈Supp(q) ∥∂ i h(x)∥ 4 2 .In other words, the overall variance to approximate the joint FIM is bounded by the gradients of the NN outputs.Classification: Categorical Distribution For multi-class classification, we instantiate our exponential family with a categorical distribution (with π(y) = 1). This corresponds to training a classifier NN with the log-loss. Let Y = [C] for T = C classes defining the possible labels. Let t(y) = ( y = 1 , . . . y = C ), where r = 1 when the predicate r is true and r = 0 otherwise. Noting our results do not depend on minimal sufficiency, this t(y) is sufficient but not minimal sufficient. In this setting, the NN outputs h correspond to the logits of the label probabilities. The resulting probabilities p(y | x) are the softmax values of h denoted by σ(x). = SoftMax(h(x)) ∈ [0, 1] T .Under this setting, we have I(h | x) = Diag(σ(x))-σ(x)σ(x) ⊤ (where Diag(σ(x)) is the diagonal matrix with its diagonal entries set to σ(x)), whose eigenvalues do not follow a convenient pattern as C increases [44] . Likewise, the maximum eigenvalue of K p (t | x) is not available in simple closed form. As such, we provide upper bounds for the maximum eigenvalues of I(h | x λ max (I(h | x)) ≤ m(x); and λmax (K p (t | x) -I(h | x) ⊗ I(h | x)) ≤ 2 • m(x), where m(x) . = min σ max (x), 1 -∥σ(x)∥ 2 2 .This upper bounds provides a tension. When the first term of m(x) is maximized, the second is minimized, and vice-versa. In particular, the dominating term depends on the uncertainty of the NN's output. When the NN's output is near random, e.g. at initialization, the first term will dominate with σ max (x) ≈ 1/C. However, as the NN becomes more certain with its prediction, the second term will start dominating: a more deterministic output p(y| x) → 1 implies that λ max (I(h | x)) → 0.Empirical Verification: Classification We examine the MNIST classification task [21] (CC BY-SA 3.0) using multilayer perceptrons (MLP) with four densely connected layers, sigmoid activations, and a dropout layer. For classification, we consider a categorical distribution with C = 10 class labels. For a random x from the test set, we compute both estimators Î1 (θ i | x) and Î2 (θ i | x) using N = 5, 000 samples. We record the variances of each estimator and compute their bounds based on Theorem 4.1. For all 20 training epochs, the Fisher information (FI) and their variances of individual parameters are aggregated via arithmetic averages over four parameter groups (corresponding to the four layers).In Fig. 2 , we present the variance scale of the estimators Î1 (θ i | x) and Î2 (θ i | x) in log-space; and the tightness of the bounds in Theorem 4.1 by consider the log-ratio log UB VI(θi | x) , where UB is the upper bounds in Theorem 4.1. In this experiment, the UB is much tighter than the lower bound (LB), which is omitted in the figures for clarity. More experimental results are given in Appendix F.We varied the NN's architecture and activation function. Across different settings, the proposed UB and LB are always valid. In Fig. 2 , one can observe that the diagonal FIM and the associated The closer to 0, the better the UB. In the right most column, the variance of Î2 vanishes:V 2 (θ i | x) = 0 ≤ V 1 (θ i | x).Thus related curves of Î2 are not shown.variances have a small magnitude. For example, in the first layer,V 1 (θ i | x) and V 2 (θ i | x)are roughly e -10 ≈ 5 × 10 -5 . The log-ratio log UB V1(θi | x) ≈ 4 means that the UB is roughly 50 times larger thanV 1 (θ i | x). Comparatively, V 2 (θ i | x)has a tighter UB which is approximately 10 times larger than itself. The UB serves as a useful hint on the order of magnitude of the variances. In Appendix D, we present tighter bounds which are more expensive to compute.In the first three layers of the MLP, V 1 (θ i | x) presents a smaller value than V 2 (θ i | x), meaning that Î1 can more accurately estimate the diagonal FIM. Interestingly, this is not true for the last layer: V 2 (θ i | x) becomes zero while Î1 presents the largest variance across all parameter groups. Due to this, one should always prefer Î2 over Î1 for the last layer. In the last two layers, Î2 is in simple closed form and, hence, does not need automatic differentiation to calculate (see Remarks 4.4 and 4.5). The shape of the variance curves are sensitive to the choices of activation functions ϕ and inputs x. In general, the variance in the first few epochs presents more dynamics than the rest of the training process. If one uses log-sigmoid activations ϕ(t) =log(1 + exp(-t)) (which is equivalent to ϕ(t) = -SoftPlus(-t), as per Remark 4.4), the variances of Î1 and Î2 only appear in the randomly initialized NN and quickly vanish once training starts, as shown in Appendix F. In this case, the learner more easily approaches a nearly linear region of the loss landscape where local optima lie. In practice, one should estimate and examine the scale of variances -which should not be neglected as per Fig. 2 -before choosing a preferred diagonal FIM estimator. 
relationship with the "empirical fisher": 
 In some scenarios, even the estimators of the diagonal FIM Î1 (θ) and Î2 (θ) can be prohibitively expensive. Part of the cost comes from requiring label samples y k for each x k , as per Eq. (3). For example, when the FIM is used in an iterative optimization procedure, y k 's need to be re-sampled at each learning step w.r.t. the current h alongside their backpropagation (accounting for sampling).As such, alternative 'FIM-like' objects have been explored which replace the samples from p(y | x) with samples from an underlying true (but unknown) data distribution q(y | x) [20, 27] . We define the data's joint distribution as q(x, y) . = q(x)q(y | x). Analogous to the FIM, the data Fisher information matrix (DFIM) can be defined as the PSD tensor I(θ) . = Eq(x)[I(θ | x)], withEQUATIONwhereI(h | x) = Eq(ŷ | x) (t(ŷ) -η(x))(t(ŷ) -η(x))⊤ denotes the 2nd (non-central) moment of (t(ŷ)η(x)) w.r.t. q(ŷ | x), and ∂h/∂θ is the Jacobian of the map θ → h. In the special case that q(y | x) = p(y | x; θ), then I(θ | x) becomes exactly I(θ | x).The DFIM I(θ | x) in Eq. ( 17) is a more general definition. Compared to the FIM I(θ | x), it yields a different PSD tensor on the θ parameter space (the neuromanifold) depending on a dis-tribution q(x, y), which is neither necessarily on the same neuromanifold nor necessarily parametric at all. The asymmetry in the true data distribution and the empirical one results in different geometric structures [5] . By definition, we have I(θ | x) ⪰ (∂KL/∂θ) (∂KL/∂θ) ⊤ , where.= q( ŷ | x) log q( ŷ | x) p( ŷ | x;θ) d ŷ is the Kullback-Leibler (KL) divergence, or the loss in a parameter learning scenario. The DFIM can be regarded as a surrogate function of the squared gradient of the KL divergence. It is a symmetric covariant tensor and satisfies the same rule w.r.t. reparameterization as the FIM. Consider the reparameterization θ → ζ, the DFIM becomesI(ζ | x) = (∂θ/∂ζ) ⊤ I(θ | x) (∂θ/∂ζ).Notice that η(x) .= Eq(ŷ | x) [t(ŷ)] ̸ = η(x) in general. As such, there will be a miss-match when utilizing I(h | x) as a substitute for I(h | x). However, as learning progresses and p(ŷ | x) becomes more similar to the data's true labeling posterior q(ŷ | x), the DFIM will become closer to the FIM.If q(x, y) = 1 N N k=1 δ(x -x k ) • δ(y -y k ) isdefined by the observed samples, DFIM gives the widely used "Empirical Fisher" [25] , whose diagonal entries areÎ(θ i ) = 1 N N k=1 (∂h a i (x k ) • (t a (ŷ k ) -η a (x k ))) 2 ,where (x 1 , ŷ1 ), . . . , (x N , ŷN ) are i.i.d. sampled from q(x, ŷ). Similar to Î1 (θ i | x), an estimator with a fixed input x can be considered, denoted as Î(θ i | x).Given the computational benefits of using the data directly -bypassing a separate sampling routine -many popular optimization methods employ the empirical Fisher or its approximation. For instance, the Adam optimizer [16] uses the empirical Fisher to approximate the diagonal FIM. However, switching from sampling y k to ŷk is anything but superficial [25, Chapter 11] -Î(θ) is not an unbiased estimator of I(θ) asI(h | x) is different from I(h | x).The biased nature of the empirical Fisher affects the other moments as well. In particular, we do not have the same equivalence of covariance and the metric being pulled back by θ → h [38] . Lemma 6.1. Given the conditional data distribution q(ŷ | x), the covariance of t given x is given byEQUATIONwhere ∆H(x) = (η(x)η(x))(η(x)η(x)) ⊤ .As a result, although the variance of the estimator Î(θ i | x) takes a similar form to V 1 (θ i | x) (i.e., Eq. ( 8)), its sufficient statistic terms do not exclusively consist of central moments. Noting the miss-match in η(x) ̸ = η(x), Lemma 6.1 reveals an additional term which shifts I(h | x) away from the 2nd central moment of t(ŷ) (w.r.t. q(ŷ | x)). Instead, these sufficient statistic terms correspond to non-central moments of t(ŷ)η(x). Some corresponding empirical Fisher / DFIM bounds are characterized in Appendix G. 
conclusion: 
 We have analyzed two different estimators Î1 (θ) and Î2 (θ) for the diagonal entries of the FIM. The variances of these estimators are determined by both the non-linearly of the neural network and the moments of the exponential family. We have identified distinct scenarios on which estimator is preferable. For example, ReLU networks can only apply Î1 (θ) due to a lack of smoothness. As another example, Î2 (θ) has zero variance in the last layer and thus is always preferable than Î1 (θ).Similarly, in the second last layer, Î2 (θ) has a simple closed form and potentially preferable for neurons in their linear regions (see Remark 4.5). In general, one has to apply Theorem 4.1 based on their specific neural network and settings and choose the estimator with the smaller variance.Our results suggest that, from a variance perspective, uniformly utilizing one of the FIM estimators Îj (θ) is often suboptimal in NNs. Our work has further extended from analyzing the conditional FIM estimators Îj (θ | x) to the joint FIM estimators Îj (θ); and we have examined the relationship between the investigated estimators and the empirical Fisher. Future directions include extending the analysis of the variance of FIM estimators to block diagonals (e.g. [26, 35] ) and adapting current NN optimizers (e.g. [16] ) to incorporate the variance of FIM estimators. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: No data or model release.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [Yes] Justification: License and credit are given for MNIST. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. 
 
 Appendix: The following section describes the data and models of Fig. 1 . In general, the toy data and models constructed consists of taking 1D output setting presented by Section 5, where the NN h θ (x) is a linear function.The 2D input data x ∈ ℜ 2 is sampled from a simple isotropic centered Gaussian N (0, I). A linear response variable a ∈ ℜ is defined by the following:a = θ ⊤ true x;where θ true = (1, 1).The outputs of y for the cases of regression and classification are differentiated by how a is used in sampling:y regression ∼ N (µ = 1, σ = 1) y classification ∼ Bern(p = σ(a)),where σ(z) = (1 + exp(-z)) is the logistic function.The model h θ (x) = θ ⊤ x consists of a linear function; and the exponential family Eq. ( 1) is chosen to be a 1D isotropic Gaussian and binary multinomial distribution (Bernoulli) for regression and classification, respectively. This corresponds to Section 5 for 1D outputs. Notice that the model exactly matches the data generating function.Natural gradient descent (NGD) is taken using both Î1 (θ) and Î2 (θ). The estimated FIM utilize only a single y | x sample for each input x. We use a learning rate of η = 0.01 over 256 epochs. A training set of 256 data points are sampled. At each iteration of NGD, we sample 4 random points from the training set for the update. The test loss is evaluated on a test set of 4096 data points sampled.Larger version of Fig. 1 with additional variance sum plotted over time is given by Fig. I. Note that variance sum is including off diagonals. Further note that the variance is calculated over joint sample in (x, y).We further present other random seed of the teaser plot in Figs. II to IV.We consider the diagonal entries of the conditional FIM I(θ i | x) and the conditional variances V j (θ i | x) of its estimators in closed form.Proof of Lemma 3.1. The proof directly follows from [37, Equation 6], [37, Theorem 4] , and [37, Theorem 6] . In what follows, we provide a proof of the Lemma utilizing the notation of this paper for completeness. We prove the statement one equation at a time. For Eq. ( 4), we consider the following computation.I(θ i | x) = E p(y | x;θ) ∂ log p(y | x; θ) ∂θ ∂ log p(y | x; θ) ∂θ ⊤ = E p(y | x;θ) ∂ t ⊤ (y)h θ (x) -F (h θ (x)) ∂θ ∂ t ⊤ (y)h θ (x) -F (h θ (x)) ∂θ ⊤ = E p(y | x;θ)   ∂h θ (x) ∂θ ⊤ t(y) - ∂F (h) ∂h h=h θ (x) t(y) - ∂F (h) ∂h h=h θ (x) ⊤ ∂h θ (x) ∂θ ⊤   = E p(y | x;θ) ∂h θ (x) ∂θ ⊤ (t(y) -η(x)) (t(y) -η(x)) ⊤ ∂h θ (x) ∂θ ⊤ = ∂h θ (x) ∂θ ⊤ E p(y | x;θ) (t(y) -η(x)) (t(y) -η(x)) ⊤ ∂h θ (x) ∂θ ⊤ = ∂h θ (x) ∂θ ⊤ I(h | x) ∂h θ (x) ∂θ ⊤ .Using Einstein notation and restricting the partial derivative to a component of θ yields the desired result.For Eq. ( 5), we shorthand δ(y) = t(y)η(x). Note that the Î1 (θ i | x) estimator can be written as follows: Thus, we haveÎ1 (θ i | x) = 1 N N k=1 ∂F (h(x)) ∂θ i - ∂h a (x) ∂θ i • t a (y k ) 2 = 1 N N k=1 ∂h a (x) ∂θ i • η a (x) - ∂h a (x) ∂θ i • t a (y k ) 2 = 1 N N k=1 ∂h a (x)∂h b (x)δ a (y k )δ b (y k ).V 1 (θ i | x) = Var Î1 (θ i | x) = Var 1 N N k=1 ∂h a (x)∂h b (x)δ a (y k )δ b (y k ) = 1 N • Var ∂h a (x)∂h b (x)δ a (y k )δ b (y k ) = 1 N • E p(y | x;θ) ∂h a (x)∂h b (x)δ a (y)δ b (y) 2 - E p(y | x;θ) ∂h a (x)∂h b (x)δ a (y)δ b (y) 2 .Let us compute each of these terms.E p(y | x;θ) ∂h a (x)∂h b (x)δ a (y)δ b (y) 2 = E p(y | x;θ) ∂h a (x)∂h b (x)∂h c (x)∂h d (x)δ a (y)δ b (y)δ c (y)δ d (y) = ∂h a (x)∂h b (x)∂h c (x)∂h d (x) E p(y | x;θ) [δ a (y)δ b (y)δ c (y)δ d (y)] = ∂h a (x)∂h b (x)∂h c (x)∂h d (x)K p abcd (t | x). And, E p(y | x;θ) ∂h a (x)∂h b (x)δ a (y)δ b (y) 2 = ∂h a (x)∂h b (x) E p(y | x;θ) [δ a (y)δ b (y)] 2 = ∂h a (x)∂h b (x)I ab (h | x) 2 = ∂h a (x)∂h b (x)∂h c (x)∂h d (x)I ab (h | x)I cd (h | x) = ∂h a (x)∂h b (x)∂h c (x)∂h d (x) (I(h | x) ⊗ I(h | x)) abcdSimplifying all term yields the result as required.Finally, for Eq. ( 6) we consider the following simplification of the estimator.i | x) = 1 N N k=1 ∂ 2 F (h(x)) ∂ 2 θ i - ∂ 2 h a (x) ∂ 2 θ i • t a (y k ) = 1 N N k=1 ∂ ∂θ i ∂h a (x) ∂θ i • η a (x) - ∂ 2 h a (x) ∂ 2 θ i • t a (y k ) = 1 N n k=1 ∂h a (x) ∂θ i • ∂η a (x) ∂θ i + ∂ 2 h a (x) ∂ 2 θ i • η a (x) - ∂ 2 h a (x) ∂ 2 θ i • t a (y k ) = 1 N n k=1 ∂ i h a (x) • ∂η a (x) ∂θ i -∂ 2 i h a (x) • δ a (y k ) = 1 N n k=1 ∂ i h a (x) • ∂ i h b (x) • I ab (h | x) -∂ 2 i h a (x) • δ a (y k ) = ∂ i h a (x) • ∂ i h b (x) • I ab (h | x) - 1 N n k=1 ∂ 2 i h a (x) • δ a (y k ) ,where the last line follows from [37, Lemma 2] (a result of p(y | x; θ) following an exponential family, see [1] ).Notice that the first quantity is a constant w.r.t. the randomness of y k . As such, we can simplify the variance calculation as follows.V 2 (θ i | x) = Var Î2 (θ i | x) = Var ∂ i h a (x) • ∂ i h b (x) • I ab (h | x) - 1 N n k=1 ∂ 2 i h a (x) • δ a (y k ) = Var 1 N n k=1 ∂ 2 i h a (x) • δ a (y k ) = 1 N Var ∂ 2 i h a (x) • δ a (y) = 1 N • E p(y | x;θ) ∂ 2 i h a (x) • δ a (y) 2 - E p(y | x;θ) ∂ 2 i h a (x) • δ a (y) 2 = 1 N • E p(y | x;θ) ∂ 2 i h a (x) • ∂ 2 i h b (x) • δ a (y) • δ b (y) -∂ 2 i h a (x) • E p(y | x;θ) [δ a (y)] 2 = 1 N • ∂ 2 i h a (x) • ∂ 2 i h b (x) • E p(y | x;θ) [δ a (y) • δ b (y)] = 1 N • ∂ 2 i h a (x) • ∂ 2 i h b (x) • I ab (h | x),where the second last line follows from the fact that η(x) = Ep(y | x;θ) [t(y)] and thus Ep(y | x;θ) [δ(y)] = 0. This yields the desired result.Lemma 3.1 shows that, for the former, V 1 (θ i | x) only depends on 1st order derivatives; while V 2 (θ i | x) only depends on the 2nd order derivatives. For the latter, V 1 (θ i | x) depends on both the 2nd and 4th central moments of t(y); while V 2 (θ i | x) only depends on the 2nd central moments.Given I ab (h | x) and ∂ i h a (x), the computational complexity of all diagonal entries I(θ i | x) is O(T 2 dim(θ)). If K p abcd (t | x) and ∂ 2 i h a (x) are given, then the computational complexity of the variances in Eqs. ( 5) and ( 6) is respectively O(T 4 dim(θ)) and O(T 2 dim(θ)). Each requires to evaluate a T × dim(θ) matrix, either ∂ i h a (x) or ∂ 2 i h a (x) -which can be expensive to calculate for the latter. This is why we need efficient estimators and / or bounds for the tensors on the LHS of Eqs. ( 4) to (6) .We consider an off-diagonal version of the bound given by Theorem 4.1. Notice that in terms of the dependence on neural network weights, the only change is splitting the "responsibility" of the i'th and j'th parameter norms.EQUATION)EQUATIONwhereγmax (M) = sup u:∥u∥2=1,v:∥v∥2=1 u a v b u c v d M abcd γ max (M ) = sup u:∥u∥2=1,v:∥v∥2=1 u a v b M ab .Proof. The proof follows similarly to Appendices J and K, where the primary difference is just swapping the regular eigenvalue-like quantities with the γ variational forms.It should be noted that the corresponding lower bounds become trivial as the additional degree of freedom of having an inf over both u and v causes the corresponding γ min definition to have negative quantities. Although it is unclear what the "tensor-like" variational quantity γmax (M) will be, for a matrix, we have the following equivalence. Lemma C.2. γ max (A) = s max (A), where s max (A) is the maximum singular value of A.Proof. The proof follows from optimizing over u and v separately:γ max (A) = sup u:∥u∥2=1 sup v:∥v∥2=1 u a v b A ab = sup u:∥u∥2=1 sup v:∥v∥2=1 u ⊤ Av = sup v:∥v∥2=1 (Av) ⊤ Av ∥Av∥ 2 = sup v:∥v∥2=1 v T (A ⊤ A)v.This is equivalent to the square root of the maximal eigenvalue of A ⊤ A, which is exactly the maximum singular value.Hence for the Î2 we have the following.EQUATIOND Bounding the Trace Variance by Full Spectrum Theorem D.1. For any x ∈ ℜ I , T t=1EQUATIONEQUATION)EQUATION)where s 2 i (A) = λ i (A ⊤ A) denotes the i-th singular values, M is the "reshaped" matrix of M defined in Theorem 4.1i.e. there exists j, k such that M jk = M abcd for all a, b, c, d,dHes(h | x) = (diag(Hes(h 1 | x)), . . . , diag(Hes(h T | x))), and vJac(h | x) = (vec(∂ 1 h(x)∂ 1 h(x) ⊤ ), . . . , vec(∂ T h(x)∂ T h(x) ⊤ )).Proof. The proof follows from a generalized Ruhe's trace inequality [24] :Theorem D.2. For A, B ∈ ℜ n×n Hermitian matrices, we have thatn i=1 λ i (A) • λ n-i+1 (B) ≤ tr (AB) ≤ n i=1 λ i (A) • λ i (B).We prove the result for each equations.For readability, we let J ia = ∂ i h a (x).For Eq. ( 22):One can notice that the trace of the FIM can exactly be expressed as the trace of two matrices.tr (I(θ | x)) = dim(θ) i=1 ∂ i h a (x)∂ i h b (x)I ab (h | x) = I ab (h | x) dim(θ) i=1 J ia J ib = I ab (h | x) dim(θ) i=1 (J ⊤ ) ai J ib = I ab (h | x)(J ⊤ J) ab = tr (J ⊤ J)I(h | x) .Thus, noting that the eigenvalue of the "squared" matrix is the matrix's singular value λ t (J ⊤ J) = s 2 t (J), with Theorem D.2, we have that:T t=1 s 2 t (∂h(x)) • λ T -t+1 (I(h | x)) ≤ tr (I(θ | x)) ≤ T t=1 s 2 t (∂h(x)) • λ t (I(h | x)) .For Eq. ( 23):Noting that M abcd = K p abcd (t | x) -I ab (h | x) • I cd (h | x). Furthermore, we have that vJac(h | x) = (vec(∂ 1 h(x)∂ 1 h(x) ⊤ ), . . . , vec(∂ T h(x)∂ T h(x) ⊤ )).Let us define the following 3D tensor withJ iab = ∂ i h a (x)∂ i h b (x) = (∂ i h(x)∂ ⊤ i h(x)) ab . V 1 (θ | x) = 1 N dim(θ) i=1 ∂ i h a (x)∂ i h b (x)∂ i h c (x)∂ i h d (x)M abcd = 1 N M abcd dim(θ) i=1 J iab J icd = 1 N T a,b=1 T c,d=1 M abcd dim(θ) i=1 J iab J icd = 1 N T 2 j=1 T 2 k=1 M jk dim(θ) i=1 vJac ij (h | x)vJac ik (h | x) = 1 N T 2 j=1 T 2 k=1 M jk (vJac ⊤ (h | x)vJac(h | x)) jk = 1 N tr M(vJac ⊤ (h | x)vJac(h | x)) .Thus, again simplifying the eigenvalue of the "squared" matrix, with Theorem D.2, we have that:1 N T t=1 s 2 t (vJac(h | x)) • λ T -t+1 M ≤ tr Î1 (θ | x) ≤ 1 N T t=1 s 2 t (vJac(h | x)) • λ t M .For Eq. ( 24):Similar to Eq. ( 22), we only need to rearrange the summation. Notice that dHes(h | x) = (diag(Hes(h 1 | x)), . . . , diag(Hes(h T | x))),thus dHes ia (h | x) = ∂ 2 i (h a | x). V 2 (θ | x) = 1 N dim(θ) i=1 ∂ 2 i h a (x)∂ 2 i h b (x)I ab (h | x) = 1 N I ab (h | x) dim(θ) i=1 ∂ 2 i h a (x)∂ 2 i h b (x) = 1 N I ab (h | x) dim(θ) i=1 dHes ia (h | x)dHes ib (h | x) = 1 N I ab (h | x) dim(θ) i=1 (dHes ⊤ ) ai (h | x)dHes ib (h | x) = 1 N I ab (h | x)(dHes ⊤ (h | x)dHes(h | x)) ab = 1 N tr I(h | x)(dHes ⊤ (h | x)dHes(h | x)) .Thus, again simplifying the eigenvalue of the "squared" matrix, with Theorem D.2, we have that:1 N T t=1 s 2 t (dHes(h | x)) • λ T -t+1 (I(h | x)) ≤ tr Î2 (θ | x) ≤ 1 N T t=1 s 2 t (dHes(h | x)) • λ t (I(h | x)) .Proof. We first notice that the exponential family density is given by, p(y | x) = exp(h y (x) -F (h(x))) and thus also haveF (h(x)) = log T t=1 exp(h t (x))The first order derivative follows as, ∂F (h)∂h i h=h(x) = exp(h i (x)) T t=1 exp(h t (x)) = σ i (h(x))As such, the second order derivatives also follow,∂ 2 F (h) ∂h i ∂h j h=h(x) = exp(h i (x))δ ij • T t=1 exp(h t (x)) -exp(h i (x)) exp(h j (x)) T t=1 exp(h t (x)) 2 = σ i (h(x)) • δ ij -σ i (h(x))σ j (h(x)).As such, we have thatI(h | x) = Diag(σ(x)) -σ(x)σ(x) ⊤ .In the following section we present additional details and results for the experimental verification we conduct in Section 5.We note that to calculate the diagonal Hessians required for the bounds and empirical FIM calculations, we utilize the BackPACK [6] for PyTorch. Additionally, to calculate the sufficient statistics moment's spectrum, we explicitly solve the minimum and maximum eigenvalues via their optimization problems. For 2D tensors / matrices, we utilize numpy.linalg.eig. For 4D tensors, we utilize PyTorch Minimize [9] , a wrapper for SciPy's optimize function.We present Figs. V to VIII which are the exact same experiment run in Section 5, but with different initial NN weights and random inputs. G "Empirical Fisher" Continued Noting Lemma 6.1's characterization of the covariance, we are able to characterize the variance of the diagonal elements of Î(θ | x), denoted as where K(h | x) the 4th (non-central) moment of (t(ŷ)η(x)) w.r.t. q(ŷ | x).V(θ i | x) . = Var( Î(θ i | x)). Corollary G.1. For any x ∈ ℜ I , V(θi | x) = 1 N ∂ih a (x)∂ih b (x)∂ih c (x)∂ih d (x) (K abcd (t | x) -I ab (h | x) ⊗ I cd (h | x)) = 1 N ∂ih a (x)∂ih b (x)∂ih c (x)∂ih d (x)K abcd (t | x) - 1 N ∂ih ⊤ (x) (Cov q (t | x) + ∆H(x)) ∂ih(x)As a result of the similarity of the functional forms of the empirical Fisher Î(θ) and the FIM estimator Î1 (θ), it is not surprising that Corollary G.1 is similar to the variance of Î1 (θ i | x). Indeed, applying Lemma 6.1 will give the exact same functional form with the 2nd central moments of t(y) w.r.t. p(y | x) exchanged with 2nd non-central moments of (t(ŷ)η(x)) w.r.t. q(ŷ | x). V(θ i | x) is therefore determined by the 2nd and the 4th moment of (t(ŷ)η(x)) up to the parameter transformation θ → h. Subsequently, the bounds presented for V 1 (θ i | x) (Eq. ( 8) and Corollary 4.6) can be similarly adapted for V(θ i | x).The extension of V(θ i | x) to V(θ i ) can also be proven in a similar manner to Theorem 4.7.Corollary G.2. Given N x samples of x ∼ q(x) and N samples ofEQUATIONwhereVar (I(θ i | x)) is the variance of I(θ i | x) w.r.t. q(x). If q(x, y) = 1 N N k=1 δ(x -x k ) • δ(y -ŷk ) for a set of observations {(x k , ŷk )} N k=1, then one can directly evaluate the DFIM without sampling and achieve zero variance, i.e., Î(θ) = I(θ). In this scenario, there is a clear trade-off between the estimators of the FIM in Eq. ( 3) and the DFIM. The estimators of the FIM are unbiased, but have a variance; while the DFIM has zero variance, but is a biased approximation of the FIM.  1 √  1 UB of √  1 √  2 UB of √  2In what follows, we derive the alternative equations for Î1 (θ i ) and Î2 (θ i ) presented in Section 3.That is, we seek to derive the following equations:EQUATIONEQUATIONWe calculate the equations separately. H.I Eq. ( 26)Proof. For Eq. ( 26), we note that∂ log p(y k | x k ) ∂θ i = ∂ ∂θ i t ⊤ (y k )h(x k ) -F (h(x k )) = t a (y k ) ∂h a (x k ) ∂θ i -F ′ a (h(x k )) ∂h a (x k ) ∂θ i = t a (y k ) ∂h a (x k ) ∂θ i -η a (x k ) ∂h a (x k ) ∂θ i = (t a (y k ) -η a (x k )) • ∂h a (x k ) ∂θ i ,where we note that F ′ a (h(x k )) = η a (x k ) which follows from the connection to expected parameters and partition functions of exponential families, see e.g. [37] .Then Eq. ( 26) follows immediately.H.II Eq. ( 27)Proof. For Eq. ( 27), we also calculate the derivative:∂ log p(y k | x k ) ∂θ i = ∂ ∂θ i t ⊤ (y k )h(x k ) -F (h(x k )) = t a (y k ) • ∂h a (x k ) ∂θ i - ∂F (h(x k )) ∂θ i . Then ∂ 2 log p(y k | x k ) ∂ 2 θ i = ∂ ∂θ i t a (y k ) • ∂h a (x k ) ∂θ i - ∂F (h(x k )) ∂θ i = t a (y k ) • ∂ 2 h a (x k ) ∂ 2 θ i - ∂ 2 F (h(x k )) ∂ 2 θ i .Now we can prove the bounds of Eq. ( 8)Proof. From Lemma 3.1, we have thatV 1 (θ i | x) = 1 N ∂ i h a (x)∂ i h b (x)∂ i h c (x)∂ i h d (x) [K abcd -I ab (h | x) • I cd (h | x)] ,where we shorthand K abcd = K p abcd (t | x) We bound two terms.∥∂ i h(x)∥ 4 2 • λmin (K) ≤ ∂ i h a (x)∂ i h b (x)∂ i h c (x)∂ i h d (x)K abcd ≤ ∥∂ i h(x)∥ 4 2• λmax (K), which follows directly from Lemma J.1We now bound the second term in a similar way, taking v a .= ∂ i h a (x) and noting thatv a v b v c v d I ab (h | x)I cd (h | x) = (v a v b I ab (h | x)) 2which directly gives us,∥v∥ 2 2 • λ min (I(h | x)) 2 ≤ v a v b v c v d I ab (h | x)I cd (h | x) ≤ ∥v∥ 2 2 • λ max (I(h | x)) 2 .which follows from Lemma I.1.Thus, together these bounds prove Eq. ( 8).K Proof of Eq. ( 9)From Lemma 3.1 we have that,V i 2 = 1 N ∂ 2 i h a (x)∂ 2 i h b (x)I ab (h L ).∥∂ 2 i h(x)∥ 2 2 • λ min (I(h)) ≤ ∂ 2 i h b (x) • ∂ 2 i h b (x) • I ab (h) ≤ ∥∂ 2 i h(x)∥ 2 2• λ max (I(h)) , which follows from Lemma I.1. This immediately gives the bound as required.Proof. The corollary holds from distributing the inf or sup and examining how the variational definition of the generalized 'eigenvalue' simplifies under tensor products. Indeed, for the minimum case,λmin (K p (t | x) -I(h | x) ⊗ I(h | x)) = inf u:∥u∥2=1 u a u b u c u d (K p abcd (t | x) -I ab (h | x) • I cd (h | x)) ≥ inf u:∥u∥2=1 u a u b u c u d K p abcd (t | x) + inf u:∥u∥2=1 u a u b u c u d (-I ab (h | x) • I cd (h | x)) = inf u:∥u∥2=1 u a u b u c u d K p abcd (t | x) - sup u:∥u∥2=1 u a u b u c u d (I ab (h | x) • I cd (h | x)) = inf u:∥u∥2=1 u a u b u c u d K p abcd (t | x) - sup u:∥u∥2=1 u a u b I ab (h | x) 2 ≥ inf u:∥u∥2=1 u a u b u c u d K p abcd (t | x) - sup u:∥u∥2=1 u a u b I ab (h | x) 2 ,where the last line holds from the fact that I ab (h | x) is PSD (thus the inner Einstein summation is always positive).Taking definitions of the types of eigenvalues, gives the statement.We note that the 'max' case follows identically.Additionally, for the lower bound, we can show the non-triviallity of the non-negativity of the minimum eigenvalue.We note that K p abcd (t | x) = Ep [vavbvcvd] , where v = t(y)η(x). Thus we have thatλmin (K p (t | x) -I(h | x) ⊗ I(h | x)) = inf u:∥u∥2=1 u a u b u c u d (K p abcd (t | x) -I ab (h | x) • I cd (h | x)) = inf u:∥u∥2=1 E p u a u b u c u d (v a v b v c v d -I ab (h | x) • I cd (h | x)) = inf u:∥u∥2=1 E p u a u b (v a v b -I ab (h | x)) 2 ≥ 0.Equality holds from simply looking at the definition of K p (t | x) and I(h | x) (as moments).M Proof of Proposition 4.3Proof. Letting v = t(y)η(y), we note that the maximum eigenvalue is given by,λmax (K p (t | x)) = sup u:∥u∥2=1 u a u b u c u d K p abcd (t | x) = sup u:∥u∥2=1 E p u a u b u c u d v a v b v c v d = sup u:∥u∥2=1 E p (u ⊤ v) 4 = sup u:∥u∥2=1 E p (u ⊤ v) 2 (u ⊤ v) 2 ≤ sup u:∥u∥2=1 E p (∥u∥ 2 • ∥v∥ 2 ) 2 (u ⊤ v) 2 ≤ B • sup u:∥u∥2=1 E p (u ⊤ v) 2 = B • λ max (I(h | x)) .N Proof of Corollary 4.6Proof. We split up the proof into the two arguments of the various min-function.For the right term:Suppose that we have a bound such that V j (θ i | x) ≤ αβ i . Then,∥V j (θ | x)∥ 2 2 = dim(θ) i=1 (V j (θ i | x)) 2 ≤ dim(θ) i=1 (αβ i ) 2 = α 2 dim(θ) i=1 β 2 i . ∥V j (θ | x)∥ 2 ≤ α dim(θ) i=1 β 2 i .Taking the appropriate α and β from Eqs. ( 8) and ( 9) proves the case for Eqs. ( 13) and (15) .For Eq. ( 14), that is takingα = 1 N • λmax (K p (t | x)) -λ 2 min (I(h | x)) ; β i = ∥∂ i h(x)∥ 4 2 .Where we note thatdim(θ) i=1 (∥∂ i h(x)∥ 4 2 ) 2 = dim(θ) i=1   T t=1 [∂ i h t (x)] 2 2   2 = dim(θ) i=1 T t=1 [∂ i h t (x)] 2 4 ≤   dim(θ) i=1 T t=1 [∂ i h t (x)] 2   4 =   dim(θ) i=1 T t=1 [∂ i h t (x)] 2   2 = ∥∂h(x)∥ 4 F .For Eq. (15) , that is takingα = 1 N • λ max (I(h | x));β i = ∥∂ 2 i h(x)∥ 2 2 .Where we note thatdim(θ) i=1 (∥∂ 2 i h(x)∥ 2 2 ) 2 = dim(θ) i=1 T t=1 [∂ 2 i h t (x)] 2 2 ≤   dim(θ) i=1 T t=1 [∂ 2 i h t (x)] 2   2 = dim(θ) i=1 T t=1 ∂ 2 i h t (x) 2 = ∥dHes(h | x)∥ 2 F .For the left term:We take the largest singular value of the network derivative term. We then further notice that s max (A) ≤ ∥A∥ F from norm ordering (of the matrix 2-norm).To further elaborate on the Eq. ( 14) case, we further need to simplify the following:s max (vJac) ≤ ∥vJac(h | x)∥ F = dim(θ) i=1 ∥∂ i h(x)∂ i h ⊤ (x)∥ 2 F = T a,b=1 dim(θ) i=1 (∂ i h a (x)) 2 (∂ i h b (x)) 2 ≤ T a,b=1 ∥(∂h a (x)) 2 ∥ 2 • ∥(∂h b (x)) 2 ∥ 2 = T a=1 ∥(∂h a (x)) 2 ∥ 2 2 = T a=1 ∥(∂h a (x)) 2 ∥ 2 = T a=1 dim(θ) i=1 (∂ i h a (x)) 4 ≤ T a=1 dim(θ) i=1 |(∂ i h a (x)) 2 | = ∥∂h(x)∥ 2 F ,where the last inequality follows from the norm ordering∥ • ∥ 2 ≤ ∥ • ∥ 1 .O Proof of Theorem 4.7To prove the Theorem, we will utilize the law of total variances. We note, that by the premise of the Theorem, we are sampling N x many samples from q(x) and N many samples from q(y | y) for each y initially sampled. To make this clear, the samples and sampling will be notated by:x k ∼ q(x) y l | x k ∼ p(y | x k )Note that using these samples, our empirical estimators for the FIM (for either estimator) will be of the form:Î1 (θ i ) = 1 N x x k   1 N y l | x k f (x k , y l | x k )   ,for an appropriately chosen f . This also gives:Îj (θ i | x) = 1 N y l | x f (x, y l | x ).Now, we simplify the variance as follows:Var   1 N x x k   1 N y l | x k f (x k , y l | x k )     = 1 N 2 x x k   Var   1 N y l | x k f (x k , y l | x k )     = 1 N 2 x x k Var Îj (θ i | x k ) = 1 N 2 x x k Var x k E y 1 | x k ,...,y N | x k Îj (θ i | x k ) + E x k Var y 1 | x k ,...,y N | x k Îj (θ i | x k ) = 1 N 2 x x k Var x k [I(θ i | x k )] + E x k [V j (θ i | x k )] = 1 N x Var x [I(θ i | x)] + E x [V j (θ i | x)] .As required.For V 1 (θ) Proof. V 1 (θ i ) = 1 N E p(x,y) ∂ log p(y | x) ∂θ i 2 -E p(x,y) ∂ log p(y | x) ∂θ i 2 .Let δ a (x, y) . = (t(y)η(x)).E p(x,y) ∂ log p(y | x) ∂θ i 2 = E p(x,y) ∂h a (x) ∂θ i ∂h b (x) ∂θ i ∂h c (x) ∂θ i ∂h d (x) ∂θ i δ a (x, y)δ b (x, y)δ c (x, y)δ d (x, y) = E q(x) ∂h a (x) ∂θ i ∂h b (x) ∂θ i ∂h c (x) ∂θ i ∂h d (x) ∂θ i E p(y | x) [δ a (x, y)δ b (x, y)δ c (x, y)δ d (x, y)] = E q(x) ∂h a (x) ∂θ i ∂h b (x) ∂θ i ∂h c (x) ∂θ i ∂h d (x) ∂θ i K p abcd (t | x)And:E p(x,y) ∂ log p(y | x) ∂θ i 2 = E p(x,y) ∂h a (x) ∂θ i ∂h b (x) ∂θ i δ a (x, y)δ b (x, y) 2 = E q(x) ∂h a (x) ∂θ i ∂h b (x) ∂θ i E p(y | x) [δ a (x, y)δ b (x, y)] 2 = E q(x) ∂h a (x) ∂θ i ∂h b (x) ∂θ i I ab (h | x) 2 = E q(x) ∂h a (x) ∂θ i ∂h b (x) ∂θ i I ab (h | x) 2 -Var X ∂h(x) ∂θ i ⊤ I(h | x) ∂h(x) ∂θ i = E q(x) ∂h a (x) ∂θ i ∂h b (x) ∂θ i I ab (h | x) 2 -Var X (I(θ i | x)) . Together: V 1 (θ i ) = 1 N Var X (I(θ i | x)) + 1 N E q(x) ∂h a (x) ∂θ i ∂h b (x) ∂θ i ∂h c (x) ∂θ i ∂h d (x) ∂θ i [K p abcd (t | x) -I ab (h | x) • I cd (h | x)] For V 2 (θ) Proof. V 2 (θ i ) = 1 N Var (η a (x) -t a (y)) ∂ 2 h a (x) ∂θ i ∂θ i + I(θ i | x) = 1 N      Var (η a (x) -t a (y)) ∂ 2 h a (x) ∂θ i ∂θ i (a) +Var (I(θ i | x)) +2 Cov (η a (x) -t a (y)) ∂ 2 h a (x) ∂θ i ∂θ i , I(θ i | x) (b)      . (a) = Var (η a (x) -t a (y)) ∂ 2 h a (x) ∂θ i ∂θ i = E p(x,y) (η a (x) -t a (y)) ∂ 2 h a (x) ∂θ i ∂θ i 2 -E p(x,y) (η a (x) -t a (y)) ∂ 2 h a (x) ∂θ i ∂θ i 2 = E p(x,y) (η a (x) -t a (y)) ∂ 2 h a (x) ∂θ i ∂θ i 2 -E q(x) ∂ 2 h a (x) ∂θ i ∂θ i E p(y | x) [(η a (x) -t a (y))] 2 = E p(x,y) (η a (x) -t a (y)) ∂ 2 h a (x) ∂θ i ∂θ i 2 -0 = E q(x) ∂ 2 h(x) ∂θ i ∂θ i ⊤ E p(y | x) (η(x) -t(y)) (η(x) -t(y)) ⊤ ∂ 2 h(x) ∂θ i ∂θ i = E q(x) ∂ 2 h(x) ∂θ i ∂θ i ⊤ I(h | x) ∂ 2 h(x) ∂θ i ∂θ i . (b) = Cov (η a (x) -t a (y)) ∂ 2 h a (x) ∂θ i ∂θ i , I(θ i | x) = E p(x,y) (η a (x) -t a (y)) ∂ 2 h a (x) ∂θ i ∂θ i I(θ i | x) -E p(x,y) (η a (x) -t a (y)) ∂ 2 h a (x) ∂θ i ∂θ i E p(x,y) [I(θ i | x)]= 0, which follows by taking the 'partial' expectation (y | x) for both terms.Thus together,V 2 (θ i ) = 1 N Var (I(θ i | x)) + 1 N E q(x) ∂ 2 h(x) ∂θ i ∂θ i ⊤ I(h | x) ∂ 2 h(x) ∂θ i ∂θ i .P Proof of Lemma 4.8Proof. The lower bound holds from just considering the non-negativity of variance. For the upper bound, we utilize the bound directly consider the bounds of Eq. ( 7),Var (I(θ i | x)) = E q(x) I(θ i | x) 2 -E q(x) [I(θ i | x)] 2 ≤ E q(x) I(θ i | x) 2 ≤ E q(x) ∥∂ i h(x)∥ 4 2 • λ 2 max (I(h | x)) .We first derive the statistics I(h | x) and K p (t | x) presented in "Regression: Isotropic Gaussian Distribution" Section 5. It follows that from the regression setting, we have that,F (h(x)) = log π(y) • exp(t ⊤ (y)h(x)) = log π(y) • exp(y ⊤ h(x)),where notably, by definition, π(y) is independent of learned parameter h(x).As such, we have that:∂ ∂h i F (h) h=h(x) = 1 π(y) • exp(t ⊤ (y)h(x)) • π(y) • exp(t ⊤ (y)h(x)) • h i (x) = E p(y | x) [y i ].Now we note that Ep(y | x)[y i ] is exactly h i (x) as the parameter h(x) specifies the mean of the (isotropic) multivariate normal distribution. As such we have that,∂ ∂h i F (h) h=h(x) = h(x) I(h | x) = ∂ 2 ∂h∂h ⊤ F (h) h=h(x) = I.Furthermore, by [37, Lemma 5], we have that,K p abcd (t | x) = ∂ 4 F (h) ∂h a ∂h b ∂h c ∂h d h=h(x) + I ab (h | x) • I cd (h | x) + I ac (h | x) • I bd (h | x) + I ad (h | x) • I bc (h | x) = 0 + I ab (h | x) • I cd (h | x) + I ac (h | x) • I bd (h | x) + I ad (h | x) • I bc (h | x) = I ab (h | x) • I cd (h | x) + I ac (h | x) • I bd (h | x) + I ad (h | x) • I bc (h | x).In summary, we have,K p abcd (t | x) = I ab (h | x) • I cd (h | x) + I ac (h | x) • I bd (h | x) + I ad (h | x) • I bc (h | x) (K p (t | x) -I(h | x) ⊗ I(h | x)) abcd = I ac (h | x) • I bd (h | x) + I ad (h | x) • I bc (h | x).Proof. The minimum and maximum eigenvalues of I(h | x) follows directly noting that the trace of a matrix is the sum of eigenvalues. As such, from the statistics presented above we have that the minimum and eigenvalue must be 1.K p (t | x) -I(h | x) ⊗ I(h | x) = I ac (h | x) • I bd (h | x) + I ad (h | x) • I bc (h | x)follows from the variational definition Eq. (10) . For instance, for the minimum eigenvalue,inf u:∥u∥2=1 u a u b u c u d (I ac (h | x) • I bd (h | x) + I ad (h | x) • I bc (h | x)) = 2 • inf u:∥u∥2=1 ∥u∥ 2 2 = 2.The maximum eigenvalue is proven identically.We first prove the following corollary which connects the maximum eigenvalues of K(t | x) to the maximum eigenvalues of I(h | x).Corollary R.1. Suppose that the exponential family in Eq. ( 1) is specified by a categorical distribution. Then,EQUATIONProof. As we are consider a categorical distribution we have thatv i = 1 -σ i (h) if y = i 0 -σ i (h) if y ̸ = i .Thus we have that |v i | ≤ 1. Furthermore, note that the maximum ℓ 2 -norm that we can have is ∥v∥ 2 ≤ √ 2. Note that this is tight when the positive and negative mass are placed only two distinct coordinates, i.e., (-1, 1, . . .).Thus using Proposition 4.3, the result follows. Now by using Corollaries 4.2 and R.1, the remainder of the proof, all we require is the bounding of λ max (I(h | x)).Proof. The first term in the maximum eigenvalue follows from,λ max (I(h | x)) = λ max Diag(σ(x)) -σ(x)σ(x) ⊤ ≤ λ max (Diag(σ(x))) -λ min σ(x)σ(x) ⊤ = max k σ k (x).The second term in the maximum follows from the trace of I(h | x) being the sum of total eigenvalues.Proof. The proof follows from the standard definition of covariance. Denoting η(x) . = Eq(y | x) [t(y)], we have:Cov q (t | x) = E q(y | x)t(y)t ⊤ (y)η(x) η⊤ (x).I(h | x) = E q(y | x) (t(y) -η(x))(t(y) -η(x)) ⊤ = E q(y | x) t(y)t ⊤ (y) -η(x) η⊤ (x) -η⊤ (x)η ⊤ (x) + η(x)η ⊤ (x).Cov q (t | x) = I(h | x) + η(x) η⊤ (x) + η⊤ (x)η ⊤ (x) -η(x)η ⊤ (x) -η(x) η⊤ (x) = I(h | x) -(η(x) -η(x))(η(x) -η(x)) ⊤ .As required.Proof. We calculate the variance:V(θ i ) = 1 N E q(y | x) ∂ log p(y | x) ∂θ i 2 -E q(x,y) ∂ log q(y | x) ∂θ i 2 .Each of the terms can be calculated: Let δ a (x, y) . = (t(y)η(x)).E p(x,y) ∂ log p(y | x) ∂θ i 2 = E q(y | x) ∂h a (x) ∂θ i ∂h b (x) ∂θ i ∂h c (x) ∂θ i ∂h d (x) ∂θ i δ a (x, y)δ b (x, y)δ c (x, y)δ d (x, y) = ∂h a (x) ∂θ i ∂h b (x) ∂θ i ∂h c (x) ∂θ i ∂h d (x) ∂θ i E q(y | x) [δ a (x, y)δ b (x, y)δ c (x, y)δ d (x, y)] = ∂h a (x) ∂θ i ∂h b (x) ∂θ i ∂h c (x) ∂θ i ∂h d (x) ∂θ i K abcd (t | x).And:E p(y | x) ∂ log p(y | x) ∂θ i 2 = E p(y | x) ∂h a (x) ∂θ i ∂h b (x) ∂θ i δ a (x, y)δ b (x, y) 2 = ∂h a (x) ∂θ i ∂h b (x) ∂θ i E p(y | x) [δ a (x, y)δ b (x, y)] 2 = ∂h a (x) ∂θ i ∂h b (x) ∂θ i I ab (h | x) 2 .Together with Lemma 6.1 proves the theorem.Proof. The proof follows identically to that of Theorem 4.7 with densities changed.NeurIPS Paper Checklist Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.This and subsequent derivations can be found in the appendix., 


Decision: All reviewers lean towards accept with at least one reviewer advocating for a strong accept. This is a highly theoretical paper that is likely to be of interest to others interested in theory. 

However, there are several weaknesses discussed by the reviewers that should be addressed in the final version. In particular, more discussion on the differences between this work and [37] should be included. Also, the authors should work to better "remark on the connections between the equations" in the technical sections to make these sections more accessible and also further discuss when one estimator should be preferred over the other.

Reviews: review: paper_summary: paper_summary: Summary
-------
The paper studies two estimators for estimating the diagonal of a Fisher information
matrix (FIM) of a parametric machine learning model. Both estimators are based on equivalent
expressions for the FIM. The first is the standard definition, expressed in terms of the
derivatives of the log likelihood, and the latter, expressed in terms of the second
derivative when the model is twice differentiable with respect to the parameters. The
estimators replace the expectation with an average over the data. The authors provide
bounds on the variance of both (unbiased) estimators in terms of the minimum/maximum
eigenvalues of the FIM, and then provide bounds on these eigenvalues under special cases.
They also provide a modest empirical evaluation.


Overall, while I think there is a fair amount of technical meat in the paper, I am 
unconvinced about the problem. I will wait to hear from more expert reviewers in this 
field before I make my decision. 


Detailed comments:
------------------

While I am theoretically inclined, I am not an expert in this sub-field. Therefore, I will
keep my comments at a high level.

Significance of problem: I struggled to understand when estimating the FIM will be useful.
- One point raised by authors is that it is used in optimizing NNs. If this was true, it
  would have been nice to see the estimators used in NN optimization empirically, or a
  theoretical analysis showing how the bounds derived here percolate to bounds on
  optimization.
- Another explanation given is that it could be used to analyze the structure fo NNs.
  But this is stated vaguely, and it wasn't clear to me that such an analysis would help
  design better neural networks.

Significance of techniques: The proof in the appendix is substantial but the authors have 
not stated which techniques used here are novel. Currently technical results are presented
one after the other. It might help the reader/reviewer
appreciate the contributions of the paper if the authors can highlight the novel tools 
they used to derive these bounds.

Presentation and writing: The technical content was clear for the most part. However,
the writing was very equation-driven. I understand that this is in part due to the nature 
of the paper, but the authors should simplify the exposition further. One option is to
cut down on some of the technical content.
strength_weakness: strength_weakness: See above
questions: questions: See above
limitations: limitations: See above
paper_summary: paper_summary: In this paper, the authors analyzed two different estimators, $I_1(\theta)$ and  $I_2(\theta)$, for the diagonal elements 
of the Fisher information matrix of a parametric neural net model. These diagonal elements are approximations for the entire matrix,
which is unfeasible to be calculated in real neural nets models. They identified the situations in which each of these estimators is preferable with respect to the other.
strength_weakness: strength_weakness: This is a highly theoretical paper and I have difficulty following all the findings presented by the authors. One of the main problems with this paper is that it has a large intersection with Soen and Sun (2021) and this becomes clear while we read the paper. In the Related work, Soen and Sun (2021) is mentioned in passing and it is a major reference work for the current manuscript. It is worth to have a better explanation of the relative merits of these two papers. Both papers analyzed two different estimators, $I_1(\theta)$ and  $I_2(\theta)$, for the Fisher information matrix of a parametric neural net model. While the older paper focuses on the matrix as a whole, the present manuscript focuses on the diagonal elements. However, the work of Soen and Sun (2021) is repeatedly referred to during the technical development of the paper. The contribution of this paper seems marginal when compared to the previous paper but still relevant.   

The results are derived assuming $\theta$ is equal to its true value. In practice, the maximum likelihood estimator (MLE) of 
$\theta is used, which introduces additional variability. The paper does not adequately address the implications of this assumption or the consequences when using the MLE, which limits the practical applicability of the results.

The paper does not discuss the impact of model misspecification on the results. In real-world scenarios, the data might not follow the proposed model exactly, and the inverse-Fisher information matrix might not represent the correct asymptotic variance of MLE estimators. The authors should comment on how their findings would be affected in such cases.
questions: questions: The results are all given using $\theta$ equal to its true value. In practice, we use the MLE estimator of $\theta$. The results 
are not valid anymore as we have this additional source of variability. The authors may comment on the consequences of this?
limitations: limitations: This is a theoretical paper with no immediate connection to societal impact.
paper_summary: paper_summary: The paper analyses two estimators of the Fisher Information matrix and specifically their variances, which (from reference 37) have closed-form but non-practical expressions. The authors show a sequence of inequalities and derive practical bounds for these variances, both element-wise and trace-wise, both conditioned to x and not. Several remarks and trade-offs on these bounds clarify when an estimator is preferable over the other.
The results hold in the general exponential family likelihood but are also further concretized for regression and classification.
strength_weakness: strength_weakness: I have found no major weaknesses in the paper. I report some minor things.

I personally dislike the notation $\hat{\mathcal{I}}_1(\theta_i)$ for the i-th element of the diagonal introduced in line 87. This notation suggest that this term only depends on $\theta_i$, but it actually depends on the whole $\theta$ vector. On the other hand, this is quite a light easy notations, and I can't think of alternatives which are not much heavier to read.

Also, I find the definition in line 120 "is defined by $\hat{\mathcal{I}}_j(\theta_i)$ with $x=x_1=...=x_N$" to be yes understandable but not super clean. The alternatives are not as compact, but I'd recommend using the extra space to clarify this a bit better. At least you should write that $y_1,...,y_N$ are i.i.d. from $p(y|x;\theta)$ and not from $p(y,x;\theta)$ as in $\hat{\mathcal{I}}_j(\theta_i)$.

It would be nice to explicitly write the proof for the variance closed expression in Lemma B.1, yet clarifying that it's not a paper contribution. This is for two reasons: (1) it would make the paper self contained, without the need for the reader to read reference 37 and (2) the notation is slightly different and it would be easier for the reader to be consistent.
questions: questions: In line 168. Shouldn't "of Eqs. (4) and (6)" instead be "of Eqs. (4) and (5)"? And consistently shouldn't "$\mathcal{V}_1(\theta_i|x)$" be "$\mathcal{V}_2(\theta_i|x)$"? \
I ask this because "small shifts in parameter space yield large changes in the output" to me refers to the magnitude of the network jacobian $||\partial_i h(x)||$, and that appear in the bound for $V_1$, not for $V_2$. Am I misunderstanding something or is this a typo?

In line 246. Isn't the bound sample complexity $\mathcal{O}(\frac{1}{N_x} + \frac{1}{N_y})$? I can't see why it should be the product, can you elaborate more on this derivation?
limitations: limitations: I think it would be nice to clarify better that the work is incremental on "Soen et all - On the variance of the Fisher information for deep learning" [reference 37] and their closed form expression reported in the appendix in Lemma B.1. (also it would be nice to move these closed form expressions to the main paper).
paper_summary: paper_summary: The authors study two popular estimators of the Fisher information matrix with respect to neural network parameters. They derive upper and lower bounds for the variance of these estimators and showcase them in applications to regression and classification problems.
strength_weakness: strength_weakness: - The paper is highly technical and the numerical examples are fairly toy.
- The paper takes strong inspiration from a previous work [Soen, Alexander, and Ke Sun. "On the variance of the Fisher information for deep learning." Advances in Neural Information Processing Systems 34 (2021): 5708-5719] on the topic. As such, I question the value the present submission adds to this line of research. In particular, I think the submission does a poor job in highlighting the central novel aspects of their work and reviewing the field.
questions: questions: - The authors should clarify clearly to what extent their submission improves upon previous works, particularly [Soen, Alexander, and Ke Sun. "On the variance of the Fisher information for deep learning." Advances in Neural Information Processing Systems 34 (2021): 5708-5719]. How do their bounds compare to the ones previously derived?

- In my opinion, the impact of the submission would be substantially improved if the authors state more clearly in what scenarios which of the two estimators will be preferred and provide a rough guideline for practitioners.

- Could the authors provide a numerical example for regression similar to Fig. 2?
limitations: limitations: There is no further need for authors to address limitations.

Comments: title: Response
Keep score
comment: Thanks for responding to my questions. A better clarification of the contribution of this paper wrt [37] is appreciated.
I keep my positive evaluation and my score.
Thanks for responding to my questions. Based on the reply, I have increased my score.
I thank the authors for their concise rebuttal and have updated my score accordingly.
Thanks for the clarification. While I do agree with the other reviewers that the improvements made on top of [37] could have been clarified better, I do think that they are overall clear and valuable. Hence, I keep my score and strongly support acceptance.