New insights and perspectives on the natural gradient method

1. Introduction and Overview

• Natural gradient descent replaces the standard gradient with the gradient times the inverse Fisher, often requiring fewer iterations but becoming impractical for very large models without approximations.Approximations make the Fisher easier to compute, store, and invert.Source: §1 · p. 3
• Natural gradient’s distribution-space interpretation yields invariance to smooth invertible reparameterizations, but finite parameter-space steps undermine this guarantee and its connection to objective progress.The invariance interpretation is exact only as the step size becomes arbitrarily small, and is only approximate for larger practical steps.Source: §1 · pp. 3–4
• The paper recasts natural gradient descent as a second-order method whose Fisher-based local quadratic model can be more informative and less conservative than gradient descent.The Fisher approximates the Hessian as an expected model-output Hessian and is often equivalent to the Generalized Gauss-Newton matrix.Source: §1 · p. 4
• The second-order perspective motivates trust-region-style damping and Tikhonov regularization, which are presented as crucial for fast, robust practical performance when quadratic approximations are inaccurate over long distances.These techniques have established uses in neural-network optimization and nonlinear regression.Source: §1 · p. 4
• The empirical Fisher is easier to implement but is substantially less useful as Fisher or curvature approximation, breaks parts of natural-gradient theory, and nevertheless preserves approximate parameterization invariance.The paper also examines diagonal uses in methods such as RMSprop and Adam.Source: §1 · pp. 4–5
• For convex quadratics, the paper derives sharper stochastic second-order convergence results that account for initialization: averaged stochastic gradient descent matches natural gradient asymptotically, while second-order methods can progress faster initially.With certain decaying learning rates, stochastic second-order methods can also achieve much better finite-iteration objective values than stochastic gradient descent; the theory does not explain their neural-network advantage over Hessian methods.Source: §1 · p. 5

• Feed-forward neural networks consist of input, hidden, and output layers, with each unit receiving the previous layer’s outputs and applying an affine map followed by a scalar nonlinearity.Source: §2 · p. 7
• Formally, the network maps x ∈ R^m0 and parameters θ ∈ R^n, comprising layer-specific weight matrices and biases, to output f(x, θ) = a_ℓ.Source: §2 · p. 7
• The layer activities satisfy a_0 = x, while each activation function φ_i produces the vector-valued nonlinear transformation at layer i, often coordinate-wise.Source: §2 · p. 7
• Most results extend beyond neural-network structure to any function f(x, θ) differentiable in both inputs and parameters.Source: §2 · p. 7

Objective

• Supervised learning minimizes an objective over training pairs (x, y) by choosing parameters θ whose network predictions match the targets.The training set consists of input-target pairs, and the network output is treated as a prediction.Source: §3 · p. 7
• The loss function L(y, z) measures the disagreement between a target y and prediction z.Source: §3 · p. 7

Prediction and loss

• Predictions may be compared directly with targets using losses such as squared error, or may parameterize a distribution scored by negative log probability.For discrete targets, f(x, θ) can provide multinomial probabilities, with the loss given by the negative log probability of the observed target.Source: §3 · p. 7

• Learning minimizes KL divergence from the target joint distribution to the learned distribution, equivalently the expected conditional KL divergence over target inputs.The input distribution is fixed, while the learned model represents the conditional distribution of outputs given inputs.Source: §4 · p. 8
• When only samples are available, replacing target distributions with empirical training distributions yields objectives based on finite input or input-output samples.If conditional targets are known, their KL divergences can be used; with one output sample per input, the formulation uses empirical conditional distributions.Source: §4 · p. 8
• With one observed output per input, the resulting objective is effectively the same as standard maximum-likelihood learning.The equivalence holds up to an additive constant that is technically infinite because of the Dirac-delta representation.Source: §4 · p. 8
• The learned conditional distribution can be expressed as a deterministic predictor composed with an output conditional distribution, and many losses can be represented this way.For losses whose exponentiated negative values have a finite, z-independent integral over outputs, the output distribution is proportional to exp(−L(y, z)).Source: §4 · p. 9

