Convergence Properties of Natural Gradient Descent for Minimizing KL Divergence
Authors: Adwait Datar
Nihat Ay
Venue: Accepted by TMLR ()

Abstract: The Kullback-Leibler (KL) divergence plays a central role in probabilistic machine learning, where it commonly serves as the canonical loss function. Optimization in such settings is often performed over the probability simplex, where the choice of parameterization significantly impacts convergence. In this work, we study the problem of minimizing the KL divergence and analyze the behavior of gradient-based optimization algorithms under two dual coordinate systems within the framework of information geometry$-$ the exponential family ($\theta$ coordinates) and the mixture family ($\eta$ coordinates). We compare Euclidean gradient descent (GD) in these coordinates with the coordinate-invariant natural gradient descent (NGD), where the natural gradient is a Riemannian gradient that incorporates the intrinsic geometry of the underlying statistical model. In continuous time, we prove that the convergence rates of GD in the $\theta$ and $\eta$ coordinates provide lower and upper bounds, respectively, on the convergence rate of NGD. Moreover, under affine reparameterizations of the dual coordinates, the convergence rates of GD in $\eta$ and $\theta$ coordinates can be scaled to $2c$ and $\frac{2}{c}$, respectively, for any $c>0$, while NGD maintains a fixed convergence rate of $2$, remaining invariant to such transformations and sandwiched between them. 
Although this suggests that NGD may not exhibit uniformly superior convergence in continuous time, we demonstrate that its advantages become pronounced in discrete time, where it achieves faster convergence and greater robustness to noise, outperforming GD. Our analysis hinges on bounding the spectrum and condition number of the Hessian of the KL divergence at the optimum, which coincides with the Fisher information matrix.

Published in Transactions on Machine Learning Research (07/2025)
Convergence Properties of Natural Gradient Descent for Minimizing KL Divergence
Adwait. Datar
adwait, datar datar&tuhh....... de
Nihat Ay1.2,3
nihat, ay @tuhh.de de
Institute for Data Science Foundations Hamburg University of Technology 21073 Hamburg, Germany
Santa Fe Institute
Santa Fe, NM 87501, USA
Leippig University 04109 Leipzig, Germany
Reviewed OpenReview: https: httpp:://openreviaau net/ forum? P id= J6F5)
openreview.
Abstract
The Kullbeck-Loibeer (KL) divergence plays central role in probabilistic machine learning, where it commonly serves as the canonical loss function. Oplimization in such settings is often performed over the probability simplex, where the choice parameterization signif- icantly impacts convergence. In this work, we study the problem of minimizing the KL divergence and analyze the behavior of gradient-basse optimization algorithms under two dual coordinate systems within the framework of information geometry the exponential family (0 coordinates) and the mixture family (1) coordinates). We compare Euclidean gradient descent (GD) in these coordinates with the coordinate-invarianl natural gradient descent (NGD), where the natural gradient is a Riemannian gradient that incorporates the intrinsic geometry of the underlying statistical model. In contimuous time, we prove that the convergence rales of GD the and 7) coordinates provide lower and upper bounds, respectively, on the convergence rate of NGD, Moreover, under affine reparameterization of the dual coordinates, the convergence rates of GD in 17 and coordinates can be scaled to 2c and 2 respectively, for any > 0, while NGD maintains fixed convergence rate of 2, remaining invariant to such transformations and sandwiched between them. Although this suggests that NGD may not exhibit uniformly superior convergence in continuous time, we demonstrate that its advantages bocome pronounced in discrete time, where it achieves faster convergence and greater robustness to noise, outperforming GD. Our analysis hinges on bounding the spectrum and condition number of the Hessian of the KL divergence the optimum, which coincides with the Fisher information matrix.
1 Introduction
The convergence properties of the natural gradient descent a algorithm, originally introduced in Amari 1996). have been extensively studied in the literature (e.g., Amari [1998); Pascanu & Bengio (2014); Martens (2020g). In particular, the natural policy gradient Kakade (2001) ) has motivated a rich body of rescarch (sec, e.g., Müller Montúfar (2024); Yuan al. (2022); Khodadadian et al. (2022).) Beyond this, the natural
gradient methods have been applied to diverse problems including Bayesian networks
2023
Ay & van

---