• The Fisher information matrix is positive semidefinite and equals the negative expected Hessian of log p(x, y|θ).Source: §5 · p. 9
• The natural gradient is defined by preconditioning the objective gradient with the Fisher information matrix.Source: §5 · p. 9
• Because q(x) is parameter-independent, the Fisher can be expressed as an expectation over Qx of the conditional model’s Fisher information.Source: §5 · p. 9
• When exact integration is difficult, Qx can be replaced by its empirical version ˆQx, as in Park et al. (2000).Analytical integration may be especially difficult for multilayer neural networks.Source: §5 · pp. 9–10
• For negative log-likelihood losses, the Fisher is an expected loss Hessian under Px,y, whereas the objective Hessian uses ˆQx,y; thus, F and H approximate each other.Source: §5 · p. 10

• The negative gradient is steepest descent under the Euclidean norm, so its direction depends strongly on the parameterization of θ.Source: §6 · p. 10
• The negative natural gradient is steepest descent in distribution space when local distance is measured by KL divergence.Source: §6 · pp. 10–11
• The Fisher matrix F supplies the local quadratic approximation to KL divergence, translating distribution-space geometry into parameter-space updates.This follows from the KL Taylor expansion, whose remainder contains terms of order O(d^3) and higher.Source: §6 · p. 11
• Although F and the natural gradient depend on the basis in θ-space, the induced direction in distribution space is invariant to parameterization, unlike the ordinary gradient generally.Source: §6 · p. 11
• A smoothly varying positive semidefinite F defines a Riemannian manifold whose metric is locally equivalent to KL divergence, though KL itself is not a valid symmetric metric.Source: §6 · pp. 10–11

• Second-order methods choose an update by minimizing a local quadratic model, solving a curvature-weighted linear system when the curvature matrix is positive definite.The update is then scaled by a step factor before parameters are changed.Source: §7 · pp. 11–12
• Gradient descent is the special case B_k = βI, using one curvature estimate for every direction; β = L safely ensures convergence for convex, Lipschitz-smooth objectives but can move unnecessarily slowly.The model upper-bounds the objective when β = L, preventing predicted decreases from masking sharp increases, while overestimating curvature in low-curvature directions.Source: §7 · p. 12
• By rescaling gradients according to curvature eigenvalues, second-order methods address parameter coupling and directional scale variation that make neural-network training sensitive to gradient-descent learning rates.These issues include strong local dependencies and large scale differences associated with phenomena such as vanishing gradients.Source: §7 · p. 12
• Newton’s method uses the Hessian for the most accurate local quadratic model and can rapidly explore low-curvature directions, but indefinite Hessians and excessive trust in the model can produce unbounded or harmful updates.Gauss–Newton approximations and damping or trust-region techniques can largely mitigate these problems.Source: §7 · pp. 12–13
• Neural-network parameter dimensionality makes forming, storing, and inverting the full curvature matrix impractical, motivating low-rank or diagonal approximations and incomplete model optimization.Approximate Newton methods address the n^2-entry matrix bottleneck, with neural-network-specific approaches reviewed by Martens (2016).Source: §7 · p. 13
• The best optimization model need not be the most locally accurate: the Fisher may define longer-view quadratic models that explain its observed superiority to the Hessian in neural-network optimization.This possibility is presented as a speculative explanation, with further discussion deferred to Section 8.1.Source: §7 · p. 13

8. The Generalized Gauss-Newton Matrix

• The GGN generalizes Gauss-Newton to any loss convex in the network output, remains positive semidefinite, and equals the Hessian when predictions are optimal or when model noise makes the approximation exact.For optimal predictions, a second-order method using GGN approaches Newton’s method near convergence; a weaker equivalence holds under intrinsic target uncertainty.Source: §8 · pp. 14–15

8.1.1 Qualitative Observations

• Unlike the Hessian, the GGN avoids negative curvature, making its quadratic models more reliable and its curvature estimates more robust on data subsets because positive contributions cannot cancel.The GGN is not an upper bound on the Hessian: it omits curvature from the network function f while retaining curvature from the loss L, which is considered more stable.Source: §8.1.1 · pp. 15–16