Published in Transactions on Machine Learning Research (07/2025)
Oostrum 2023)
over-parametrized neural networks Zhang et al. (2019 van Oostrum & Ay [2021) van
Oostrum et al. (2023) and infinitel-wide networks Karakida et al. (2019) Karakida & Osawa 2021 to
name a few. Related to our focus, recent work has also investigated convergence rates of natural gradient
flows and their discrete counterparts (sce Zhang et al. (2019); Xiao 2022 Yuan et al. 20222; ; Khodadadian
et al. 2022) : Müller úonr (2024)). ) A commonly observed phenomenon is that natural gradient descent outperforms Euclidean gradient descent, albeit at higher computational cost. In this work, we revisit this comparison in simple yet illuminating setting: minimizing the Kullback-Leibler (KL) divergence over the probability simplex. The KL divergence is fundamental loss function in probabilistic machine learning, arising naturally from the maximum likelihood principle and information-theoreii considerations (Mohri et al. 2018 Section 12.1.1). Despite the apparent simplicity of the problem, we observe that natural gradient flows do not universally outperform standard Euclidean gradient flows.
Specifically, we consider two dual parametrizations of the probability simplex: the exponential family repre- sentation (the coordinates) and the mixture family representation (the n coordinates) Amari (2016). We prove that the natural gradient flow converges faster than the Euclidean gradient flow in the 0 coordinates (the 0-gradient flow), consistent with results in the literature. However, the natural gradient flow (despite yielding straight-line trajectorics) converges more slowly than the Euclidean gradient flow in 1j coordinates (71ggadient flow). This demonstrates that the often-reported rapid convergence of natural gradient flow cannot be simplistically attributed to the straightness of its trajectories. Furthermore, leveraging the invari- ance of the canonical divergence under affine transformations, we show that the convergence rates of the Euclidean gradient flows in 1 and coordinates can be adjusted to 2c and e* respectively, for an arbitrary <> 0, while the natural gradient maintains fixed convergence rate of 2, sandwiched between them. Thus, by setting c 1, we can construct pair of dual coordinates, (77. 0); that match the convergence rate of the natural gradient flow. Since the advantages of the natural gradient are not immediately apparent in the contimuous- time setting, we extend our analysis to the discrete-time case, where the natural gradient demonstrates both faster convergence and greater robustness to noise, outperforming Euclidean gradient descent. We show that the fundamental reason behind the superiority of natural gradient lies in the optimal conditioning of the loss landscape: the natural gradient updates are equivalent to minimizing the loss func- tion 1l17 - 11, whose Hessian has a condition number equal to The main contributions of this paper are summarized as follows:
1. We analyze the convergence rates of Euclidean gradient flows in 1) and 0 coordinates, and of the natural gradient flow (Thcorem 3 and Theorem 9h. We show that while the natural gradient flow converges faster than the 0-gradient flow, it is slower than the r-gradient flow. This result builds upon bounds on the spectrum of the Hessian of the loss function established in Lemma 2 These theoretical findings are supported by illustrative numerical experiments.
2. Exploiting the duality and the invariance of the canonical divergence under affine transformations, we demonstrate in Theorem g that the convergence rates of Euclidean gradient flows in 7 and 0 coordinates can be adjusted to 2c and 2. respectively, for an arbitrary e >> 0, while the natural gradient maintains a fixed convergence rate of 2, sandwiched between them. This shows that there exists pair of dual coordinates, (77 9) such that the convergence rate of 7 and f gradient flows matches the convergence rale of the natural gradient flow.
3. We analyze the discrete-time dynamics in Section 4 where Theorems and 3 establish the superior robustness properties of natural gradient dynamics compared to their Euclidean counterparts. The core reason for this superiority is the optimal conditioning of the underlying loss landscape. Tn particular, natural gradient updates can be interpreted as minimizing the loss function _ fllm - 112, whose Hessian exhibits the ideal condition number of 1.
4. We complement our theoretical results with empirical studics in Section 5 extending the analysis to practical settings where only finite samples from the larget distribution are available. Specifically, we consider optimization of the empirical KL divergence in both full-batch and stochastic gradient descent (SGD) scttings. Our results show that NGD consistently outperforms standard GD when learning rates are optimally tuned, in alignment with Theorems 5 i and
g Furthermore, when using sufficiently small and equal learning rates across all methods, the sandwiching behavior predicted
2

---

Published in Transactions on Machine Learning Research (07/2025)
by Theorem 8 for continuous- time dynamics persists in the discrete-time, sample-based setting thus validating the relevance of our theory beyond idealized assumptions.
Notation
Let I denote the set of real numbers. For function of two variables I= and E IR",, g Rd IRm R. we use the notation
(z, v) Úz
Zg
0' (x.
8'
(I, y)
021
V vzg(z,v)~
3g 32
(x), y)
(I, v)
8 (z, y) 3u
For functions of single variable I, suppress the subscript and simply write ff(x) and V f')(). For manifold M with two global charts @m .M - &m (M) C R" and &. M + 6c(M) C R" with coordinates (MM and f E de (MM, we slightly abuse notation and write for any smooth function L: M R
[(n) == f(
[()) == [()-
V == v (Co ( (n),
VC(0) = ((d-
CL) == V (Lo ) (1).
V²L(0)= == V (E) (0).
For a point p M, we write 7pp = (m)p) and 0p = (p) to denote the point p in the 7 and 0 coordinates, respectively. For a symmetric matrix Q, we write Q (resp. = 0) to denote that Q is symmetric positive definite (resp. positive semi-definite) Building
this notation, we wrile Q > P (resp. Q > P) to mean Q Q > (resp. Q P = 0). For symmetric matrix Q, let Amin (Q) and Amax (Q) denote the minimum and maximum eigenvalue of Q. Since Q > implies that all eigenvalues of are positive, A@sr(Q)
can define the condition number of Q as cond(Q) = (Q) For any R", let ||>|| denote the standard Fuclidean norm, and define the closed norm ball of radius centered at x by (x) == {v E R" Ily - x|| I =}. For any real matrix M, let || = suP2#0 eigenvalue of M when M 0. ||=||
|Mall be the induced matrix 2-norm, which coincides with the maximum
2 Information Geometry Preliminaries and Gradient Flow Dynamics
In this section, we review the information geometric preliminaries and arrive at the continuous-time gradient flow dynamics which are then analyzed the following section. For further details on the underlying concepts
of information geometry, the reader is referred Amari (2016); Amari & Nagaoka (2000): Ay et al. (2017)
2.1 Discrete Distributions in Mixture and Exponential Coordinates
Consider the family Sn of probability distributions over a discrete random variable X with sample space Ω (1,2. n,n n1n+11. Let p; be the probability that X takes the value i. Then, any p E S, can be written
p(x) Σ p.),(x),
where (x(x) is the delta distribution over Ω concentrated at i. Thus, S, can be identified with the n-dimensional simplex i..., S, = ( P2. Pn. Pu+1) ) E RR*1|| p; 0, E Pt = =1 This family ad- mits representations both as mixture family and an exponential family Amari (2016). This can be seen by noticing that any p Sn can be written
pxx) Σ nd((x) + -1-
Σ Tk 8u+1(2) exp -+(0) Σ (1)
Mixture family representation
Exponential family representation
Note that definition the simplex excludes the boundary.
3

---

Published in Transactions on Machine Learning Research (07 (07/202) 2025)
m
122
Tee
m
1/1
01 02
c
+
4
14 00 + c0
1/1
72
Figure 1: Left: Coordinate system with the natural parameters = (71. of the mixture family represen- tation of S2. Right: Coordinate system with the natural parameters 0 = 02) of the exponential family representation of S2.
where (0) log (1+E e is the log partitiio function ensuring the normalization constraint ExEn (o( 1 for the exponential family representation. With nj (m, "2, ,7,),, we obtaín ordinate system for the simplex, with serving the natural parameter of the mixture family. We let om Sn = //, ... nm) E R= denote the global chart for S, in the mixture family coordinate
system. This is depicted in Fig. I (left). Similarly, with 0 == (0,, 12,00,,... 0n), we obtain an alternate coordinate system for the simplex with being the natural parameter of the exponential family. Define de Sn 3 p + = (01, ... ..,) R. as the global chart for S, in the exponential family coordinate system. This is depicted in Fig. (right).
2.2 Convex Duality and Bregman Divergence
For the family Sa, there exists dual relationship between the coordinates n) and 0. Since v0) log (1+ Ei= eat is strictly convex, it is possible to define its convex conjugate .(π) maxo (n 0 (0)). Optimality condition on the maximizer dopt = 0 yields the relationship Vab(0) = 7,, which can be solved to
obtain 0, log
This results in ((n) (n x())) = Σn log n. the negative of Shannon
entropy. Conversely, & the convex conjugate of & leading to = V)((). Since V, is the identity map, application of the chain rule gives
v2p(n) [V2x(0]]
(1)
The convex conjugate functions and yo define pair of Bregman Divergence D& and D, satisfying
(0p || 0g) = (((p))) - (l,)) Vo(0,) - 0g)
(2)
=p(())) (p(1)) vep(n)) mp) =: Di 11 1p).
(3)
4

---

Published in Transactions on Machine Learning Research (07/2025)
In our setting of Sn this Bregman divergence equals the canonical KL-divergence D(g between probability distributions and p, Le.,
a+1 log
De (0p || 0g) = ( Il mp) ((all))
Σ qi
qi
(4)
Pj
i=1
where (Ng Tp) and (0,, 0p) are the coordinate representations of (q.p) p) in the 7j and 0 coordinates, respectively, For further details, the reader is referred to Amari & Nagaoka (2000).
2.3 Gradient and Natural Gradient Dynamics
For given target distribution Su, let the loss function Lg S, - R be defined by Lg(p) = D(glp)). As discussed in Section | we abuse the notation slightly and interchangeably use q. 0g = (<) or = (m) to denote probability distribution E Sn, the same distribution in coordinates and in 1j coordinates, respectively. Thus, we write L, to mean Lo (mn)) and (0,), to mean L,(0=' (0,)). The gradient of the loss function can be computed the different coordinate systems using equation 2 and equation as
(lg) - (xp(p)) ( mp)) + Vy((p)= =- (pp))))) mpp)
(5)
(0p) = V ve(0,) - (o,))
(6)
Analogously, for given target distribution p E Sn, let the loss function C : Sn - be defined by ()(g) D(gl (p). The gradient of this loss function can be computed in the different coordinate systems
using equation 2 and equation 4
(mn) = Vp(nn)) - (p((pp)
(7)
(0,) = (v,()) - (20(0,)(0,-- - 0gg) + VV((g) =0 -gy)
(8)
Building on the pair of conjugate dual functions, we define a Riemannian metrie 9 on Sn, which assigns to each point p E S, an inner product (, on the tangent space TyS.. The tangent space T,S,, can be identified with { (1) (,,,,,,... u2,... Uyy 9x+1) E R ETH 0; =} (sce Ay et al.
(2017)) In coordinates, this
Riemannian metric is represented by the matrix V²p(()) in the 7 coordinates and v=v(0) in the coordinates. Importantly, this Riemannian metric coincides with the Fisher metric (Amari) 2016 Theorem 2.1). Using this Riemannian metric, we define the Riemannian gradient grad (.)(p) also called as the natural gradient, at point p Sn through the relation
(grad Lv(,), u) = (p)pl
(9)
for all in the tangent space of S, at the base point p. This defining property allows us to compute the natural gradient in n coordinates using equation as
grad (mp) [V>(()))] (c, (ng -mp)
Similarly, grad (g) can be computed in the 0 coordinates using equation g as grad (=,0))=-00 0,). Interestingly, the natural gradients when represented in appropriate coordinates take on particularly simple linear forms- - they directly point towards the target distributions. Since the situation with the loss function Cp is analogous to that of Ly. for brovity, we will focus our analysis on Lg for the remainder of the paper
We now introduce the gradient flow dynamics in both n and coordinates, as well as the natural gradient flow which will be analyzed in the subsequent sections. Although the natural gradient flow dynamics are coordinato-invariant, we express them in 1) coordinates to exploit the particularly simple lincar structure.
For a given larget distribution q E S, and an initial distribution po E Sn, consider the gradient flow dynamics described by equation 10 and equation II and the natural gradient flow dynamics described by equation
12
(t() VC,(((t)),
(7(0) ?pp
(10)
0() (O()),
((0) = Opo:
(11)
ing (t) = grad Eg (Tig (t)). (0) = 1pp
(12)
We include the convergence rate ana analysis of gradient flows for CG in Appendix for completeness.
5

---

Published in Transactions on Machine Learning Research (07 (07/202) (025)
We will analyze these dynamics in the following sections.
3
Convergence Analysis in Continous Time
We start the convergence analysis with general result which is at the core of the analysis. This result frequently appears in various forms, typically emphasizing the upper bound in inequality 14 (see (Wensing & Slotine 2020 Proposition 1) for example). We present an adaptation of this result to our sctting.
Proposition 1 (Convergence of general gradient flows). Consider the gradient flow dynamics
+(l) - ())((()) x<lo) -o.
(13)
where sufficiently smooth function f U > IR, with U C R being open neighborhood of xo, satisfies the following properties:
(a) has unique minimizer v. E U satisfying Vf(x.)= 0.
(b) there exist positive constants m and L such that m 1 f2)(z) 3 L- for all x in the sublevel set Ulf(z) 2 f(xo)}
Then the solution I [oo, 20) + U of equation 13 satisfies
(1) a(t) E S for all 2 to and
(i) With c = f(xo)= - f (I.),
we get that
2m((t-to) e for all 2 to. c. e 2L(t- s f(x(t)) - f(x.) S c.c
(14)
i.e., f(x(t)) converges exponentially to fxx) with a rate larger than 2m and smaller than 2L.
Proof. See Appendix C.1
To facilitate the application of Proposition 1 to the dynamics given in equation 10 and equation establish bounds on the Hessian the loss function in the following Lemma.
we
Lemma 2 (Bounds on the Hcssian of the loss function). Let p. 9 E Sn. Then the following statements hold:
(1) (Global bound) The Hessians of by and C sutisfy
(0) < V (g)
VO (O om (Sa)
and
(15)
< (n)
Vn (S.).
(16)
(8) (Local bound at optimum) The Hessians of Cy and Cg. evaluated at the optimum satisfy
I v (Tg) = V2 (g()))
(17)
I C = V2 (,0p)))
(18)
Proof. See Appendix C.2
The positive definiteness of the Hessians established in Lemma 2 (scc equation 15 and equation 16 shows that E, is convex in the 0 coordinates as well as in the 1 coordinates, and C is convex in the coordinates. One might conjecture that C; might likewise be convex in the coordinates. However, this turns out not to be the case, as illustrated by counterexample. Figure shows sample contour plot of (7 for n 2. The plot clearly shows that the sublevel sets of (0) are non-convex, disproving the conjecture. To summarize, while the KL divergence is geodesically convex along m-geodesics in both of its arguments, it is geodesically

---

Published in Transactions on Machine Learning Research (07/202) 2025)
0.9
0.8
2
0.7
0.6
0
0.5
0.4
-1
-2
0.2
-3
0.1
-4
-2
0
2
4
θ1
Figure 2: Contour plot of (0) showing that the sublevel sets are non-convex.
convex along e-geodesics only with respect to its second argument (sce Boumal 2023 Definition 11.3) for definition of geodesic convexity and (Amari 2016 Section 2.4) for definitions of and m- geodesics). With these uniform bounds on the Hessians in place, we are now equipped to control the exponential decay rates of gradient flows through Proposition This is presented in the next result which is the main result of this section.
Theorem 3 (Convergence analysis). Let E Sn be the larget distribution and po E Sn be the initial distribution. Suppose n, 0 and Thg be the solutions to dynamics described by equation 10 equation 11 and equation respectively. Then
ii) there erist positive constants < my S Leg Cgy En and T such that
- < (n((1)) S - 2mat -2mpt S E 25 CGe 2ba
VIT
(19)
i.e., (n(()) converges to zero exponentially with rate higher than 2. Furthermore, if (,((()) is
sufficiently small, the result holds with T =.
ii) there erist posilive conslants mg S La < 1 and cy such that
Cae 21 < c&e 2Le!
-Lat < (,((t)) S cge 2m Vt 1+20 0,
(20)
i.e. Lg (0(1)) converges zero erponentially with rate lower than
(tti) there crist positive constants C1 and c2 such that
cie 2 < (y(gg)()) S oze
Vt 02
(21)
i.e,, (t))((t)) converges to zero exponentially with rate 2.
Proof. See Appendix
C.3
Theorem 3 shows that gradient dynamics in the mixture family coordinates exhibit faster convergence rates than natural gradient dynamics, which, in turn, outperform gradient dynamics in the exponential family coordinates. On one hand, this supports the generally obscrved superiority of the natural gradient dynamics over gradient dynamics in the exponential family coordinates. On the other hand, demonstrates that natural gradient dynamics slower than gradient dynamics in the mixture family coordinates. Note that

---

Published in Transactions on Machine Learning Research (07 (07/202) /025)
en(( *- (t) g-((t)
e D(gll((t) D(gllnag()) D(gll((t))
KL divergence
0.2
0.1
t
10 0 gradient flow de- scribed by equation 11 and the natural gradient flow described by equation 12 for = and E [0, 1.5] superimposed on the level curves of KT. divergence. The markers on the curves show equal time intervals for each curve. Right: KL divergence evaluated along the solutions plotted as function of time t, The intervals between markers the left figure are unrelated to the intervals between markers in the right figure.
Figure 3: Left: Simulation trajectories of 7/ gradient flow described by equation
although we choose to represent the natural gradient dynamics in the coordinates, the obtained convergence rate bound in equation 21 is independent of this choice. Furtheremore, that the convergence rate bound for the 1) coordinates established in equation 19 Theorem asymptotic nature and formally holds only after some time 7 > 0, whose value is not explicitly characterized. This is common feature of asymptotic convergence rate analyses, where rate guarantecs apply beyond some unknown T. Finally, stated in the last sentence of item (i), Theorem 3 the result can be interpreted as describing local behavior: if the initial distribution lies sufficiently close to the target, the bound holds from the start (i.e., with T = 0). More precisely, sufficient condition for the bound to hold with T = 0 is that (n) > 1 for all belonging to the sublevel set S, == {17 E 6m (S,) I (n)))) Despite the absence of a global estimate on T, our empirical results (see Figure for n and Figure for n 10) indicate that the ordering of convergence rates predicted by theory emerges early in the flow. This supports the practical relevance of the asymptotic comparison.
We now present numerical experiments with = 2 to illustrate the theoretical results developed so far. Figure 3 (left.) depicts the trajectories of the 77 graaddientt flow described by equation 10 0 1-graddient flow described by equation 11 and the natural gradient flow described by equation 12 superimposed on the level curves of KL divergence. Although the natural gradient flow follows straight trajectories, it is slower than the 77 -rgradient flow but faster than the A 8-gradient flow, as scen from the unit time markers along the curves. Figure (right) confirms this by plotting the KL divergence over time along these trajectories.
logarithmic scale, revealing exponential convergence. Best-fit lincar curves are superimposed to estimate the slopes, which correspond to the convergence rates. The -ggrddient flow exhibits the fastest convergence (slope = 7), the 0-gradient flow is slowest (slope 2 0.475), and the natural gradient flow lies in between (slope ~ 2.04), closely matching the theoretical prediction of rate 2. This experiment repeated over 100 randomly chosen initial conditions and one randomly chosen target distribution, as shown Figure 4 (right). The empirical convergence rates align well with the theoretical bounds from Theorem 3 confirming that n r-grrddient flows exceod rate 2, natural gradient flows converge at rate 2, and 8-gradient flows fall below rate 2.
To better highlight the convergence rates, Figure (left) presents the KL divergence
Finally, note that the nat ural gradient dynamics governed by equation 12 can be equivalently described
ing (1) -grad (ynn (1)) = "ng (1) -V (1
where (t)) = 11 Bag (t)l The natural gradient dynamics thus correspond to minimizing a convex quadratic function with identity Hessian.
8

---

Published in Transactions on Machine Learning Research (07/202) 2025)
empirical convv rate
8
log KT. divergence
@o0am o
20
r) gradient
natural gradient
0 gradient
log(D((lln((t)) log(D(allnn) (())
2
-40
(og(D(g))0(t))
20
40
60
80
100
4
data
Figure Left: KL divergence evaluated along the solutions to 17 -rgrrdient flow described by equation O O-gradient flow described by equation 11 and the natural gradient flow described by equation 12 for n plotted semi-log scale. The dashed lines show the best-fit linear function used to estimate the slope which gives the convergence ratc. Right: Empirical convergence rates for 100 randomly chosen initial distributions and a randomly chosen larget distribution for n 2. The dotted line shows the theoretical lower bound for the convergence rate of --ggradient. flows and the dashed line shows the theoretical upper bound for the convergence rate of 0 -rgradient flows.
10
-10-2
0.4
4
f(s) = 82
f(s) = 5
(g + 24)
Cule + s w)
0.3
(0, +
(m +
2
0.2
0.1
0
0.1
10
0
5
10 2 0.1
0.5
0.5
8
Figure 5: Left: Local sections of KL divergence around the minimizer g plotted as C. + 8. --,,, where all U, unit norm vectors distributed evenly on the unit circle. A quadratic function f(s) = 1 is also shown for reference. Right: Local sections of KL divergence around the minimizer plotted L. + w:). where all unit norm vectors distributed evenly the unit circle. A quadratic function s(s) 1/8 is also shown for reference.

---

Published in Transactions on Machine Learning Research (07 (07/202) 2025)
0
KT. divergence
0.6
log KL divergence
0.4
- D(alla(!))
-20
D(gllngg())
D(qll((t))
0.2
- log( ((glll(t))
40
log(D(gllin (gggDgll)ng((tt))
0
log (()()())))
3
0.5
1.5
2.5
Figure 6: KL divergence evaluated along the solutions to 7/ggrddeet flow described by equation 0 0-gradient flow described by equation plotted linear scale (left) and [ and the natural gradient flow described by equation [2 for semi-log scale (right). The intervals between markers in the left figure are unrelated to the intervals between markers in the right figure. The dashed lines in the right figure show the best-fit linear function used to estimate the slope which gives the convergence rate.
Figure 5 offers insight into convergence behaviors by plotting local sections of the KL divergence near the optimum along different directions. This is achieved by plotting the functions (1g + +.. for different L, + w:) on the unit circle for the 0 coordinates. This illustrates the local curvature of
directions t; the unit circle for the ") coordinates and by plotting the functions for different directions the function around the optimum. Since the natural gradient dynamics correspond to minimizing a convex quadratic function with identity Hessian as discussed above, we overlay the plots with a quadratic function (()) s for reference. The plots reveal that the functions 8 + L, + 8 v.)) exhibit higher curvature than the quadratic reference function f, while the functions s ~ (0, s-w;) appear flatter. This provides the core intuition behind the fast convergence of the gradient flow in coordinates in comparison to gradient flow in the coordinates and shows why the natural gradient flow falls in between the two extremes.
To extend our numerical study to a higher-dimension and | for the case n = 10. Figure G shows the KL divergence evaluated along gradient flow trajectories: the left panel shows the divergence on a linear scale, while the right panel uses semi-logarithmic scale before, best-fit lines are overlaid to estimate the slopes, which correspond to the convergence rates. The graddientt flow shows the fastest convergence (slope 16.239), the 8-gradient flow is the slowest (slope = 0.099), and the natural gradient flow lics in between (slope ss 1.974), closely matching the theoretical rate of 2. For greater confidence, repeat this experiment over 100 randomly chosen initial distributions and one randomly chosen target distribution; the results are shown in Figure 7 The empirical convergence rates align well with the theoretical predictions from Theorem 3 confirming that 7 gradieet flows exceed rate 2, natural gradient flows converge at rate 2, and 0 gradient flows fall below rate 2. Moreover, based on these empirical studies, we also observe that the difference convergence rates more pronounced for n 10 as compared to n
setting, we repeat the experiments from Figures 3
highlight exponential convergence.
Finally, the analysis presented in this section raises an interesting question: Since the dual pairing between the coordinates and 0 is preserved under appropriate affine transformations (as discussed in the following section and in Amari (2016), how do the convergence rates of the resulting dynamics change under such affine coordinate transformations? This question addressed in the following sul subsection.
3.1 Convergence Rate Analysis Under Affine Coordinate Transformation
We first review the effect an affine transformation of coordinates on the duality pairing between 0 and 77 (or equivalently between and ip). Since the convergence rate analysis from the previous scction hinges on bounding the Hessian of the loss function, we investigate how the Hessian transforms under an affine change of coordinates. Consider new coordinates 0 that related to the original 0-coordinates via affine
10

---

Published in Transactions on Machine Learning Research (07 (07/202) 2025)
-10-
empirical cony. rate
empirical conv. rate
gradient
gradient natural gradient
gradient
5
0
20
40
60
80
100
0
20
40
60
80
100
data
data
Figure 7: Left: Empirical convergence rates for 100 randomly chosen initial distributions and randomly chosen target distribution with 10. The dotted line shows the theoretical lower bound (12.94) for the convergence rate of --gradientt flows and the dashed line shows the theoretical upper bound (0.31) for the convergence rate of A- grardiinnt flows. Right: Same data, but with the y-axis restricted to highlight. the variation in convergence rates of --ggrddient flows.
transformation: 0= +b, where A E Rax= is an invertible matrix and b E RR is an arbitrary vector. Let & be defined as (0) = <AA0 + b). A simple application of the chain rule shows (2ú(0)) = ( Therefore, (0) > 0 if and only (2)((0) > 0, since A is non-singular. This implies the strict convexity of v, and it possible define its convex conjugate ((n) max, (n) (()).. Optimality condition on the maximizer Dops = yields the relationship 7 Viy(0) = A (1) = n.n. It is straight-forward to verify that with these newly defined convex conjugate pairs of functions &i and c the Bregman divergence still gives the original KL divergence, i.c.,
(0,,, 0,) Dy The Hessians of the loss function when evaluated the optimum, transform as follows:
D (0, || 0,) D& hp) = (0,
|l mp) D(lll)).
v voh Fg (Tay v2 = V ((g)4, (o,) V (0,, = A V (0,))
(22)
(23)
This calculation immediately gives us the following theorem.
Theorem 4 (Dual coordinates with identity Hessian). Let e be a positive conslant, E Sn be the target distribution and po E Sn be the initial distribution. There exists a pair of convex conjugate functions vp and & inducing the pair of dual coordinates and for Sn, with coordinate maps om and de such that
vCCy( [ (0, I.
(24)
Consider the gradient flow dynamics:
iit (Tt)),
7(0) Npu
0(t) ((()),
((0) 0pp-
Then for any > 0, there erist positive constants C12 02 C32 C₁ and T such that for all t 2 T,
-(-te)) < (n((t) Cge 2(c-s)) e3e 2(]++) < (,t((t) < C4 2(2-))
(25)
(26)
i.e., (T() and (()) converge exponentialyy with rate 2c and ? respectively.
11

---

Published in Transactions on Machine Learning Research (07/2025)
Proof. See Appendix C.4
Note that by plugging c= in Theorem see that there exists dual pair of coordinates that achieves the convergence rate of the natural gradient dynamics. However, the affine transformation that leads to this convergence rate depends on the target distribution q and thus cannot be known in advance. Furthermore, Theorem illustrates that the convergence rate of gradient flows in the transformed coordinates can be made arbitrarily small or arbitrarily large by scaling the coordinates. In contrast, the natural gradient flow is independent of the choice of coordinates, and therefore, has a coordinate independent convergence rate.
The superiority of the natural gradient method in terms of convergence rates is not immediately clear from the continuous-time analysis presented so far. In order to facilitate a meaningful discussion of convergence rates in continuous time, Muehlebach & Jordan (2020) propose a particular time-normalization approach that can be deployed in our setting. Alternatively, a more direct comparison between the Euclidean gradient and the natural gradient can be made by studying discrete-time gradient descent iterations. To claborate this, we turn our attention to the discrete-time setting in the next section, where the advantages of the natural gradient method become evident.
4 Convergence Analysis in Discrete Time
For a given target distribution g E Sn and an initial distribution po E Sn, the discrete-time gradient dynamics are given by
(k + 1) = (nkk= on (.g(((()) = (kk) + a. a. (()((())) - -(k)),
77()) = M/
(27)
(kk + 1) = (Kk) - Oxg (O(:)( = 0(k) - 0g ((())) + 02 (v))))))
((0) = 0 0p
(28)
Tgg (k: + 1) = mmgg(t ang grad (Tmg = mmg (k) ang (Thg (k) - na). Mag (0) = 1p
(29)
where Ony, 0g and Ong are the learning rates. Unlike in the continuous- setting, the choice of coordinates
used to represent the natural gradient dynamics in discrete time influences the analysis convergence rales,
primarily due to the presence of the learning rate a in the update equations (see Martens] (2020): [Song
et al. (2018)). In what follows, restrict the analysis to the representation of the natural gradient in the 77 coordinates.
To simplify the analysis, let us linearize these dynamics around the equilibrium points and examine the local convergence rates of the linearized dynamics. Owing to the already linear natural gradient dynamics, these do not need to be linearized. These linearized dynamics are given by
7)kh + 1) = (I - a, v2 ((a)) (1k) + a, p(ng)an (kk + 1) = (I - ag v (b₂)) 0(k) + V (0,)0
(7) = 1p
(30)
((0) = 6 p
(31)
Mag (k: + 1) = (1- (nnann)) (k) + Ong" 779
Ung (0) = 7po
(32)
It turns out that the discrete-time natural gradient dynamics in the coordinates, when linearized about the equilibrium 0g, leads to update equations that are identical to equation 32 This elaborated in Appendix
32 invariant to any affine transformation of the coordinates. Therefore, the update equation 32 represents local linearized dynamics for all dual pairs of coordinates.
B Furthermore, also note that the update equation
The dynamics described by equation 30 equation 31 and equation 32 can be written in the general form
x(A 11) (1 aQ) xkk) + +oQx*
where Q is symmetric positive definite matrix. Notice that these dynamics result from gradient descent iterations when optimizing the convex quadratic function f(z) = [xx "aa=-=""Tee =T= 10---- x²). In this discrote-time setting, we say that a sequence [(=(k) converges exponentially to f(x.) with rate p E [o, 1) if there exists a positive constant e and an integer ko such that [f(=(k)) - f(x))|5 ep holds for all 2 kg- smaller value of p corresponds to faster convergence. The convergence rates of gradient descent algorithms for minimizing convex quadratic functions have been extensively studied. For example, the following result from Nesterov 2018)) shows that the condition number of Q determines the convergence rates.
12

---

Published in Transactions on Machine Learning Research (07/2025)
Theorem 5 Nesterov (2018)). Let f(x) 2======"x===000000 Q(Q x²) with Q 0, and let denote the condition number of Q. Consider the gradient. descent. ilerations z(k + 1) = zzk) aV [(x())) Then:
(0) [x((k) converges to zero at rate (1-1)² when a (Q) (standard choice).
(i) f(((k)) converges to zero at rate (1
2 2 )² when a Amnis ( (Q) (optimal choice). x+1
Note that the gradient descent dynamics in and 0 coordinates correspond to setting Q = (M,,) and Q Q= ( respoctively. By directly applying these results to the discrcto-time gradient descent dynamics described by equation 30 and equation 31 we observe that poorer conditioning of Q leads to larger convergence rate P. and thus slower convergence. The condition numbers associated with the gradient descent dynamics in the and coordinates can be bounded away from as stated in Lemma 6
Lemma 6 (Bounds on the condition number of the Hessian). Let q E Sn. Them,
< Kig S Lo) (o,
(33)
where Kig
with Tmii and 7min being the smallest and the second-smallest element of [1ao1. . Ingn, respectively.
1imin.2
Vmni
Proof. See Appendix C.5
The natural gradient dynamics correspond to setting Q = I which yields optimal conditioning. Observe that the discrete-time natural gradient descent dynamics achieve a convergence rate of | - for a E (0, 2). Furthermore, achieves an optimal convergence rate 0, ie., convergence in a single step for the optimal learning rate a Note, however, that since this analysis pertains the linearized system, the actual natural gradient descent does not converge in single step. Finally, with the goal of studying the properties of the stochastic gradient descent (SGD), measurements. Although the noise models studied next do not exactly model the stochastic behavior of the SGD, they take us a step closer to it and provide valuable insight. Furthermore, practical implementations of the natural gradient method involve approximating the Fisher information matrix by an empirical version of it Martens 22200p. This can also be captured to some degree by the noise models studied next. Specifically, we study two noise models motivated by (Polyak] 1987 Chapter 4): relative deterministic noise (multiplicative) and absolute random noise (additive). These and other similar noise models have been studied in the optimization literature and they evidently show that the condition number plays central role in these
examine the robustness of these dynamics imperfect gradient
analyses (see Guille-Escures et al. 2021 h: Lessard et. al. 2016 Van Scoy & Lessard 2021
4.1 Robustness Analysis with Relative Deterministic Noise
Let us first consider the relative deterministic noise (multiplicative) model which replaces the gradient vector v by (I+ (bk)) where A (l)) E Ra captures the noise at time instant k. This leads to dynamics
(1k + 1) r7kk) a. (T A(+)) ((n() mg),
(7()) = 1pp
(34)
[ 1) 0(k) ax - - ((:) V (O,) 0g),
6(0) 0po
(35)
Thg (k: + 1) TTg (k) ang (1 (k))) ma). mag (0) T/pa
(36)
where the learning rates ans 0g and ang are chosen optimally assuming the noisc-free conditions (A = 0).
Theorem 7 (Robust stability under relative deterministic noise). Consider a largel distribution E Sa, an inittal distribution p E Sn and the discrete- dynamics described by equation 34 equation 33 and equation 30 respectively, where the learning rates Onp 0g and ang are chosen optimally for each case assuming the absence of noise (A =). Let k cond(V2ip(n)) cond(V?+((,))) Then the following statements hold:
13

---

Furthermore, for Published in Transactions on Machine Learning Research (07 (07/2025)
(i) If the sequence of perturbations ((k is such that for some > 0, [| < <1-- for all l 2 0, then the natural gradient dynamies described by equation 36 are stable, i.e., limg->oc Mng (k) 7/ 0.
(8) There erist time-tneariant perturbations An and As with ||n | 1.2. that destabilize the gradient descent dynamics described by equation and equation 35 respectively.
Proof. See Appendix C.6
The above result shows that the natural gradient dynamics exhibit larger robustness margin in comparison to the robustness margin of the y and 0 gradient dynamics which depend on the condition number This shows that the superiority of the natural gradient dynamics be again attributed to the optimal conditioning (K 1). Also note that the above noise model includes time-varying perturbations to the learning rate and shows that the natural gradient dynamics tolerate much higher deviation from the optimal learning rate, Furthermore, note that statement (i) of the above theorem proves convergence for the situation where the inverse of the Fisher information matrix G is replaced by (1+ -(+))G- with | < 1 for all & 2 0. This result thus also makes progress towards the more practical implementations of the natural gradient involving an empirical version of the Fisher information matrix Martens (2020).
4.2 Robustness Analysis with Absolute Random Noise
Now let us now consider the absolute random noise (additive) model which perturbs the original dynamics by adding an independent and identically distributed noise signal 6(k) for k E (0,1,...). This leads to dynamics
17Kk 1) mT(k) a, (ng) n4) + o(k), 0(k 1) 0(k) as (0,) 0,)) 8(k),
n70) Tpo
(37)
0(0) Opo
(38)
7ng (k 1) (k)) ang (nag (k) na) 5(k).
mag (0) Upp
(39)
where learning rates ag ag and ang are chosen optimally for each case assuming the absence of noise (8(k) = 0). Furthermore, assume that 8(k) an independent and identically distributed stochastic process satisfying B]](k)] 0 and E[ō(k)((k)*] = for all k 2
Theorem 8 (Robustness against additive noise). Consider a target distribution q E Sn, an initial distribution Po E Sn and the discrete-time dynasiic described by equation 37 equation 388 and equation 39 respectively, where the learning rates Ony ag and ang are chosen optimally for each case assuming the absence of noise (8(k) 0). Then
(i) oo B[(n(k) (n)((()))))))))))))))))))))))( - 11 = En,
(x+1)² dnc
(1) lim B[(((k)) (,)(O(k)) -0,) = Eg (x+1)² I
Ax
( B[(Ang(A) (b)) nyy (Tng (k) mT for all lo 2 0.
The upperbound in (i) and (ii) is tight, i.e. Mn and Eg have eigenvalues equal to (x+1)² n = 2, we get equality in (1) and (ii).
4x
Proof. See Appendix C.7
The above result explores the elfect of adding an independent and identically distributed (i.i.d.) noise signal at every iteration of the dynamics. It establishes that the largest eigenvalues of the steady-state error covariances are given by (x+1)² I for the 7 and 9 gradient dynamics whereas the error covariance with the 4x
natural gradient dynamics equals the noise covariance which corresponds to optimal conditioning (== 1).
14

---

Published in Transactions on Machine Learning Research (07 (075)
(kg)l((k))
0.6
0.2
Daglimng (k))
KL Divergence
KL Divergence
e (kglln(k))
D(a))((k))
Diglinng mag (k)) ((g))((k))
0.4
0.1
0.2
20
40
60
80
200
400
600
800
1,000
k
k
Figure 8: KL divergence evaluated along optimization trajectories generated by discrete-time gradient descent dynamics described by equation 40 equation 41 and equation 42 All methods use the same and sufficiently small learning rate. The left. panel shows results for n = 2 with Og = 0g = Ony = 0.01; the right panel shows n = 10 and On = ag = ang = 0.001. Convergence behavior is consistent with its contimous-time counterpart.
5 Numerical Experiments with Empirical KL Divergence
The gradient expressions analyzed in previous sections involve the target distribution q, which is typically unknown in practical machine learning applications. In reality, one only has access to a data sequence D (I1, NN) consisting of samples drawn i.i.d. with respect to q. This section aims to bridge the gap between our theoretical analysis and practical implementations by investigating empirical versions of the KL divergence and their associated optimization dynamics. Recall that the KL divergence can be decomposed
(p) = Σ 1=1 pi
qs log
(1 - Σ ai log (p:),
where (q) is the (negative) Shannon entropy of independent of and the second term is the tross-entropy equal to EEwg []og(P=)] As a result, minimizing (p) is equivalent to minimizing the cross-entropy term, since the entropy term is constant with respect to p. Given a data sequence D sampled with respect to g. we can estimate this expectation using the empirical mean
N log (px,)
n+1
Σ
q. log (p4),
where d the relative frequency of outcome in D. This defines the empirical target distribution (qi...-, (q.....,(q-1) 1n- Sn. This motivates the definition of the empirical KL divergence, defined for a given data sequence D. as (x(p) Lq(p). Since this differs from the empirical cross-entropy only by constant (independent of p), minimizing (p( is equivalent to minimizing the empirical estimate of (p). Therefore, our continuous-tine analysis in Section applies directly in this empirical setting by substituting q with i, assuming lies in the interior of the simplex S,, which typically holds for large enough data sequences. As N + 20, we have that q + q, ensuring that the empirical KL divergence converges to the true KL divergence. This provides a direct connection between our thcorctical results and their applicability in data-driven settings.
In what follows, we outline the setting described above more explicitly. Consider an arbitrary target dis- tribution and generate the data sequence D of independent and identically distributed samples drawn with respect to q- Analogous to the discrete-time updates given in equations equation 27 equation 28 and
15

---

Published in Transactions on Machine Learning Research (07 (07/202) /025)
convergenec time
60
convergence time
60
convergence time
60
40
40
40
20
20
20
0
0
2
3
0.5
1
1.5
5
10
15
20
Ong
.10
Ong
00
Figure 9: Convergence time as a function of learning rate for discrete-time gradient descent dynamics governed by equation 40 (left), equation 41 (right) and equation 42 (center). Results correspond to dimension n = 10. For each learning rate, convergence time defined as the maximum number iterations (across 100 random initializations) required to reach a fixed tolerance near the optimum. Natural gradient descent achieves the fastest convergence when appropriately tunod, as evidencod by comparing the lowest convergence times in each plot. These results are sunmarized in Table
equation 29 consider the empirical gradient descent dynamics described by
(/k: + 1) = n(k) - an vLD(n(k)),
(70) 1lp
(40)
0(k: 1) 0(k:)
(lD(O()))),
((0) pp
(41)
mng (I: + 1) ( - ang grad (1)))))))) "mg (0) 17p-
(42)
For sufficiently small and equal learning rates 0x1 continuous-time behavior. This is illustrated in Figure where we set an 0g = Ang = 0.01 for n = 2 (left) and Og = Okg 0lng = 0.001 for n = 10 (right). The sandwiching property derived in the continuous-time setting (see Theorem 3 thus extends, as expected, to the discrete-time setting when the learning rates are chosen be equal across the different methods and sufficiently small.
(nng the discrete-time dynamics approximate the
Drawing motivation from Theorem which analyzed convergence rates under different choices of learning rates, we next explore how varying the learning rates individually for each method influences empirical convergence properties. specifically, we sample 100 random initializations and, for cach algorithm and cach learning rate, run simulations from these 100 initializations. We then measure the convergence time defined as the maximum umber iterations (across the 100 initializations) required to reach a fixed tolerance from the optimuu. Figure plots convergence time as a function of learning rate for discrete-time gradient descent dynamics governed by equation 40 equation and equation 42 Results correspond dimension 10.
The optimal learning rates and convergence times shown in Table 1
Table 1: Optimal Learning Rates and Convergence Times from Figure
optimal learning rate optimal convergence time
17 coordinates natural gradient
On 0.0036
k = 29
Ang E [0.7141, 1.16]
k = 2
0 coordinates
Og E [1.99, 14.24)
k: = 12
With optimally tuned learning rates, natural gradient descent consistently converges faster than standard gradient descent in cither coordinate system.
5.1 Numerical Experiments with Stochastic Gradient Descent
Recall that the key motivation behind Theorems and R is to go beyond gradient descent potentially taking step towards better understanding the behavior of stochastic gradient descent (SGD). To this end, we extend our empirical study to the SGD setting, where at each iteration k, mini-batch Dz. D is
16

---

Published in Transactions on Machine Learning Research (07 (075)
(kg)l((k)
0.6
0.2
Dagl (k)) ((k))((k))
KL Divergence
KL Divergence
e (kglln(k))
Diglinng (k))))) ((g))((k))
0.4
0.1
0.2
20
40
60
80
0
200
400
600
800
1,000
Tteration
Iteration
Figure 10: KL divergence evaluated along optimization trajectories generated by discrete-time gradient descent dynamics described by equation 43 equation 44 and equation 45 All methods use the same and sufficiently small learning rate. The left panel shows results for n = 2 with On Ag Ang = 0.01; the right panel shows n = 10 and Og = ag = ang 0.001. Convergence behavior consistent with its contimuuss-time counterpart.
drawn uniformly at random, and gradients are computed with respect to the empirical loss LD. To satisfy standard conditions for the almost sure convergence of SGD, adopt time- time-varying learning rate of the a where the parameter a is chosen large enough to delay the onset of decay until sufficiently form a. k+aa
large i and is kept equal across all the different algorithms. The resulting dynamic equations are then described by
7/k + 1) = =(k)- an
VLo (n(k),
17(0) = Npp
(43)
k +
0(k: + 1) = 0(k)) ag
(L)x((())))
((0) = 0
(44)
k: k+ a
Mng (k + = (k) ang
grad (k)) Mng (0) = Tp
(45)
h k
Analogous to the full-batch setting, Figure 10 shows that for sufficiently small and equal learning rates On = og = Ongg, the sandwiching property derived in the continuous-time setting (see Theorem 3 extends to the setting of stochastic gradient descent. We next plot convergence times as a function of learning rate in Figure [1 and summarize the optimal values in Table 2 Once again, stochastic natural gradient descent
Table 2: Optimal Learning Rates and Convergence Times from Figure
optimal lcarning rate optimal convergence time
17 coordinates natural gradient 0 coordinates
ag 0.0025
k = 41
Ong E [0.541, 0.695]
k = 3
ag E 6.93
k = 20
outperforms its Euclidean counterparts when the learning rate is optimally tuned. Figure provides example trajectories illustrating the convergence behavior under both full-batch and stochastic updates with optimally chosen learning rates. Tn contrast. to sandwiching property observed in Figure g and Figure TO can be seen in Figure 12 that the natural gradient trajectorics converge faster than their Euclidean counterparts the full-batch well stochastic settings.
5.2 Summary
These empirical studies validate that our theoretical insights apply practical scenarios involving sampled data from the target distribution. In both full-batch and stochastic optimization settings, natural gradient
17

---

Published in Transactions on Machine Learning Research (07/202) 2225)
convergence time
100
convergence time
100
convergence time
100
80
80
80
60
60
60
40
40
40
20
20
20
0
0
2
4
5
0.5
1
1.5
2
5
10
15
-10 3
Ong
Okg
Figure Convergence time as function of learning rate for discrete-time stochastic gradient descent governed by equation 43 (left), equation 44 (right) and equation 45 (center). Results correspond to dimension 10. For each learning rate, convergence time is defined the maximum number of iterations (across 100 random initializations) required to reach a fixed tolerance near the optimum. Natural gradient deseent achieves the fastest convergence when appropriately tuned, as evidenced by comparing the lowest convergence times in each plot. These results are summarized in Table
- (k)l)n(k))
D(alln(k))
0.3
0.3
Diglling (lg))((k))
D(gllnng(() ((ql(((k))
KL Divergence
KL Divergence
0.2
0.2
0.1
0.1
10
15
20
20
40
60
80
100
Figure 12: KT. divergence evaluated along the trajectories generated by equation
40 equation and equa-
tion 42 (left) corresponding the full-batch setting and by equation 43 equation 44 and equation 45 (right) corresponding to the stochastic setting. Results correspond to dimension set to their optimal values obtained from Figure for the left panel and from Figure I for the right panel. Natural gradient descent exhibits faster convergence across both settings with appropriately tuned learning rates.
10 and the learning rales are
18

---

Published in Transactions on Machine Learning Research (07/2025)
descent (NGD) consistently outperforms standard gradient descent (GD) under optimally tuned learning rates. This corresponds to Theorems 5 and In contrast, when learning rates are chosen to be equal for all methods and sufficiently small, the sandwiching property established in Theorem 3 for contimuus--time dynamics holds also in the discrete-time setting.
Table 3 summarizes the correspondence between theoretical results and empirical findings. These empirical observations thus closely align with the theoretical results developed in the preceding sections, reinforcing the practical relevance of our analysis.
Table 3: Summary of Theoretical and Empirical Results
Sandwiching property
NGD outperforms GD
Setting
(a. = Ong Og sufficiently small)
(agg Ong: 0gg individually optimized)
Theory
Theorem
Theorems
Empirical studies (GD)
Figure
Table
Fi Figure
(left)
Empirical studies (SGD)
Figure
Table
F Figure
12 (right)
6 Conclusions and Outlook
In this work, we revisited the convergence properties of natural gradient flows in comparison to their Eu- clidean counterparis, focusing on the minimization of the KL divergence over discrete probability distri- butions. Our analysis revealed a more nuanced picture than the commonly held belief in the universal superiority of natural gradient methods. In the continuous-time sctting, while the natural gradient flow in- deed outperforms the Euclidean gradient flow in the 0 coordinates, consistent with traditional expectations, we showed that it converges more slowly than the n-gradient flow, despite following straight-line trajec- tories in these coordinates. This demonstrates that the commonly observed rapid convergence of natural gradient flow cannot be simplistically attributed to the straightness of its trajectories. Our discrete-time analysis of gradient descent dynamics further clarified that the fundamental reason behind the superiority of natural gradient methods lies in their optimal conditioning: natural gradient updates effectively minimize an optimally conditioned loss landscape, leading to consistently better performance compared to their Eu- clidean counterparts. Furthermore, our empirical results reinforce the theoretical results by demonstrating that natural gradient descent (NGD) maintains its convergence advantages over standard gradient descent (GD) even in practical, sample-based settings. In both full-batch and stochastic optimization, NGD ex- hibits faster convergence when learning rates are optimally luned, and satisfies the theoretical sandwiching property under sufficienlly small and equal learning rates. These findings highlight the practical relevance of the theoretical analyses and support the use of NGD in rcal-world applications where access to the full distribution is limited to finite samples. Overall, our findings refine the understanding of natural gradient methods and highlight the subtle, yet important, nuances that govern their behavior.
The theoretical results presented in this paper primarily concern exact natural gradient, with robustness analyses taking initial steps toward incorporating noise in gradient measurements and system dynamics. A natural next step is to extend this analysis to practical settings where the Fisher information matrix must be approximated. For instance, Theorem 7 already accommodates structured perturbations of the form (I + A(k))G-1, where G is the Fisher information matrix, provided [|A())]a < 11- E for all k 2 O, ensuring the stability of the resulting updates. This suggests that if approximations such as K-FAC can be modcled as structured perturbations of G-1 satisfying the above condition, then it may be possible to establish convergence guarantees for the resulting approximate natural gradient algorithms. A systematic exploration of the trade-off between the quality of the approximation of the Fisher matrix and the resulting convergence guarantees would be an important direction for bridging theory and practice.
While we focused on the probability simplex equipped with dual coordinate systems, an important extension would be to general dually flat statistical manifolds (Amari & Nagaoka (2000); [Ay et al. (2017)) ) where similar tools may be employed to study optimization dynamics in broader settings. Such manifolds are induced by general Bregman divergences going beyond the setting of KL divergence. Notably, the identity (2p((n) [V2 (0))- holds for any pair of dual Bregman divergences induced by strictly convex functions y
19

---

holds at the optimum. One can then proceed to derive the Published in Transactions on Machine Learning Research (07/2025)
and b, respectively Amari & Nagaoka (2000] So it is plausible that one can derive analogue of Lemma for this general setting provided (1)(()) >
"sandwich" ordering as in Theorem
Additionally, our framework may be extended to richer families of probability distributions, such as general exponential families derived from Boltzmann machines without hidden units, and more intricate mixtures of exponential families associated with Boltzmann machines with hidden variables (Amari et al. (1992)) These models exhibit more complex geometries that may reveal deeper interactions between parametrizations and optimization dynamics.
While our analysis in Section relies on linearization of the gradient flows around the optimal solution (i.e., the larget distribution 4), we acknowledge that extending these results to the fully nonlinear setting remains an interesting direction for futurr work. A potential avenue for future research involves characterizing a neighborhood around the optimum in which the Hessians of the loss functions can be uniformly bounded, m L. This would allow from optimization theory for the class S(m, L) of strongly convex functions with Lipschitz-contine dients. In such settings, it may be possible to rigorously extend convergence guarantees and comparative analyses to the nonlinear regime. More generally, for any convergent optimization algorithm, the nonlinear dynamics necessarily approach the linearized ones asymptotically near the optimum. Thus, for any con- vergent algorithm, the trajectory of the nonlinear dynamics must eventually enter a neighborhood of the optimum where the linear approximation becomes accurate. Therefore, our results can also be interpreted as describing either the local or asymptotic behavior of the optimization dynamics.
i.e,, their eigenvalues lie within [m, L] for some
draw on classical results
gra-
Finally, although discreto-time analysis highlights robustness advantages of natural gradient methods, it does not fully capture the stochasticily inherent in stochastic gradient descent (SGD). Developing more precise theoretical model that explicitly incorporates the stochastic dynamics of SGD remains an important avenue for future work.
References
Shun-ichi Amari. Neural Learning in Structured Parameter Spaces Natural Riemannian Gradient. In Advances in Newral Information Processing Systems, volume 9. MIT Press, 1996.
Shun-ichi Amari. Natural Gradient Works Efficiently in Learning. Neurol Computation, 10(2):251-276, 1998.
ISSN 0899-7667. doi: 10.116200889776003300077777
Shun-ichi Amari. Information Geometry and Its Applications, volume 194 of Applied Mathematical Sciences. Springer Japan, Tokyo, 2016. TSBN 978-4-431-55977-1 07171-555777717--44----557777777 doi: 10.1007/988-44431-66777-
Shun-ichi Amari and Hiroshi Nagaoka. Methods of Information Geometry. American Mathematical Soc., 2000. ISBN 18--0-8218-400-44
Shun-ichi Amari, Koji Kurata, and Hiroshi Nagaoka. Information Geometry of Boltzmann Machines. IEEE Transactions on Neural Networks, 3(2):260 271, 1992. ISSN 1941-0093. doi: 10. 10.1109/72.122861
Nihat Ay. On the Locality of the Natural Gradient for Learning in Deep Bayesian Networks. Information Geometry, 6(1):1 49, 2023. ISSN 2511-249X doi:
020-00038-y.
Nihat Ay and Jesse van Oostrum. On the Risher-Rao Gradient of the Evidence Lower Bound. arXiv preprint arXiv:2307.11249, 2023. URL https: httpss://axxvvvoroee wttp:///aabb//..112222222 11249
Nihat Ay. Jürgen Jost, Hông Vân Lé, and Lorenz Schwachhöfer. Information Geometry, volume 64 of Ergebnisse Der Mathematik Und Ihrer Grenzgebiete 34. Springer International Publishing, Cham, 2017. TSBN 978-3-319-564777 378-3-319-56477-44 doi: 1110079978-33319-566
Rajendra Bhatia Perturbation Bounds for Malrix Eigenwalues. SIAM, 2007.
Nicolas Boumal. An Introduction Optimization on Smooth Manifolds. Cambridge University Press, first edition, 2023. ISBN 978- 978-1-00-916616-- 978-1-00-916617111 978-1--0--116615-7
doi: 51011000111000777777771000101110
20

---

Published in Transactions on Machine Learning Research (07/2025)
Charles Hallle- Manuela Girotti, Baptiste Goujaud, and Ioannis Mitliagkas. A Study of Condition Numbers for First-Order Optimization. In Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, pp.
PMLR, 2021.
Sham M Kakade. A Natural Policy Gradient. In Advances in Neural Information Processing Systems, volume 14. MIT Press, 2001.
Ryo Karakida and Kazuki Osawa. Understanding Approximate Fisher Information for Fast Convergence of Natural Gradient Descent in Wide Neural Networks. Journal of Statistical Mechanics: Theory and Experiment, 2021 (12):120000 December 2021. ISSN 1742-5468. doi: 10.1088/1742-55 0555//77722-444//3///5 5416///33////..... /ac3ae3.
Ryo Karakida, Shotaro Akaho, and Shun-ichi Amari. Universal Statistics of Fisher Information in Deep Neural Networks: Mean Field Approach. In Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics, pp. 1032-1041. PMLR, 2019.
Sajad Khodadadian, Prakirt Raj Jhunjhunwala, Sushil Mahavir Varma, and Siva Theja Maguluri. On Lincar and Super-Linear Convergence of Natural Policy Gradient Algorithm. Systems & Control Letters, 164: 105214, 2022. ISSN 01676911. doi: 10.1016/j.syscoole.. 202..105211..
Laurent Lessard, Benjamin Recht, and Andrew Packard. Analysis and Design of Optimization Algorithms via Integral Quadratic Constraints. STAM Journal on Oplimization, 26(1):57-95, 2016. ISSN 1052-6234, 1095-7189. doi: 10.1137/15M100990
James Martens. New Insights and Perspectives the Natural Gradient Method. Journal of Machine Learning Research, 21146)):1-76, 2020.
Mehryar Mohri, Afshin Rostamizadeh, and Ameet Talwalkar. Foundations of Machine Learning. Adaptive Computation and Machine Learning. The MIT Press, Cambridge, Massachusetts, second edition, 2018. ISBN 978-0-262-03344--6.
Michael Muehlebach and Michael Jordan. Continuous- Timei Lower Bounds for Gradient-Based Algorithms. In Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pp. 7088-7096. PMLR, 2020.
Johannes Müller and Guido Montúfar. Gcometry and Convergence of Natural Policy Gradient Methods. Information Geometry, 7(1):485 523, 2024. ISSN 2511-249X. doi: 1011077/418884-22-0000---
Yurii Nestcrov. Lectures Conver Optimization, volume 137 of Springer Optimization and Its Applications. Springer International Publishing, Cham, 2018. ISBN 313-311-9117777777-----------109
doi: 10.1007/
978-3-319-91578-4
Razvan Pascanu and Yoshua Bengio. Revisiting Natural Gradient for Deep Networks, 2014.
Boris T. Polyak. Introduction to Optimization. Translations Series in Mathematics and Engineering, 1987.
Yang Song, Jiaming Song, and Stefano Ermon. Accelerating Natural Gradient with Higher-Order Invariance. In Proceedings of the 35th International Conference Machine Learning, pp. 4713 4722. PMLR, 2018.
Eduardo D. Sontag. Mathematical Control Theory: Deterministic Finite Dimensional Systems. Springer Science Business Media, 2013. ISBN 978-1-4612-0577-7.
Jesse van Oostrum and Nihat Ay. Parametrisation Independence of the Natural Gradient in Overparametrised Systems. In International Conference on Geometric Science of Information, pp. 726-735. Springer, 2021.
Jesse van Oostrum, Johannes Müller, and Nihat Ay. Invariance Properties of the Natural Gradient in Overparametrised Systems. Information Geometry, 6(1):51 67, 2023. ISSN 2511-249X. doi: 10.1007/ 4418844002-00066..
Bryan Van Scoy and Laurent Lessard. The SpeeddRRbustae Trade-off for First-Order Methods with Addi- tive Gradient Noise. Xiu preprint arXiv:2109.05059 2021. URL htttp::/arxiv..et wwwp:///ww.wr.ab/b//////// 05059.
21

---

Published in Transactions on Machine Learning Research (07 (025)
Patrick M. Wensing and Jean-Jacques Slotine. Beyond Convexity Contraction and Global Convergence of Gradient Descent. PLOS ONE, 15(8)=e0236611, 2020. ISSN 1932-6203. doi: 10.137 (0.1371/jouna)) /oso///aaall.poe.022200 0236661..
Lin Xiao. On the Convergence Rates of Policy Gradient Methods. Journal of Machine Learning Research, 23(282):1 36, 2022. ISSN 1533-7928.
Rui Yuan, Simon Du, Robert M Gower, Alessandro Lazaric, and Lin Xiao. Linear Convergence of Natural Policy Gradient Methods with Log-Linear Policies. In The Fleventh International Conference Learning Representations, Kigali, Rwanda, 2022. URL https:: fttt::://teeecompPeee parioom--arii..... hal. sclence/baal-0418266
Guodong Zhang, James Martens, and Roger B Grosse. Fast Convergence of Natural Gradient Descent for Prer-Parameterized Neural Networks. In Advances in Newrol Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
A
Gradient Flows for D
For a given target distribution p E Sn and an initial distribution qu E Sn, consider the gradient flow dynamics described by equation 46 and equation 47 and the natural gradient flow dynamics described by equation 48
given by
fi(t) vCz(((t)
n(0) 77
(46)
((t) (O)((()),
8(0) 0qp
(47)
e (t) grad (,(0ng(t)) = (ng + 8p. Ong (0) = 6
(48)
The following is an analogue Theorem 3 applied to the above dynamics.
Theorem 9 (Convergence analysis). Let p E Sn and E Sn such that ((qullp) is sufficiently small.
Suppose 4, 0 and Ong be the solutions to dynamics described by equation 40 equation 47 and equalion 48
respectively. Then
(i) there erist positive constants m² S L3 cg and ezsuch that
cee S cGE 24 < ()(()) S che 2m3
V8 1200 0,
(49)
i.e., ()(() converges exponentially with rate lower than 2.
(i) there erist positive constants
< m. S Li and en such that Ge - 24 < ((((() cGe G 2m0t co.e
V4 0
(50)
i.e., (n(tt) converges exponentially with rate higher than 2.
(il)) there erist positive constants ci, ca and T such that
ce < (Ong (t)) S cse -2 Vt V20
(51)
i.e., (t)g converges exponentially with rate 2.
Analogous to the proof of Theorem 3 (see equation [6 and the discussion thereafter), it can be shown using the Pinsker inequality (Mohri et al. 2018 Proposition E.7) that np is the unique minimizer of 2 and S := {n E om (S,)),,4 - (n()) is compact. Furthermore, compactness of S along with bound
Proof. Consider the gradient flow dynamics described by equation 46
15 implies that
I m
I 4 (m) = L
En S;
where m. min Amin (V (V'C,()) and L; max Amax (V²C;())). Applying Proposition we
get the desired inequality given in equation 50
22

---

Published in Transactions on Machine Learning Research (07/2025)
Now consider dynamics described by equation 47 implies that there exists an 0 such that
Continuity of C along with equation 18 from Lemma
< m. 1 13 C, 3 L. 1
40 E (0p)
(52)
for some positive constants mă _ Lg 1. Using this local strong convexity condition, it can be shown
that if (010) is sufficiently small, S == {0 (s,)) < (,()())) is contained in (0p). Applying
Proposition 1 we get the desired inequality given in equation 49
Finally consider the gradient flow dynamics described by equation 48
which can be solved exactly to obtain
ng (t) 0g + e (0o - 0g)
(53)
Since lim Ong (1) = Og: we can use equation 52 along with (Nesterov 2018 2.1.8) to show that if (0(0)) is sufficiently small, then
Theorem 2.1.5 and Theorem
Lg m 10ng (t) ngll3 ()(ong(t) S [0n (t) 0yl1
(54)
for t 2 0. Plugging in the exact solution from equation 53 get the desired inequality given in equation 51
B Linearized Discrete-time Natural Gradient Dynamics in 0 Coordinates
In this appendix, we show that the discreto-time natural gradient dynamics in the coordinates given by (ng(k + 1) Ong (k) Ong grad L. (k))
0ng(0) Opo
when linearized about the equilibrium 0g. lead to update equations that are identical to the ones in the n coordinates. From the defining property given in equation of the natural gradient, and equation 6 we get that
grad - -10 (V)
(0,)) (0 0g),
where denotes a first-order approximation obtained by linearizing around 0g. The linearized dynamics the 0 coordinates are thus described by
0ng(k: (k + 1) (k) ang (Ong (1:) 0g).
0n (0) 8pg
which has the same form as in in equation 32
C Proofs
C.1 Proof of Proposition
Proof. [Proposition
Let x be a solution of equation 13 and define E(t) = f(x(t)) - f(x.). Note that I
B(() = (VS(x(()),
- (V)(((()))())))(())))))
= ([V]/((((()1 < 0.
(55)
Therefore, for all I 1 lo, [(x(I) (x.) E() < E((o) = s(x(lo)) - f(x.) which implies statement (i).
Furthermore, we get from Nesterov 2018 Section 2.1) that m 1 1 f2 = L- for all x E S implies
2m(f(x) (((())))(xx)))>)))) < [[
2l(((x) fxx.))
Vx E S.
(56)
Since z(t) S for all t 12 to, inequalities given in equation 56 and equation 55 imply that for all 2 to.
-2 F(() F(()) < -2m F(())
23

---

Published in Transactions on Machine Learning Research (07 (075)
Integrating from to to get that
E(t) = E(to) / E(s)ds s E(to) + [ (-2m) E(s)ds, to
E(t E(to) + f( -E(s) -E(to) ' (-2))-(-EEss)d
Finally applying the Bellman-Gronwal. Lemma (Sontag 2013 gives us
Lemma C.3.1) to the above two inequalities
E(t) - E((p)e 2m(t-ta)
-E(t) S c 2L(t-ta)
which diretly gives us the desired inequality 14
C.2 Proof of Lemma 2
Proof. [Lemma 2 The Hessians of Lu can be evaluated using equation 5 and equation 6 as
v2l,()) = v2 - Di)p(n)) nl. (0) (2)(
(57)
where D2(((n) R²n + Raxn is the third order derivative of ip whose action on a vector R" is given by (D&p(n)(e)),s Ek (()()) on.onome Uk In particular, note that the inverse relationship given in equation give
V (n)) = (rip(m)) IV (0,)) [,
(58)
Similarly, we can compute the Hessians of C using equation and equation as
(² = v²y(1))
(59)
V² = <²v(()) 12 0)
and use equation to obtain the inverse relationship
v (,,,)) ((p(pp))) [V' [V3
(60)
Recall that
p()) = Σ log Mi
+
1
Σ
n
log
-
Σ
==1
1=1
and it's Hessian V2 Φ can be explicitly computed to be
1
1
Vpp(n)
(61)
2i=1
7n
Since the second matrix the right hand side is positive semi-definite, get
14 max, 7h 1 1 V2
(62)
This together with equation 59 gives us
I 12V2
Vn E (Sm)
(63)
24

---

Published in Transactions on Machine Learning Research (07 (07/2025)
and together with equation 57 and the inverse relationship given in equation gives us
0 V (0) -
v0 (S.)
(64)
proving equation 15
Furthermore, it can be shown by direct computation that
1 mi
(1-
1- 11 Indls )
V (7)
+ (1-
((1- Σイ mi)
-0
vn E @m))
(65)
proving equation L6
Finally, ovaluating the global bounds from equation 15] at the optimum points and using the inverse rela- tionships given in equation and equation 18
58 and equation 60 yields the desired local inequalities described by equation 17
C.3 Proof of Theorem 3
Proof. [Theorem First consider the gradient flow dynamics described by equation 11 and note that (gg) = O. Using the Pinsker's incquality (Mohri et al. 2018 Proposition E.7) along with the fact that @n o 1 is bijective, we have that for any Op 7 0gg L0p0p) 0. Thus 0, is the unique minimizer of b, giving condition a) of Proposition It can be shown that Cg has bounded sublevel sets in the I coordinates. To this, observe that |[]]| + oo implies that for some i, 10:] + 00, i.e., either 0; -00 or O, co. This implies that
et Σラ= +0
Or
Pat1 = 1 Σソ= ea 7 0.
Since the Lg blows up to infinity the boundary of the simplex, get that Lg is coercive, i.e, L, p) + 00 as ||@|| + co. Observe that if some sublevel set of Ly is unbounded, there exists a sequence of points (2)) ... such that [[0()]] - oo but (0))) stays bounded. This a contradiction. Therefore, we have that Lg has bounded sublevel sets in the coordinates. Finally, continuity of Ey implics that the sublevel sets are closed which implies that Sg is compact. Using compactness of S. = (0 E (S,)) (0(0)) and the global bound given in equation 15 from Lemma 2
we get that
< m4 1 Lo I
V0 Sa
where my minees, Amin (0) 0 and Le maxess, Amax ()), < 1. Applying Proposition gives us the desired conclusion in the form of equation 20
Analogous to the previous case, now consider the gradient flow dynamics described by equation 10 We can use Pinsker's incquality (Mohri et al. 2018 Proposition E.7) to obtain
0111-
Inall ))
(y(
2
1
m Ingli
+
+
1
Σ
T -(1-1
2
==1
2
1
1
1 2
Σ n [nglil
Ing - n/li
|17g 1//22
(66)
2
This shows that By is the unique minimizer of by and together with the continuity of L implics that
S = {n E dmm (S,)) < (x(()))
is compact. Using the global bound given in equation 16 we get that
0 < m, 1 - v2 = In
Vn E Sn
(67)
25

---

Published in Transactions on Machine Learning Research (07/2025)
where myy ss,, Amin (V2Lg(n)) 0 and L, maxnes, Amax (V)L)()))) Applying Proposition see that there exists a positive constant such that
ce 2f & Ln((t)) I ce 2m,
v2 0.
(68)
This implies that (n)(t) converges to O. Thus, for any > 0, there exists a T > 0 such that (n((t)) - 25 for all 2 T, which, using the lower bound obtained in equation 66 implies that || - 7((1)l - for all 1 2 7.
This means that for any E > O, there exists a T > 0 such that the set ST := {n E om (S.) C, (1) _
(,(((T))) (ng). Furthermore, continuity of V Eg along with equation 17 from Lemma 2 implies that there exists > 0 such that
I K myy I 13 V Im
Vn E E (1g)
for some positive constants < m, I Lg. Putting everything together, we get that there exists a T 2 0 such that X myy I (n) I Ln holds for all E Sr. Thus, applying Proposition with to = T, we get that there exists a positive constant such that
( 2bn = že -L, T) S (,(()) & ōe -mm = (mm) de 2m,
holds for all f 2 T which is the desired form in equation 19 with Cn 20, and E, = 2m, Finally, note
that if Eg (n7()) is sufficiently small, then So c (1g) implying that the above bound holds with T = 0. This completes the proof for statement (i).
Finally consider the gradient flow dynamics described by equation which can be solved exactly to obtain
Bag (f) = 77 + e (Nho 7n)
(69)
Using
Nesterov 2018
Theorem 2.1.5 and Theorem 2.8 along with equation 67 we get that
L
my 2
Mma (t) null2 S (Yig (t)
[1 (t) nyl12
(70)
2
Plugging in the exact solution from equation 69 we get the desired incquality given in equation 21
C.4 Proof of Theorem
Proof. [Theorem Since ())((,,) is symmetric positive definite, we can consider its symmotric matrix square root Dy. such that (),) = Dg Dg- Setting A= = Ve-D, in equation 22 and equation 23 and using we obtain equation 24 Using the continuity of Hessians, we get that for any e O, there exists d > O such that
the inverse relationship (2 = [22 (0,)),
(c (c- 8) (T) 1 (c
Vi E (s
(-- . C = ( +
+
.I
vB E B3 (62))
Analogous to the proof Theorem 3 we can apply Proposition 1
to obtain the desired incqualitics given in
equation 25 and equation 26
C.5 Proof of Lemma
Proof. [Lemma
6
Recall that
(n = =²p(())= -V' = =VV2,,00,)1~
This establishes the equality cond 2 (ma)) condd (Q,))
Recall that
7
v²v(n)=
+
(71)
26

---

Published in Transactions on Machine Learning Research (07 (07/2025)
Applying Weyl's inequalities Bhatia
2007 to the above rank one perturbation matrix, get
Amax
(n))
and
Amin (V²((n))
Muiin
1/min,2
This diroctly implics the desired inequality given in equation 35
C.6
Proof of Theorem 7
Proof.
[Theorem
Consider perturbed dynamics of the form x(k+ + 1) <<k) a(I + (z))) =²),,
(72)
where Q is symmetric positive definite matrix. This encompasses the dynamics described by equation 34 equation 35 and equation 36 by setting Q equal to V) (2()(f,) and 1, respectively and sctting x equal to 77.. 0g and f7g respectively.
Let A We will now prove that the dynamics described by equation 72 are stable, i.e.,
Amin
limg-+so (k(kk) - -=1== O, if for some & > 0, ||a())]22 < - for all k 2 0. This would directly imply
statement (i) by plugging in Q I.
By defining the error variable e(k) = r(k) x,, we get that e(k + 1) (I aQ)e(k) (A(K))e(A))
(73)
2 Note that a is chosen optimally assuming no noise (A = 0), iee, (Q) (see Theorem 5 With this choice of a, we get that - ||/-a = and ||aQ||2 x+I where the ||-1|2 is the induced 2-norm which coincides in our case to the largest eigenvalue magnitude owing to symmetry. This can scen most directly by an cigenvalue decomposition of the involved matrices. Therefore, using the triangle inequality and the submultiplicative rule of induced matrix norms, we get that,
2x
2x R x+1) [|e(A>]|
(1 + 2/ (1)))2) K |e(k)]l-
[le(k + 1)ll s
+ + ||())]
+ 1
Note that if [1())]2 < - for all #20 0, then (++)])(()]l))<=] of inequalities
K K +
< which gives us the chain
x+1
|e(k + 1)ll s pl|e({)]] S p²l|e(k - 1)ll ... ***te(0))
Since p < we get that limg le (k)]]=== 0. Since Q implies 1, this implies statement (i).
Since the above argument only derives a sufficient condition for stability, we still need to prove statement (ii) separately. We now construct a time-invariant perturbation A such that || = and the dynamics described by equation 72 are unstable, i.e,, =(k) does not converge to x². To this end, let Q = UAUT be the cigenvalue decomposition of the symmotric matrix Q where is the diagonal matrix containing eigenvalues A1 A2 5 ... < An in ascending order along the diagonal. Construct A as
0 0
A = U
UT.
Plugging this in equation 73 and using
Ar+An we get that
e(k: + 1) = Me(k),
(74)
where M (1 A1+ (1 + A)Q) contains an eigenvalue at -1. Since stability of dynamics described by equation 74 requires the spectral radius of M be less than 1, and since M contains an cigenvalue -1, the dynamics are unstable. This completes the proof for statement (it) by applying the constructed the choices (2y(1)) and Q (2b((0)) respectively.
27

---

Published in Transactions on Machine Learning Research (07 (025)
C.7 Proof of Theorem 8
Proof. [Theorem of the form
8
Following the same stralegy as in the proof of Theorem consider perlurbed dynamics
1)) (Q(((k)) x²) + o(k),
(75)
where Q is symmetric positive definite matrix. This encompassses the dynamics described by equation 37 equation 38 and equation 39 by setting Q equal to vln(n,, (2v)(0)) and I, respectively and setting x equal to Ml. 0g and f/g respectively. By defining the error variable e(k) == x(k) - xx, we get that
e(k 1) (K) 8(k).
(76)
Note that is chosen optimally assuming noise (6(+) = 0), i.e., x + Aman (Q) (see Theorem
k-1 With this choice of a, get that ||/ 1 K-I. where the 11-112 is the induced 2 norm which coincides in our case the largest magnitude eigenvalue owing to symmetry. Define P(k) == (le(k)e((k)) where the expectation taken over the different realizations of the noise process 8(k). Using the fact that ō(k) and e(k) are independent random variables along with [[((k)] O and ([Ō(k)((k))] I, get that
P(k 1) oQ) -Q) +
(77)
Since aQ) has all eigenvalues in (-1, 1), it can be shown that limg P(k) where solves the
P = (I(a +1.
(78)
To solve this equation, let Q = UAUT be the eigenvalue decomposition the symmetric matrix Q where A is the diagonal matrix containing eigenvalues è S A2 S S An in ascending order along the diagonal.
Note that equation 78 can be solved obtain
I-#
P=U
UT
1-pa
where H4 are the cigenvalues of (I aQ). Therefore, the largest eigenvalue of P is 1-(
(x+1)²
This proves statements (i) and (ii) by plugging in equal to (²re(n)) and (2y) respectively. Finally, plugging Q and Amnn Amace = 1 in equation 77 directly gives statement (iii).
28


Reviews:
[{"review_id": "GWHJ8IyTFe", "reviewer": "Reviewer_DB1o", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "First and foremost, I would like to express my sincere appreciation for the authors’ efforts in developing this paper.\nI would also like to note that the following review is written with the understanding that TMLR places a strong emphasis on accuracy, persuasiveness, and clear evidence, as well as capturing the interest of readers, rather than solely focusing on novelty or impact.\n\n\n## Strengths\n- **Rigorous, transparent theory:** The paper derives exponential-rate bounds for the three flows and proves the “sandwich” ordering in Theorem 3, resting on uniform Hessian eigenvalue bounds (Lemma 2).\n- **Fresh coordinate-invariance perspective:** The affine-reparameterisation result (Theorem 4) shows how Euclidean GD can be made arbitrarily faster or slower while NGD stays fixed at rate 2, clarifying that NGD’s edge is invariance, not intrinsic speed.\n- **Discrete-time insight and robustness:** By linking each method to a quadratic with Hessian condition number $κ$, the authors show NGD enjoys $κ = 1$ and therefore the best linear rate and the largest stability margin; Lemma 6 and Theorem 7 back this up for multiplicative and additive noise.\n- **Empirical corroboration:** Trajectories and log-slope fits for $n = 2$ mixtures match the analytic predictions almost exactly (slopes ≈ 7 / 2.0 / 0.5 for η-GD/NGD/θ-GD in Figure 3.), giving the theory convincing face validity.\n- **Clear presentation:** The exposition is largely self-contained, with helpful figures, consistent notation and full proofs in the appendix, making the paper accessible to both optimisation and information-geometry audiences.\n\n\n## Weaknesses\n- **Limited empirical scope:** All experiments use $2$-dimensional simplices; it remains unclear whether the continuous-time ranking or the $κ = 1$ advantage of NGD persists in higher-dimensional or real-world models.\n- **Theoretical evaluation relying on single metric:** Results are specific to KL divergence (I understand that this should place as a future work); extending the analysis to other Bregman or $f$-divergences—or even squared loss on the simplex—would will push the better understanding of convergence of NGD.\n- **Limitation of noise model analyses:** In robustness analysis, the authors assume a GD scenario with either deterministic relative noise or i.i.d. additive noise. Therefore, the theoretical results presented in this paper do not yet extend to more realistic and commonly used settings—such as those involving stochastic gradients—where both noise types may coexist and exhibit state-dependent variance. Extending the analysis to such contexts remains an important direction for future work.\n- **Limitation of theoretical understanding:** The theoretical results presented in this paper do not address the convergence behavior when using practical natural gradient approximations such as K-FAC, nor do they provide detailed analysis for alternative geometry-aware optimization methods, such as mirror descent. This could be also an important direction for future work.\n- **Typos:** A few minor issues, e.g., “discrete-tine” → “time” (in Section 3).\n\n## Concerns regarding the interest of TMLR readers\n- While the paper delivers a neat information-geometric analysis, its scope is narrowly focused on KL divergence over the probability simplex and on the natural-vs-Euclidean gradient debate.  Many TMLR readers—especially those working on large-scale model or applications—may not immediately perceive actionable take-aways because the toy-scale examples stop at $n=2$ components. Of course, using the setting $n = 2$ is entirely reasonable for verifying the theoretical results. However, including results for multiple values of n would make the findings more intuitive and enhance their credibility.\n\nThese points aside, the paper offers a crisp theoretical contribution with a nuanced take on when natural gradients help, and with a bit more empirical breadth and practical discussion it could be a strong addition to TMLR.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper revisits natural-gradient descent (NGD) for minimising KL divergence on the simplex and delivers four take-aways:\n- Continuous-time theory: NGD’s convergence rate is sandwiched between two Euclidean gradient flows—slower than mixture-coordinate GD but faster than exponential-coordinate GD—so it is not universally optimal.\n- Coordination insight: An affine re-parameterisation can always make a Euclidean gradient flow match NGD’s rate, underscoring that NGD’s edge lies in its coordinate invariance, not raw speed.\n- Discrete-time advantage: Once discretised, NGD enjoys perfect conditioning (Hessian condition number = 1), giving faster linear convergence and greater robustness to gradient noise than the Euclidean updates.\n- Empirical confirmation: Simple 2-component mixture experiments align closely with the analytic rates and visualise why NGD’s straight geodesic paths are sometimes slower in continuous time yet superior in discrete settings.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "The following list is what I believe it is important to clarify and correct. If I have misunderstood anything, I would greatly appreciate your guidance in identifying the inaccuracies.\n\n- Extend the empirical study beyond $n=2$.\n  - Please include experiments on at least one higher-dimensional simplex (e.g., $n=5 \\ \\textrm{or} \\ 10$) and report the observed linear-rate constants and stability margins, and discuss whether the “η-GD ≻ NGD ≻ θ-GD” ordering and the $κ = 1$ conditioning advantage of NGD still hold.\n\n- More detailed limitation of this work and future direction\n  - The manuscript would benefit from a more detailed discussion of the applicability and limitations of the presented theoretical findings, as well as important future directions. For example, it would be valuable to clarify that the current analysis does not account for practical approximations such as K-FAC, and to speculate on how such approximations might affect the convergence behavior or robustness properties. Such a discussion would help a broader TMLR readership better understand the scope and implications of the work.\n  - Briefly comment on other Bregman or f-divergences: would the sandwich ordering or $κ = 1$ property plausibly generalise? Even a speculative remark with citations would help position the work.\n\n- Correction of Typos.\n  - Please carefully check for typographical errors throughout the paper (e.g., Section 3).", "broader_impact_concerns": "I believe this work does not raise any ethical concerns, as it is a theoretical study focused on the fundamental issue of convergence of (N)GD in the context of machine learning.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "xJg8wR5ZAl", "reviewer": "Reviewer_rhRm", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths: 1. The research problem is fundamental in nature. While one can argue that NGD is natural, its convergence properties are less clear. 2. The paper provided interesting and rigorous observations in those regards.\n\nWeaknesses: 1. The paper focuses on learning a Categorical distribution. While it is an interesting scenario, it is to some extent limited, and it is unclear whether the results apply for more general scenarios.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "From the viewpoint of information geometry, Categorical distributions can be parameterized using either the exponential family parameterization or the mixture family parameterization. The paper studies the convergence properties of gradient descent (GD) under these two different parameterizations and natural gradient descent (NGD). It was shown that, in the continuous time setting NGD not necessarily outperforms GD; instead, its convergence rate may lie between two parameterizations of GD. However, under discrete time setting, NGD performs better and yields more robust updates.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Can the authors comment on the applicability of the analysis to more complex cases, e.g. for distributions that are not in the exponential family or when training a non-linear neural network?\n\nMinor: \n\nPage 9, third to last time of the last paragraph of Section 3: \"discrete-tine\" should be \"discrete-time\"", "broader_impact_concerns": "Not applicable as the paper is theoretical", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "5x5Z5BtDiT", "reviewer": "Reviewer_Ccrj", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strengths:\n\n- The paper provides an interesting insight that the good performance of the natural gradient descent cannot easily be attributed to its straight trajectories and Euclidean gradient flows can be faster than natural gradient flow if the coordinate is well-chosen for the problem at hand. This is expected to have an impact and open up a new research direction: given an objective function, design coordinates so that Euclidean gradient flow is better than the natural gradient flow.\n\n- The writing is of high quality, and the mathematical exposition is generally precise. Illustrative figures are provided along the way to back the claims, which I found quite pleasant to follow.\n\n## Weaknesses:\n\n- The faster convergence of the Euclidean gradient flow only exhibits after some positive time T that we do not have any quantitative knowledge about. This makes the conclusion \"the Euclidean gradient flow w.r.t. the mixture coordinate is faster than the natural gradient flow\" a bit qualitative, i.e., we do not really know when it is actually faster (can be at the very end of the flow).\n\n- The paper lacks a discussion of how the studied gradient flows relate to machine learning contexts. For example, in Section 4, applying these flows in practice requires full knowledge of the target distribution (as reflected in $\\eta_q$ and $\\theta_q$), which is typically unavailable in real-world machine learning scenarios.\n\n- The analysis in Section 4 relies on linearizing the gradient flows around the optimal solution (i.e., the target distribution). Again, it assumes access to the very distribution we aim to approximate. Furthermore, it remains unclear how one can rigorously derive concrete results for the original (nonlinear) dynamics based on this linearized analysis.\n\n- Some mathematical arguments need a careful scrutiny, as detailed in my Requested Changes below.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper revisits gradient flows of the (both exclusive and inclusive) KL divergence over the probabilistic simplex. To be specific, the authors studied the convergence properties of two Euclidean gradient flows as the probabilistic simplex admits two different coordinates (mixture and exponential family) and the natural gradient flow (Riemannian gradient flow) which is coordinate-independent. They showed that, for the objective being the inclusive KL, the Euclidean gradient flow w.r.t. the mixture coordinate exhibits faster convergence than the natural gradient flow, which, in turn, exhibits faster convergence than the Euclidean gradient flow w.r.t. the exponential family coordinate. They further showed that this \"sandwich\" relation also holds after an affine transformation of the coordinates. Analogous results were derived for the exclusive KL. \n\nThis main result shows that the common argument in the literature which is \"Natural gradient descent is faster Euclidean gradient descent due to having straighter trajectories\" is quite hand-waving and is not always the case. Indeed, in this study, nature gradient flow has a straight trajectory towards the optimal solution, yet converges more slowly than the Euclidean gradient flow w.r.t. the mixture coordinate. \n\nAnother contribution of the paper is the demonstration that (discrete) natural gradient descent is more robust than both Euclidean gradient descents under the presence of multiplicative and additive noise. This finding partially explains the empirically observed superior performance of natural gradient descent compared to Euclidean gradient descent in imperfect environments, such as those involving noise and discretization.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- Proof of Theorem 3: equation (14) only implies that $L_q $ is strictly convex instead of strongly convex. Therefore, it seems to me that it is not enough to conclude the set $S_{\\theta}$ is compact (mass can escape to infinity). Please clarify this point and make necessary changes with appropriate citations for the results used. Similarly, please provide more details for the argument $S_T \\subset B_{\\epsilon}(\\eta_q)$, particularly the transition from the smallness of $L_q(\\eta(T))$ to the smallness of $S_T$. \n\n- Lemma 2 provides a global bound for the Hessian of ${L}_q$ in the $\\theta$ coordinate and a local bound in the $\\eta$ coordinate. In the perspective of a reader, I implicitly assume that there is no global bound in the $\\eta$ coordinate in this case. However, a global bound in the $\\eta$ coordinate exists  in the proof of Theorem 3 (eq. (55)). I suggest incorporating this result into Lemma 2 to make it more complete.\n\n- Proof of Theorem 7: Here $\\Delta(k)$ is time-varying so we can not define the universal $\\rho$ at the end of page 20. \n\n- Figure 4 and the sentence \"...and shows why the natural gradient flow falls in between the two extremes\". I do not think the reference plot of  $f(s)=(1/2) s^2$ can represent the local section of the KL divergence for the Riemannian geodesics. I suggest the authors to draw a third subfigure for the local section of the KL along geodesics.\n\n- Page 5: the language could have been made more precise: $\\langle \\cdot \\rangle_g$ is defined as a dot-product in the tangent spaces of $S_n$, not on $S_n$ itself. You can write \"a Riemannian metric $g$ on $S_n$\" instead. By the way, the tangent space of $S_n$ is not explicitly defined anywhere, please identify it more precisely.\n\n- In Theorem 3 and other places, please say that $c_{\\eta}, \\bar{c}_{\\eta}$, etc. are positive.\n\n- The notations $\\eta_{qi}$, etc. in Lemma 6 and other places in the Appendix are very confusing. I suggest using $[\\eta_{q}]_i$, etc. instead.\n\n### Typos: \n\n- Equation (20), $\\nabla \\phi(\\theta_q)$ should be  $\\nabla \\phi(\\eta_q)$.\n\n- Page 18, line 5 in Section D.3: $S$ should be $S_{\\theta}$.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)