8.1.2 A More Detailed View of the Hessian vs the GGN

• The Hessian decomposes into loss curvature, inter-layer interaction terms, and activation-function curvature; GGN retains the positive-semidefinite loss term while discarding components whose signs can change unpredictably.Activation and interaction curvature can vary rapidly with parameters and predictions, whereas loss curvature remains positive semidefinite and is expected to be more stable when averaged over training data.Source: §8.1.2 · p. 17

8.1.3 Some Insights From Other Works

• Prior work shows that GGN-based quadratic approximations preserve important global properties of nonnegative losses, coincide with Hessian diagonal blocks for piecewise-linear networks, and are supported by wide-network linearization analyses.These results suggest that modeling neural-network curvature omitted by GGN may be unnecessary or even counterproductive for optimization.Source: §8.1.3 · pp. 17–18

9.1 Computing the Fisher (and Matrix-Vector Products With It)

• Directly evaluating or inverting the Fisher is generally infeasible for high-dimensional neural networks, except in unusually simple model and distribution settings.Replacing Q_x with its empirical counterpart makes the Fisher expression easy to compute and closely related to the generalized Gauss-Newton matrix.Source: §9.1 · pp. 18–19
• Matrix-vector products with the Fisher can be computed efficiently using a linearized forward pass, multiplication by F_R, and standard backpropagation.This approach is practical when the predictive-distribution Fisher F_R is easy to compute.Source: §9.1 · p. 19

9.2 Qualified Equivalence of the GNN and the Fisher

• The Fisher and generalized Gauss-Newton matrices are equivalent when the predictive distribution has suitable curvature, notably for exponential-family models parameterized by their natural parameters.This includes Gaussian mean parameterizations and multinomial-softmax models, corresponding to squared-error and cross-entropy losses.Source: §9.2 · pp. 19–20
• When the loss Hessian is indefinite or the generalized Gauss-Newton matrix is undefined, using the Fisher or natural-parameter reparameterization yields the same positive-semidefinite curvature matrix.For Gaussian losses parameterized by mean and log variance, the loss Hessian can be indefinite, whereas the Fisher remains usable; placing softmax within the loss also restores Fisher–GGN equivalence under the stated conditions.Source: §9.2 · pp. 20–21

Step-size selection

• Natural-gradient step-size selection is difficult: heuristic adaptive schemes lack firm guarantees, while α_k = ρ/k has stochastic convergence guarantees but may perform poorly in practice.Practical optimization therefore cannot rely on a theoretically valid schedule alone.Source: §10 · p. 21

Discrete-update limitation

• Straight-line natural-gradient steps can veer far from their target in distribution space because the parameterization’s path is generally nonlinear over non-negligible distances.Figure 1 illustrates this failure mode for large discrete updates in the original parameter space.Source: §10 · p. 22, Figure 1 · p. 22

Second-order interpretation

• When the Fisher equals the GGN, the negative natural gradient with α = 1 is the optimal update for a GGN-based local quadratic model, but that model can become inaccurate over long distances.Small Fisher eigenvalues can also produce very large proposals, making uniform step-size reduction too crude because it rescales every eigen-direction equally.Source: §10 · p. 22

Damping methods

• This second-order view motivates update damping techniques that constrain or penalize the local quadratic-model solution, including Tikhonov regularization, trust regions, structural damping, and Krylov Subspace Descent.These methods address update proposals more selectively than simply reducing the global step-size.Source: §10 · pp. 22–23

Practical damping

• Tikhonov-damped GGN updates with step-size 1 made substantially more optimization progress than undamped updates, which required carefully chosen step-sizes to remain feasible.This practical evidence comes from Hessian-free optimization using the exact GGN, equivalent to the Fisher in that work.Source: §10 · p. 23

11. The Empirical Fisher

• The empirical Fisher replaces the model-distribution expectation with a target-distribution or empirical expectation, producing a positive-semidefinite matrix that is essentially free when gradients are already computed.It can also be applied to objectives without an obvious probabilistic model, but it is not generally equivalent to the Fisher or Gauss–Newton matrix.Source: §11 · p. 23, §11.1 · p. 23

11.1 Comparisons to the Standard Fisher

• The empirical Fisher has rank at most |S|, often simplifying diagonal or block computations, whereas stochastic estimators can also make corresponding standard-Fisher computations practical.Nevertheless, the true Fisher is preferable when possible because it approximates the Hessian in important cases and can equal the generalized Gauss–Newton matrix under stated assumptions.Source: §11.1 · pp. 23–24
• In the quadratic example, the empirical Fisher is θ^2 while the exact Fisher is 1, causing failure to converge unless ξ < 1 and step sizes decay sufficiently fast.Even when convergence occurs, its rate is limited by the step-size decay, typically O(1/k) or slower.Source: §11.1 · p. 24
• By contrast, the exact Fisher achieves fast linear convergence with rate |1 − α| for any fixed step size satisfying 0 < α < 2.This deterministic result cannot generally hold in realistic stochastic settings, but it remains desirable when gradients are accurately estimated.Source: §11.1 · pp. 24–25

11.2 A Discussion of Recent Diagonal Methods Based on the Empirical Fisher

• Diagonal AdaGrad, RMSProp, Adam, and related methods estimate empirical-Fisher diagonals, but historical averaging changes convergence behavior and may avoid failures of the current-iterate empirical Fisher.Exponentially decayed averages generally perform better pre-asymptotically because they forget stale contributions, whereas averaging all past gradients can shrink too slowly.Source: §11.2 · pp. 25–26

11.3 The Constants λ and ξ

• The damping parameter λ is best understood as Tikhonov regularization that keeps updates within a region where the local quadratic model remains reliable, requiring adaptive adjustment during optimization.The exponent ξ has been chosen as 3/4 for conservative diagonal preconditioning or 1/2 for regret proofs, while ξ = 1 is supported by an approximate-Newton method under stronger assumptions.Source: §11.3 · pp. 26–27
• Hazan et al. achieved O(log(k)) regret with ξ = 1, compared with AdaGrad’s O(·) bound, but their analysis used stronger hypotheses and the key advantage concerned hidden constants.Their approximate-Newton method’s constant avoided dependence on the Hessian condition number that affected the analyzed SGD version.Source: §11.3 · p. 27

12. A Critical Analysis of Parameterization Invariance

• In the infinitesimal-step limit, equivalent update directions produce the same path after mapping between smoothly reparameterized coordinates.For affine reparameterizations, the equivalence extends exactly across arbitrary finite sequences; more generally, sequences converge to equivalent smooth paths as α →0.Source: §12 · pp. 28–29
• Finite-step natural-gradient methods are not exactly invariant to smooth invertible reparameterizations, but their error decreases with step size and reparameterization curvature.The error is zero for affine reparameterizations and scales as α^2 under the theorem’s conditions.Source: §12 · pp. 29–30

12.1 When is the Condition J⊤ζ BθJζ = Bγ Satisfied?

• Curvature matrices built from parameterization-invariant output Jacobian forms satisfy the required transformation condition, including the Generalized Gauss–Newton, Fisher, and empirical Fisher matrices.The construction assumes the reparameterization function is independent of x and y.Source: §12.1 · p. 30
• The Hessian generally fails the sufficient transformation condition, except in special cases such as affine reparameterizations or zero gradients across all training examples.Thus Hessian-based updates do not generally inherit the same parameterization-invariance property.Source: §12.1 · pp. 30–31

• The natural gradient minimizes Fisher-metric squared distances between output changes optimal for the loss and those induced by parameter updates.The optimal output changes follow the natural-gradient direction for a KL objective and, under the GGN interpretation, the loss’s second-order Taylor approximation.Source: §13 · pp. 31–32
• The nonnegative discrepancy terms yield an upper bound on improvement predicted by the local quadratic model.When FR = HL, the bound equals the optimal improvement predicted by the second-order loss model if each case’s output z can change independently.Source: §13 · p. 32
• The natural gradient remains definable without F^-1 when F and ∇h use the same data and each FR is invertible.It can be selected as the minimum-norm minimizer of M(δ), computed in practice with F’s pseudoinverse, provided parameterization-invariant choices are maintained.Source: §13 · p. 32

14.1 Amari’s Fisher Efficiency Result

• Amari’s Fisher-efficiency result reaches the Cramér-Rao variance bound for shrinking 1/k stepsizes, including mini-batches, but assumes convergence, full-distribution Fisher estimation, realizability, and applicable density conditions.Mini-batch gradients are allowed, whereas mini-batch Fisher estimates generally invalidate the proven result; faster training-set convergence can represent over-fitting.Source: §14.1 · pp. 32–33
• With Polyak averaging, the dominant O(1/k) objective term matches the Cramér-Rao bound for any curvature matrix, while learning-rate decay can improve its constant but not the overall rate.The starting-point-dependent term scales as O(1/k^2) or better and may matter more under limited iteration budgets.Source: §14.1 · pp. 35–36

14.2 Some New Results Concerning Asymptotic Convergence Speed of General Stochastic 2nd-order Methods

• For the stochastic iteration on a quadratic objective, Theorem 3 gives mean and variance expressions for constant rates satisfying αλ1(B−1H∗) ≤ 1 and rates αk = 1/(k+a+1), with a ≥ 1.For the latter rate, the variance includes a matrix-valued error Ek that shrinks as 1/k2.Source: §14.2 · p. 37
• Theorem 3 shows that on convex quadratics, E[θk] follows deterministic optimization independently of gradient-noise covariance, while the covariance affects only variance around that mean.This follows because quadratic gradients are affine and therefore commute with expectation.Source: §14.2 · pp. 37–38

14.2.1 Some Consequences of Theorem 5

• Under decaying stepsizes, Hessian preconditioning can have a smaller noise-dependent asymptotic term than first-order scaling, especially when H∗ is ill-conditioned, though ill-conditioning alone is insufficient.For fixed stepsizes, the expected objective approaches a noise- and stepsize-dependent constant.Source: §14.2.1 · pp. 39–40
• With B = H*, noise scaling is no worse than twice the B = λ_n(H*)I case and can be much smaller for broad spectra.For λ_i(H*) = n − i + 1, the B = λ_n(H*)I noise scales as Ω(n^2/k).Source: §14.2.1 · pp. 40–41

14.2.2 Related Results

• Prior analyses are more general but less detailed, whereas Theorem 5 characterizes asymptotic behavior specifically for convex quadratics.The paper states that no known result totally subsumes Theorem 5 for convex quadratics.Source: §14.2.2 · p. 41

14.3 An Analysis of Averaging

• Exponential moving-average weighting forgets early iterates faster than basic averaging, but its variance does not converge exactly to zero.The practical variant uses β_k = min{1 − 1/k, β_max}, with β_max close to 1.Source: §14.3 · p. 42

14.3.1 Some Consequences of Theorem 6

• Under realizability, stochastic gradient descent with averaging achieves n/(k + 1) + o(1/k), matching Fisher-efficient methods despite using no curvature matrices.Second-order methods improve the noise-independent term, especially when H* is ill-conditioned and the initialization has components along low-eigenvalue directions; this can matter before asymptotic noise dominates.Source: §14.3.1 · pp. 43–44

14.3.2 Related Results

• Compared with Theorem 6, related averaging bounds may fail to establish convergence or distinguish noise-dependent from noise-independent contributions.The paper presents Theorem 6 as more detailed for convex quadratics, including non-asymptotically dominant terms.Source: §14.3.2 · pp. 44–45

15. Conclusions and Open Questions

• Natural gradient should be treated as a GGN-based second-order method, using trust regions, damping, Tikhonov regularization, and Levenberg-Marquardt heuristics for robust performance.This interpretation assumes equivalence between the Fisher information matrix and the GGN.Source: §15 · p. 45
• Rigorous global convergence for GGN optimization remains unresolved, while existing quadratic local bounds favor the Hessian and cannot explain GGN superiority or damping’s practical importance.These bounds assume constant curvature and therefore exclude benefits from alternative curvature matrices.Source: §15 · p. 45
• The Fisher–GGN–Hessian relationship remains difficult to interpret in non-realizable settings, including bounds whose constants become opaque after substituting Fisher curvature.The report leaves open whether GGN advantages over Hessians, Fisher–GGN differences, and stochastic second-order pre-asymptotic gains can be rigorously characterized.Source: §15 · p. 46

Appendix B. Proofs of Convergence Theorems

• Appendix B proves Theorems 3, 5, and 6 by deriving a variance recursion for stochastic second-order optimization and analyzing its operator properties.The recursion for V_k is central to the convergence analysis.Source: Appendix B · pp. 47–48

B.1 The Constant αk = α Case

• For constant and diminishing step sizes, the appendices establish positive-spectrum and PSD-preservation properties, derive variance formulas with error terms, and bound those errors.The constant-step-size case uses Λ and (I − Λ)^−1, while the diminishing schedule uses Ξ − I and an error recursion.Source: §B.1 · pp. 48–49, §B.2 · p. 53, Proof · p. 56

Proof

• The proof establishes upper and lower bounds for the expected objective value using matrix spectral inequalities and iterative recurrences.The derivation uses PSD/NSD matrix properties, eigenvalue bounds, and propositions to obtain the stated bounds.Source: Proof · p. 58

B.3 Proof of Theorem 6

• Theorem 6 analyzes averaged iterates from a constant-step-size stochastic iteration on a quadratic objective under αλ1(B−1H*) < 1.The proof uses unbiased stochastic gradients and conditional independence to derive the variance of the averaged parameter.Source: §B.3 · p. 61, Proof · pp. 61–63

Appendix C. Derivations of Bounds for Section 14.2.1

• Appendix C bounds quadratic-form terms using λn(H*) and tr(Σg), with one bound scaled by the condition number κ(H*) = λ1(H*)/λn(H*).The displayed bounds include ≤ 1/2 λn(H*) tr(Σg) and a corresponding condition-number expression.Source: Appendix C · p. 66

Appendix D. Some Self-contained Technical Results

• Lemma 17 derives a trace identity from A^T P + PA + Q = 0 when P is invertible.Premultiplication by P^-1 and trace cyclicity reduce the equation to a relation involving tr(A) and tr(P^-1Q).Source: Appendix D · p. 67
• Lemmas 18–20 establish trace inequalities for PSD or NSD matrices and commuting PSD matrices ordered by the Loewner relation.The NSD result follows by applying the PSD case to -X, while the commuting-matrix result uses positivity of A^1/2(T-S)A^1/2.Source: Appendix D · p. 67
• Lemma 21 shows that Φ(X) = X − DXD^T is invertible when D has real eigenvalues strictly between −1 and 1, and preserves positive semidefiniteness under inversion.The proof represents Φ as I − D ⊗ D; its eigenvalues are positive, and Φ^-1(X) is a convergent sum of PSD matrices for PSD X.Source: Appendix D · p. 68
• Proposition 22 bounds harmonic numbers using logarithms, the Euler–Mascheroni constant, and finite-size correction terms.The resulting inequalities support later bounds on differences of harmonic numbers, with equality in the identified equal-index case.Source: Appendix D · pp. 68–69
• Proposition 23 converts products involving nonnegative parameters into exponential and polynomial bounds depending on a, b, i, and k.The derivation combines a termwise bound with Proposition 22 and introduces ν(a) = (a + 2)^3/(a(a + 1)^2).Source: Appendix D · pp. 69–70

Appendix E. Proof of Corollary 2

• Appendix E proves that additive optimizers using B^-1-scaled gradients follow the same path under θ- and γ-parameterizations, given equivalent starting points and the stated regularity conditions.The affine case follows from Theorem 1; the non-affine case assumes uniform continuity and bounded normalized updates as α → 0.Source: Appendix E · pp. 70–71
• For k ≤ T/α, the parameterization discrepancy converges to zero as α → 0, because its recursion is α^2E + (1 + αU) times the previous discrepancy.Starting from matched initial points, the geometric-series bound and (1 + αU)^(T/α) → exp(UT) yield convergence over continuous paths of finite non-zero length.Source: Appendix E · p. 71