Adaptive gradient descent on Riemannian manifolds and its applications to Gaussian variational inference
Authors: Jiyoung Park
Jaewook J. Suh
Bofan Wang
Anirban Bhattacharya
Shiqian Ma
Venue: ICLR 2026 Poster (2026)
Keywords: Adaptive method
Riemannian optimization
Variational Inference
Primary area: optimization
TL;DR: We propose RAdaGD, a novel family of adaptive gradient descent methods on general Riemannian manifolds, and its applications to Gaussian variational inference.

Abstract: We propose RAdaGD, a novel family of adaptive gradient descent methods on general Riemannian manifolds. RAdaGD adapts the step size parameter without line search, and includes instances that achieve a non-ergodic convergence guarantee, $f(x_k) - f(x_\star) \le \mathcal{O}(1/k)$, under local geodesic smoothness and generalized geodesic convexity. A core application of RAdaGD is Gaussian Variational Inference, where our method provides the first convergence guarantee in the absence of $L$-smoothness of the target log-density, under additional technical assumptions. We also investigate the empirical performance of RAdaGD in numerical simulations and demonstrate its competitiveness in comparison to existing algorithms.

Published as a conference paper at ICLR 2026
ADAPTIVE GRADIENT DESCENT ON RIEMANNIAN MANIFOLDS AND ITS APPLICATIONS TO GAUSSIAN VARIATIONAL INFERENCE
Jiyoung Park
Jaewook J. Suh
Bofan Wang
Texas A&M University widydd1 55100t151000a00 awm.tammm..mmmmmmmmmm edu
Rice University jacksuherice.eet
Nanjing University
wangbi LangbE@smail smail..mm.......m.mmm..mm. nju edu, en
Anirban Bhattacharya
Shiqian Ma
Texas A&M University
Rice University sqma@rice.edh
anirbanbestat tamu.edu
ABSTRACT
We propose RAdaGD, a novel family of adaptive gradient descent methods on general Riemannian manifolds. RAdaGD adapts the step size parameter without line search, and includes instances that achieve a non-ergodic convergence guaran- tee, f(*+) fxx) < O11/k), under local geodesic smoothness and generalized geodesic convexity. A core application of RAdaGD is Gaussian Variational Infer- ence, where our method provides the first convergence guarantee the absence of T-smoothness of the target log-density, under additional technical assumptions. We also investigate the empirical performance RAdaGD in numerical simula- tions and demonstrate its competitiveness in comparison to existing algorithms.
INTRODUCTION
Riemannian optimization is an actively studied area due to its various applications relating to ma- chine learning, computer vision, and statistics (Fillard et al., 2005; Bačák, 2014a; Hosseini & Sra, 2015; 2020; Kim al., 2025a). In this paper, we consider the Riemannian optimization problem
min f(x) zEN
where NCM is a geodesically convex subset of a Riemannian manifold M, and f: N + R is continuously differentiable geodesically convex function. A standard approach to solve this problem Riemannian gradient descent (RGD): Fk++ exp2, 8k. Grad f(xx))
Here, exp(-) denotes the exponential map at x, is the step size at iteration k. and Grad represents the Riemannian gradient. Asymptotic properties of RGD have been extensively studied in the liter- ature; see, for example, (Absil et al., 2008; Bačák, 2014a). When f is additionally geodesically L-smooth function, RGD is known to achieve convergence rate f(xx) f(x+) O(1/k) with constant step size
Z (Zhang & Sra, 2016).
However, the of RGD with fixed step size of } has limitations since:
assuming an L-smooth function with geodesic convexity is arguably a strong assumption, requires prior knowledge of the smoothness parameter L.
To overcome this issue, the study of adaptive algorithms in Euclidean space has emerged recently (Malitsky & Mishchenko, 2020; 2024; Latafat et al., 2024a;b; Oikonomidis et al., 2024; Zhou al., 2025; Li & Lan, 2025; Suh & Ma, 2025). However, to the best of our knowledge, the extension of these methods to the Riemannian optimization setting remains largely unexplored at this stage.
Equull contribution, alphahetically ordered.

---

Published as a conference paper at ICLR 2026
In this work, we propose Riemannian Adaptive Gradient Descent (RAdaGD), inspired the by Eu- clidean adaptive method proposed in Suh & Ma (2025). RAdaGD is a family of novel line-search- free adaptive algorithms
Riemannian manifolds, which includes instances that guarantee the convergence rate f(xx) - [(x.) S (11k), similar to RGD with constant step-size. We relax the usual I-smoothness assumption to local smoothness and show that this condition is satisfied by all Ca functions.
Building on this local smoothness setup, we propose new algorithm with provable guarantees. Our main focus is Gaussian variational inference (GVI), which approximates target distribution и on R", most prominently a Bayesian posterior distribution, by non-singular Gaussian. Interest in GVI intimiately connected to the celebrated Bernstein- von Mises theorem, which states that under suitable regularity conditions on the likelihood function and prior distribution, a Bayesian posterior asymptotically approaches Gaussian shape (van der Vaart, 1998, Chapter 10). Yet, existing GVI algorithms provide convergence guarantees only for L-smooth target log-densities (Lambert al., 2022: Quiroz et al., 2023; Diao al., 2023; Kim et al., 2023). To the best of our knowledge, ours is the first method with provable guarantees the non-L-smooth setting. Moreover, numerical results show it is competitive with outperforms prior algorithms.
1.1 CONTRIBUTIONS
We summarize the contributions of this paper as follows.
We provide RAdaGD, a family of novel adaptive gradient descent methods for Riemannian opti- mization. RAdaGD automatically tunes the step size by approximating the local smoothness pa- rameter, without an expensive calculation such as line search. We provide instances that achieve a non-ergodic convergence rate fxxx) - f(x.) S (1/kk) To the best of our knowledge, this is the first of its kind for Riemannian adaptive algorithms.
While RAdaGD inspired by Suh & Ma (2025), its extension to Riemannian manifolds requires highly nontrivial arguments. Our algorithms differ from Suh & Ma (2025) in the step size rule and the analyses, specifically to handle the parameter Ç related to manifold curvature. We first provide the analysis when is bounded and known, and also consider cases when it is unbounded or unknown.
Our results provided for locally geodesically smooth and generalized geodesically convex functions. Although generalized geodesic convexity can considered a somewhat strong as- sumption, local geodesic smoothness is weaker assumption compared to the conventional L- smooth assumption, allowing a broader class of functions. In particular, we show that all C2 functions complete Riemannian manifold are locally geodesically smooth.
Leveraging the advanced properties of the locally smooth setup discussed above, we show that RAdaGD can be applied Gaussian variational inference (GVI) under an additional technical assumption. To the best of our knowledge, our algorithm is the first to guarantee convergence for this application without assuming L-smoothness of target den-itysity. While having a supe- rior theoretical convergence rate, we also show that our algorithm performs among the best in practice. We also validate the effectiveness of our algorithm through numerical experiments in possibly negatively curved spaces.
2 PRELIMINARIES
We provide the organization of prior works in Appendix A. A brief introduction to the Rieman- nian notation is given in this section, while the detailed explanation is deferred to Appendix B.I. Throughout this paper, we denote M as n-ddimensional differentiable manifold and (: as the inner product defined on each tangent space T,M, where M. We also denote || := ((u,, v). and omit the subscript when there is no ambiguity. For x, y E M, we denote the distance by d(x, v). We denote geodesic curves (curves with zero acceleration) by y. A geodesic that satisfies 7() x, ; dt is called a minimizing geodesic. The exponential map exp,: T.M + M is defined by exp,, (v) = 7(1). where (0) = x and 7((0) = v. Here, T.M is the tangent space at We call the locally well-defined inverse the logarithmic map and denote it
(11) = y. and d(æ, y) S/ [+(
2

---

Published as a conference paper at ICLR 2026
by loge. The parallel transport : T.M + TyM transports tangent vectors parallelly along the geodesic curve y such that 7(0) = z and 7(1) = =y
Given differentiable function f M R. the Riemannian gradient of f at æ, denoted by Grad f(x) E T.M is the tangent vector satisfying (Grad fz),, v) = (v) fod all v e T,M. We say N M is geodesically convex subset of M if for all
Here, afz is differential of x. E N there exists a geodesic such that 7(0) = I, 7(1 y. and 7(t) E for all E [0, 1].
We say differentiable function N R is geodesically convex if
f(U) 2 (x) (Grad ((x) log, v),, V
We impose the following assumptions the manifold and the function, which standard in Riemannian optimization literature (Alimisis et al., 2020; 2021; Kim & Yang, 2022).
Assumption 2.1. The manifold M and the target function f satisfy the the following conditions.
(i) M a smooth complete manifold, and N S M a geodesically convex subset of M. N has sectional curvature lower bound Kwin -oc.
(ii) Exponential maps, logarithmic maps, and parallel transports are all well-defined via the minimizing geodesic on N. In particular, the exponential map and parallel transports are computationally tractable.
(iii) The function f ::N - R admits a minimizer on N.
Based on the first assumption, we define a parameter ( [0, 00) - [1, 50) related to K min as
if K Kmin 2 0, otherwise.
((p) (ov lov- coth (pv-Kmin)
(2)
In most cases, we will consider p d(r, x.) for some I E N. Lastly, we follow the arithmetic convention that 1 oo and
0.
LOCAL GEODESIC SMOOTHNESS & GENERALIZED GEODESIC CONVEXITY
In this section, we define local geodesic smoothness and generalized geodesic convexity, which are central our analysis. We begin with local geodesic smoothness. This assumption allows for wider range of functions compared to the conventional geodesic T-smoothness assumption usually considered in gradient descent. We discuss the details in the next subsection.
Definition 3.1 (Local geodesic smoothness). function f M -> R. is called locally geodesically smooth if for every compact set K C M there exists a constant Lg such that
I Grad (xx Grad f [w)ll, S Lxd(a, v), Vx,y E K.
We next define generalized geodesic convexity (Ambrosio al., 2008, Definition 9.2.4), also known as convexity along generalized geodesics. While the definition was initially formulated in Wasser- the extended definition for general Riemannian manifold introduced in Park et al. (2025, Definition C.1). As mentioned in Altschuler et al. (2021, $ A.2), the added flexibility of generalized geodesics is sometimes crucial for the study of RGD.
stein space,
Definition 3.2 (Generalized geodesic convexity). A function f N IR is called generalized geodesically convex with base z E M for all x. y E N, we have,
f(v) 2 (x) + (TS Grad (x), log, y log, m)=-
(3)
f is called generalized geodesically convex if (3) holds for all E N.
We clarify that generalized geodesic convexity implies geodesic convexity (1) by taking 2 = x. and therefore constitutes stronger assumption. Still, note that various applications have been studied under this assumption (Salim et al., 2020; Chewi et al., 2020; Ahn & Chewi, 2021; Altschuler et al., 2021; Diao et al., 2023; Park al., 2025).
3

---

Published as a conference paper at ICLR 2026
3.1 IMPORTANT PROPERTIES OF THE SETUP
We now show that every (M) function on a complete Riemannian manifold is locally geodesi- cally smooth. Even though this statement is somewhat trivial in Euclidean space, we were unable to find its formal extension to Riemannian manifolds. Since this property a core ingredient that enables fruitful applications of our setup, we state it as theorem emphasize its importance. Theorem 3.3. Let M be a smooth complete Riemannian manifold. Assume that f E (M). Then, is locally geodesically smooth.
We provide the proof in Appendix C.I. Theorem 3.3 shows that the local geodesic smoothness is nice relaxation of the usual geodesic smoothness. For example, the squared distance function x + d' (Lx,))) p) is not globally geodesically L-smooth on non-Euclidean Hadamard manifold (Criscitiello & Kim, 2025), but it is locally geodesically smooth it is a C² function.
Now we prove the local version of the co-coercivity condition on a Riemannian manifold. This extends the Euclidean version proposed in Suh & Ma (2025, Lemma 1). The proof is deferred to Appendix C.3.
Proposition 3.4. Let K be compact set on N. Assume that f N + R is locally geodesically smooth and generalized geodesically convex. Suppose K C B(c, R) and K == BR)) S N for some E N and R > Denote the local smooth parameter of K as L Then, for all x, y E K,
1
f(x) f1) (Grad f(y), fog,, 2),, a),
[ry Grad fly) - Grad f(x) 0.
(4)
2L
CONVERGENCE ANALYSIS OF RADAGD
Our Riemannian Adaptive Gradient Descent (RAdaGD) method is described in Algorithm 1. Our algorithm is tractable as long as the exponential maps, Riemannian gradients, and parallel transports computationally tractable.
Algorithm Riemannian Adaptive Gradient Descent (RAdaGD)
1: Input: 2o E R", so > 0, {Axhe20 (Bxta2o [Bx A-1 = 0
2: for k 0, 1.... do 3:
Tk+1 expx, (-8k Grad f(x))
rex 0kt1
Grad (*R+11) Grad (x*x
1
Li+1
26
(5)
2 {(*++1) - ((xx s8 (Tikes Grad ((*++1), Grad (ax))
[A=1++ +1 As
Bz+1
sk
1
r = min
As
BB+1
r
(
Au
Ba11
-1
if sp 2
(6)
Br
Le+
8k+1 min
(BSSk) rk
Lkt1
4: end for
Throughout this section, with the definition (2) in mind, we denote
Ç := sup ( (d(x, zEN
Ck C (d(xx,,,
The algorithm above is defined by the sequences As, BR. and Bg. To start with the simplest case, we assume Ç is upper bounded and known. Such an assumption commonly appears in the Rieman- nian optimization literature (Zhang & Sra, 2016; Zhang al., 2016: Alimisis et al., 2021; Kim & Yang, 2022). We provide a representative choice of Ag. BR, and B, that achieves the non-ergodic convergence rate f(x) - f. I O (1/k), which matches the non-adaptive rate. To the best of our knowledge, this rate is the first of its kind for Riemannian adaptive methods.
4

---

Published as a conference paper at ICLR 2026
Theorem 4.1. Assume that Ç and known. Suppose f N + R is localty geodesically smooth and generalized geodesically convex. Suppose { is generated by RAdaGD with A = a(k=1)+1+$, BR a(k 1) and B& a(k + 1) + (, where (0, 1]. Define Lo as
r Grad f f0o) - Grad f(xo)) Zg
Lo
(8)
||ogm Toll.o
with arbitrary To and initial point 2o such that 2o # To. Letr = ( B
80 - To
1
R² = d( (= Boo+Ca Co 2 |[Grad [x)))] and I. be a smoothness parameter of f on Ban( B3d(Ev.2.) (T.). Then the following convergence rates hold for function value:
L
L
f(xx) ff 2r a:k + +1++++
R² =
(9)
k
Also, the following convergence rates hold for gradient norm square:
2L
L2
min (E(k....))
||rrad (x;))1 S x2 a(k + 1)(k + 2)
182 = 0
1.2
4.1 PROOF OF GENERALIZED CHOICE OF Ak. B₁ AND BN
To prove Theorem 4.1, we establish the convergence rate for a more general choice of An. Bg and By. We start by establishing a Lyapunov analysis. While our analysis is inspired by Suh & Ma (2025, Proposition 12), B3 is introduced as a new parameter to handle Çu. We provide the proof in Appendix D.I.
Proposition 4.2. Suppose f : N + R is locally geodesically smooth and generalized geodesi- cally convex. Suppose {#*]k20 is generated by RAdaGD with {Ak]kzo {Baha20 and { (Ba x20 satisfying
Ag +
As B&+1
Ak
Bx+1
2 1,
21,
2r,
Ak 2 0,
(10)
Ak+1
Av
with some r E (0, Moreover, let 80 = r max
AAo Ao S Ta where Lo is defined in (8). Define
A 0, B-1 Bo + So. x Fo. 8-1 = So, and the Lyapunov function as
V₂: 5k+1Ax - f.) s Bz ||Grad fxx))1
1
|ogg. I,
Vk 2 11,
(11)
Fix k 2 -1. Suppose that the following inequality holds:
B 2 BB41 + Sx+1.
(12)
Then the following inequality holds:
Va+1
Vv (s + 8k+1 $642A0+1((f. f((x+1))
A 9k+1
BB 2
(13)
min
2 sk [(x))]
f(xx)]l S.
Lk
We immediately obtain the following boundedness result from the above Lyapunov analysis.
Corollary 4.3. In addition to (10), suppose that (12) holds for all k 2 Ko - 1, for some ko 2 0.
Then, {2x Jezo is bounded. Specifically, D& BB Bx(x.) for all k 2 Ko + 1, where R 2V Voo--1
Next, we establish lower bound on the step size 8k- We can obtain this lower bound by leveraging the boundedness of the iterate 2 E BR(x.). The proof is deferred Appendix D.2.
Lemma 4.4. Let if be geodesically smooth and generalized geodesically convex. Suppose {sk] generated by RAdaGD with {Ak}a20 {Bx]xzo {Bx] k20 and so satisfying the assumptions in Corollary 4.3. Let L be a smoothness parameter of I on Ban(x.) Bed(Eo. (xx)) where in is the vector in (8) and R is defined as in Corollary 4.3. Then, we have:
r
Sk 2
Vk 2.
(14)
L
5

---

Published as a conference paper at ICLR 2026
Leveraging Proposition 4.2 and Lemma 4.4, obtain convergence rates for generalized choice of Ag. By and By- The proof is provided in Appendix D.3.
Theorem 4.5. Suppose f M - R is locally geadesically smooth and generalized geodesically convex. Define Lg in (8). Suppose {2z) is generated by RAdaGD with {Ax) {Bk [2k)
] x20 and Sy satisfying the assumptions in Corollary 4.3. Define R as in Corollary 4.3 and L as in Lemma 4.4.
(i) The following holds, which indicates a non non-ergodic convergence rate when limg Ag = oo:
R² = (1))
f(xx) - S Ar
A4
As result, applying Proposition 3.4 with y - Tg. and I = Ix we have
L R² O
L
||Grad (+x)]l V
A
A
(ii) Also, (1k)220 achieves minimum selection convergence rate
min (f(zz) - fx) S
R²
0 (
iE (k)
2 2 2i 8i
and
L
L
min
|Gradd (xx)]]
<
12 0
A (ko.....)))
p2 Ei min (Ai, Bi}
Eisko As
(iii) Lastly, we have lim Zk F. for some minimizer T., when lim infę min {AR. B& >> 0.
Readers may notice that the convergence results above depend on ko, which was mentioned in (12). However, we can set Ko = 0 by properly defining Bx. as demonstrated the corollary below. Note that the choices of {Ax [22 {Bu and [B
k20 in Theorem 1 are based on Corollary 4.6.
Corollary 4.6. Assume that Č < oc and known. Define Bg B& + for k 2 Suppose {Ax
and { B&& )xo are chosen to satisfy (10). Then Theorem 4.5 holds true with
ko = O, R² = d(zo, x., 2 + (Bo + Go)s |/rrad f(xa))1²
In the above corollary, we assumed that the value of Č is bounded and known. This is standard assumption in previous literature (Zhang Sra, 2016; Zhang al., 2016; Alimisis et al., 2021; Kim & Yang. 2022) and is unproblematic in the non-negatively curved setting, where { = 1. However, if manifold is negatively curved, such an assumption can be restrictive. In the next subsection, we introduce two corollaries that extend to the case where each assumption is dropped.
4.2 EXTENSIONS TO POSSIBLY NEGATIVELY CURVED SPACES
Assuming Č < oc can be ensured when the whole manifold is bounded, but such an assumption can be thought of as too restrictive. The following corollary demonstrates that this boundedness condition can be weakened through an inductive argument. Specifically, it suffices to assume that the initial distance d(zo, x.. is bounded above by a known constant do. from which we can derive an explicit upper bound To for Go- The boundedness of subsequent ÇK (for & > 0) then follows inductively from the algorithm's update rules, and the choices of (Aa). (B&), ( (Be)- The proof is provided in Appendix D.5.
Corollary 4.7 (Extension for unbounded }). Assume that d(zo, xx)) is bounded above by a known constant dg- Suppose (Auhe2 and (Bx) Ju20 are chosen to satisfy (10). Denote To = (RR) where R² da (Bo {( |Grad f(x) 11? Set Br BB + To all k 2 0. Then
CA S Co=
v 0,
and Theorem 4.5 holds true with koo 0, R2 d(zn, xx 2 (Bo Go)o [[Grad [x00)]

---

Published as a conference paper at ICLR 2026
Corollary 4.7 tells us that the boundedness of in Corollary 4.6 is not mandatory. While Corol- lary 4.7 offers a weaker assumption compared to Corollary 4.6, it suffers from the requirement of prior knowledge of Zn. The following corollary states that, as a trade-off, if Ç is known to be bounded, we still achieve the same asymptotic convergence rate without any prior knowledge of Ç. A similar idea-replacing the known Ç condition with a boundedness assumption also explored in Dodd et al. (2024).
Corollary 4.8 (Extension for unknown (). Assume that Ç < oc. Suppose that {Akkz2 {Bk k20 and {Bx) satisfy (10) and limg-=oo (Bx- - Bx) = oc. Then there exists ko 2 0 that satisfies (12), and the convergence results in Theorem 4.5 hold true.
Remark 4.9. Note that we can find proper Ag. BR. and B& that satisfy Corollary 4.8 without knowing ? explicitly. For example, we may consider alternative choices As = a(k + 1) + 1, Ba ā(k 1), and B& a(k + 1) for Theorem 4.1, with ā E (0, a),, when { is unknown. We provide the details in Appendix D.6.
Although Theorem 4.5 and its corollaries establish convergence rates under the conditions of Propo- sition 4.2, some readers may still wonder about the respective roles and the design mechanism behind Ag. Bg, and Bg. Basically, all the sequences As. Bg. and B3 are carefully designed to sat- isty the incqualities in (23), thereby ensuring the validity of the Lyapunov analysis in Appendix D.I, which forms the core of the convergence analysis. We provide further discussion in Appendix D.7.
One may naturally ask how the theoretical guarantees of RAdaGD depend on the curvature of the underlying manifold. To illustrate this dependence, we provide an analysis in the simplified setting where is geodesically L-smooth, which makes the curvature dependence of our method explicit. Consider the parameters of Theorem 4.1 with a = 1. Under this setup, careful calculations yield
3 3+5
(1+ (2 (1++
((x) - S+
Le "dl.oo....
2(k+2+
(3
On the other hand, for the classical RGD with fixed stepsize (Zhang & Sra, 2016, Theorem 13), one gets
Č - f 1+) d( x2,², for all k 2 1.
When comparing this non-adaptive RGD with our method, convergence rate includes extra constant factor of order ((). This factor is harmless when the manifold has nonnegative curvature (in which case Ç 1, and only the 17/2 constant appears), but may become large on manifolds with strong negative curvature. At this moment, it remains an open question whether this additional factor is an inherent cost of adaptivity or merely an artifact of the analysis.
Lastly, we conclude this section with providing remark regarding the computational complexity of RAdaGD. general, it difficult to provide a universal computational complexity comparison be- tween RAdaGD and standard Riemannian gradient descent (RGD). The additional cost of R AdaGID depends on the underlying manifold, as each iteration requires evaluating two Riemannian metrics and performing one parallel transport in the computation of L Consequently, the overall com- putational complexity is governed by the cost of these geometric operations. In our main application in Section 5, these additional operations have the same time complexity the basic steps of RGD. Therefore, the asymptotic time complexity remains unchanged.
5 APPLICATIONS
In this section, provide core application of our method, Gaussian variational inference (GVD). GVI can be formulated as an optimization problem on a non-negatively curved Riemannian mani- fold. Additional applications under different geometries and tasks can be found in Appendix G.
The parameter choices we focus on in the experiments are as follows, which we found useful after testing various candidates:
Lemma 5.1. Suppose that the manifold has non-negative curvature. Then A =((k: 2) 1 + 2,
B₁ 2(k 2)0.1, and By 2(k + 2)111 + satisfy the conditions in Corollary 4.6.
7

---

Published as a conference paper at ICLR 2026
We note that these choices achieve non-ergodic convergence rate f(x)) f(x)) O (1/19 The proof of Lemma 5.1 can be done by simply checking that (10) holds, as we provide in Ap- pendix E.1.1. While the parameters in Theorem 4.1 may be theoretically sound, they suffer from conservative stepsizes, as discussed in Suh & Ma (2025, $ 5.3).
5.1 GAUSSIAN VARIATIONAL INFERENCE
The goal of GVI is find the Gaussian distribution closest (in Kullback-Leibler (KL) sense) to the target distribution which has density with respect to the Lebesgue measure given by x(x) x e- e-((x) for some measurable function V _R + R. Denoting by (WR)") the space of non-singular Gaussians on IR", this optimization problem can be written as
min F(u)= == Ex-ulx(x)] + Ex (oomallog((())) AE (R
(15)
V=V(a)
:=((a)
for the functional F (W") R. Here, V is called the potential function (or simply the potential) of π. V denotes the potential energy, and H. denotes the negative entropy. It is well known that F(H) equals (KI (1 ||
up to an additive constant. For readers unfamiliar with GVI, one may view F as the objective function f and as the variable x used in Section 4.
favorable aspect of this problem is that BW(R)) admits Riemannian structure, namely the Warssssnsssrttninn (BW) geometry (Bhatia et al., 2019; Chewi et al., 2020; Altschuler et al., 2021; Han et al., 2021; Thanwerdas Pennec, 2022a;b; Thanwerdas, 2022)². Any non-singular Gaussian measure can be parameterized by its mean and covariance, i.e., (m, Σ) E R" x SPD(n), where SPD(n) denotes the space of positive definite matrices in _= Under the BW metric, (W))) becomes product Riemannian manifold R= x SPD((n) with non-negative curvature (Takatsu, 2009, Theorem 1.1). Consequently, in this setting we have < = 1 (Takatsu, 2009, Theorem 1.1). We throughout use i and (m, Σ) interchangeably to denote points (R*). Further details on the Bures Wasserstein geometry are provided in Appendix B.2.
This geometric structure naturally suggests applying RGD methods, in particular RAdaGD, solve (15) with objective function f= F. To this end, we first check whether F is generalized geodesi- cally convex and locally geodesically smooth, properties determined by the target distribution through its potential V. Lemma 5.2 gives a sufficient condition for F to be generalized geodesically convex. Although this result is well known Wasserstein geometry (Ambrosio et al., 2008; Salim et al., 2020; Diao et al., 2023), we include a proof in Appendix E.1.2 for completeness.
Lemma 5.2 (Generalized geodesic convexity of F). F generalized geodesically convex on (WR" V is conver in IR
We next establish the local geodesic smoothness F. To the best of our knowledge, this property has not been identified in either the Riemannian optimization or Wasserstein geometry literature. The result is nontrivial, since the geodesic incompleteness tion of Theorem 3.3 (Thanwerdas & Pennec, 2022b). Nevertheless, the same conclusion holds in BW(R"), as shown in Proposition 5.3. The proof is given in Appendix 1.1.
BW (IE") prevents a direct applica-
Proposition 5.3 (Local geodesic smoothness F). The following facts hold:
(i) Every (BW((())) function is locally geodesically smooth on (WRR"".
(ii) V E C*0 (BW(()) if the following holds for some C, > 0 and B E (0,
[V(X) < C(1 [Xx1') )
(16)
(iii) H E C (BW((())))
(iv) Accordingly, F. locally geodesically smooth on BW(R") if (16) holds.
Lemma 5.2 and Proposition 5.3 together imply that J is locally geodesically smooth and generalized geodesically convex whenever the target potential V is convex, i.e., π is log-concave, and satisfies
'There exists another natural Riemannian metric on the space of Gaussian measures, namely the Fisher-Ran Rao metric. However, for GVI, the Barres-Wasseree geometry is thcoretically preferable. See Appendix G.1.1 for detailed comparison.
8

---

Published as a conference paper at ICLR 2026
the rather mild growth condition. Many statistical models involve such potential V. ranging from classical L-smooth potentials like Gaussian and logistic regression to more complex non non-L -m--ohhooth potentials. In particular, Lemma E.2 highlights examples including generalized linear models such as Poisson regression and Gamma regression with log-link (McCullagh & Nelder, 1989).
With this background, apply the RAdaGD algorithm to solve GVI (15) when V is convex and satisfies (16). The extension of RAdaGD and Corollary 4.6 to GVI problem is nontrivial, due to the geodesic incompleteness of BW (. Nevertheless, with minor algorithmic adjustments and an additional technical assumption, the same convergence guarantee holds, as shown in Corollary 5.4. The modified algorithm WAddGVVI and its proof are provided in Appendix E.1.4.
Corollary 5.4 (Adaptive Bures- Wasssesennn gradient descent). Consider the optimization (15) where V is convex and satisfies (16). Choose {At) {Bx) and (Bx ) to satisfy the assumptions in Corollary 4.6 with Ç 1. Let 8 (0.1) be an arbitrarily small number. Run RAdaGD to solve (15), with an additional minimization step on the step size (which corresponds to BWAdaGVI);
1-8
Sk+1 f- min Sk+1,
(17)
max; Ax xp4+1 V VV)
If there exists 6 > such that Amin ( 2 c foor all O, then Corollary 4.6 holds with N M BW (TR) and r 5 1 5.
Remark 5.5. The additional minimization of the step size in (17) and the uniform eigenvalue lower bound in Corollary 5.4 are introduced to address geodesic incompleteness. Bounds extreme cigenvalues commonly appear in the Wasserstein geometry literature; see discussions on Xu & Li (2025, Assumption 2). Particularly, Lambert al. (2022, Lemma 6) enforced such condition via maximum eigenvalue clipping under strongly convex V, while we did not pursue this direction. We also note theoretically ideal cases where the condition holds (Appendix E.1.4). Empirically, the assumption appears benign, as eigenvalues of the iterates stay bounded away from zero (Figure 1),
The ingredients of RAdaGD for BWR"") space- the exponential map (Definition B.1 BW- gradients (Definition B. 13), and parallel transports (Definition B.15)-are explicitly known, ensuring tractability. Both the exponential map and parallel transport are fully computable. A technical con- sideration is that F and Gradew F involve expectations (Remark B. 14), which may lack closed- form expressions. However, many practical statistical problems admit closed-form BW-gradients (Lemma E.2); when they do not, stochastic approximations can be used (Ranganath et al., 2014; Kim et al., 2023; Lambert et al., 2022; Diao et al., 2023; Luu et al., 2025), although we defer ex- tending our theory to this generality to future work. Thus, the algorithm remains tractable for wide range of applications.
Corollary 5.4 improves upon previous GVI algorithms by replacing the global L-smoothness quirement on V with a weaker growth condition. This relaxation is crucial, as existing methods es- tablish convergence only under the T-smoothness assumption (Domke, 2020; Lambert al., 2022; Diao et al., 2023; Kim et al., 2023), whereas many practical problems lack global smoothness. Thus, our algorithm is, to the best of our knowledge, the first to provide provable convergence guarantees for GVI in this broader setting under mild assumptions.
Bayesian Poisson regression is representative case where global L-smoothness fails. For response Y, E N and predictor X; E R" (i = 1,..., (), a Poisson regression model assumee Y. X; ~ Poisson(exp(X? 0)) independently across i. Assuming an improper prior on 0, the potential V in (15) takes the form
8
V( = Σ -
Because of its exponential dependence on 8, this V is not globally L-smooth. However, by Lemma E.2, it is convex, satisfies the growth condition (16), and both V and its BW-gradient admit analytic forms. Thus, RAdaGD is theoretically guaranteed and fully tractable for this problem.
We conduct numerical experiments on Poisson regression GVI, comparing method with the state-of-tth--at algorithm of Diao et al. (2023) (denoted FBGVI). Since FBGVI requires global L- smoothness of V. it lacks theoretical guarantees for this task. For our algorithm, we chose Au, Br. and Bg as in Lemma 5.1, with i 50 and 25, using independently generated X; N(0, T.))

---

Published as a conference paper at ICLR 2026
Smallest cigenvalue tracking
Poisson regression CVI
value
1/21
FBGVI (s = 1/1.01
10
111
FRGVI 1/4]
RAdaGD OIK!
RAdaGD 01K"
10
Pil-Fuu
10h7
TBCV [:2]
*BoV1 (s 1:1.2)
10-11
0.7
RA4ADD OV!
MMdeGG 0X/
Laplace Apsroximation
0.6
10'
10¹
meration
Iteration
Figure 1: Comparison between FBGVI (Diao al., 2023) and our methods on Poisson regression target with & = 50, 25. Left: F(ua) Fmin. where Fmin is the minimum value among all experiments. Right: Minimum eigenvalues Amin (Ex) over the iteration.
and Y; ~ Poisson(exp(02 X:)) for = I...., l. Figure 1 shows that our method outperforms FBGVI across step sizes 1,1/2, 1/4.
Our algorithm extends beyond the preceding example. Within GVI, the growth condition in Propo- sition 5.3 covers many practical, non-L-smooth potentials (see Lemma E.2), where our method remains applicable and effective. Beyond GVI, there are tasks under other geometries that either fully satisfy our assumptions (Appendix G.3.1) or violate them but still see strong empirical per- formance (Appendix G.2, G.3). Additional applications, including further GVI tasks and broader Riemannian optimization problems, are presented Appendix G.
Codes for experiments be found at https: htt///tttuuuu..t///g com ww http//llldddl dyddl 5510 RAdaGD.
CONCLUSION
In this work, propose RAdaGD, family of line-search-free adaptive gradient descent meth- ods on Riemannian manifolds. Instances of RAdaGD achieve a non-ergodic convergence rate f(xx) ((x) < O(1/k) and extend possibly negatively curved manifolds, which is the first of its kind. We establish convergence guarantees for functions that are locally geodesically smooth and generalized geodesically convex. Local geodesic smoothness broadens applicability beyond the standard L-smooth assumption, yielding the first algorithm with theoretical guarantee for GVI under locally smooth targets (BWAdaGVD.
We believe this work opens several directions for future research. Extending Euclidean adaptive acceleration methods (Li & Lan, 2025; Suh Ma, 2025) or non-adaptive accelerated Riemannian methods (Zhang & Sra, 2018; Ahn & Sra, 2020; Kim & Yang, 2022) to develop Nesterov-type ac- celeration on manifolds would be natural next steps. A proximal variant of our method could further broaden its applicability (Chen et al., 2020: Huang & Wei, 2022; Wang et al., 2022; Diao et al., 2023). Another direction is to clarify the fundamental cost adaptivity. As shown Section 4, our method incurs additional (()) factor in the convergence rate compared to non-adaptive RGD. Whether this factor is intrinsic or merely an artifact of the analysis remains open. Extensions in- corporating retraction maps inexact parallel transport, which are relevant for practice (e.g., Appendix 1.1.1), are also of interest. Finally, developing stochastic variants with controlled error in estimating LA-1 would be meaningful (Bonnabel, 2013; Ranganath al., 2014; Zhang al., 2016; Kim et al., 2023; Diao et al., 2023; Luu et al., 2025).
ACKNOWLEDGMENTSS
Authors thank anonymous reviewers for their insightful feedback during the review process. Re- search of Shiqian Ma was supported in part by NSF grants CCF-2311275 and ECCS-2326591.
10

---

Published as a conference paper at ICLR 2026
REFERENCES
P.-. Absil, R. Mahony, and R. Sepulchre. Optimization AeaimiiattionAlgoitttmm on Matrix Manifolds. Princeton University Press, 2008.
Kwangjun Ahn and Sinho Chewi. Efficient constrained sampling via the mirror-Langevin algorithm. In Neural Information Processing Systems, 2021.
Kwangjun Ahn and Suvrit Sra. From Nesterov's estimate sequence to Riemannian acceleration. International Conference on Machine Learning. 2020.
Foivos Alimisis, Antonio Orvieto, Gary Becigneul, and Aurelien Lucchi. A continuous-time per- spective for modeling acceleration in riemannian optimization. Artificial Intelligence and Statis- tics, 2020.
Foivos Alimisis, Antonio Orvicto, Gary Becigneul, and Aurelien Lucchi. Momentum improves optimization on Riemannian manifolds. Artificial Intelligence and Statistics, 2021.
Jason Altschuler, Sinho Chewi, Patrik Robert Gerber, and Austin J Stromme. Averaging on the manifold: Dimension-free convergence of gradient descent. In Neural Infor- mation Processing Systems, 2021.
Shun-ichi Amari. Information Geometry and Its Applications. volume 194 of Applied Mathematical Sciences. Springer, 2016. ISBN 978-4-431-55977-1 doi: 10. 1107/978-44431-557778
Luigi Ambrosio and Nicola Gigli. Construction of the parallel transport in the Wasserstein space. Methods and Applications of Analysis, 15(1):1-30, 2008
Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré. Gradient Flows in Metric Spaces and in the Space of Probability Measures. Lectures in Mathematics ETH Zürich. Springer Science & Business Media, Basel, 2ed edition, 2008.
Attila Andai. On the geometry of generalized gaussian distributions. J. Multivar Anal., 100:777- 793, 2007.
Manuel Amese and Daniel Lacker. Convergence of coordinate ascent variational inference for log- concave measures via optimal transport. arXiv:2404.08792, 2024.
Miroslav Bačák. Convex Analysis and Optimization in Hadamard Spaces, volume 22 of De Gruyter Series in Nonlinear Analysis and Applications. De Gruyter, 2014a.
Miroslav Bačák, Computing medians and means in Hadamard spaces. SIAM Journal on Optimiza- tion, 24(3):152--1566 2014b.
Jonathan Barzilai and Jonathan M. Borwein. Two-point step size gradient methods, IMA Journal of Numerical Analysis, 8(1):141-14: 1988.
Glaydston C. Bento, Orizon P. Ferreira, and Jefferson G. Melo. Iteration- Complex-- of Gradient, Subgradient and Proximal Point Methods on Riemannian Manifolds. Journal of Optimization Theory and Applications, 173(2):548-56) 2017.
Rajendra Bhatia, Tanvi Jain, and Yongdo Lim. On the wassersin distance between positive definite matrices. Expositiones Mathematicae, 155-), 2019.
Anirban Bhattacharya, Debdeep Pati, and Yun Yang. On the convergence of coordinate ascent variational inference. The Annals of Statisties, 53(3), 2025.
David M. Blei, Alp Kucukelbir, and Jon McAuliffe. Variational inference: A review for statisti- cians. Journal of the American Statistical Association, 112(518):859-877 2017.
Stéphane Bonnabel. Stochastic gradient descent Riemannian manifolds. IEEE Transactions Automatic Control, 58(9):2095-2107, 2013.
Nicolas Boumal. An Introduction to Optimization on Smooth Manifolds. Cambridge University
Press, 2023.
11

---

Published as a conference paper at ICLR 2026
Stephen Boyd and Lieven Vandenberghe. Convex Optimization. Cambridge University Press, 2004.
Yann Brenier. Polar factorization and monotone rearrangement of vector-valued functions. Commu- nications on Pure and Applied Mathemalics, 44, 1991.
Timtor-Emmanue Brunel and Jordan Serres. Concentration of empirical barycenters in metric spaces. International Conference on Algorithmic Learning Theory, 2024.
Yang Cai, Michael Jordan, Tianyi Lin, Argyris Oikonomou, and Emmanouil- Vasileios Vlatakis- Gkaragkounis. Last-iterate convergence of adaptive riemannian gradient descent for equilibrium computation. arXiv:2306. 16117, 2025.
José A. Carrillo, Yifan Chen, Daniel Zhengyu Huang, Jiaoyang Huang, and Dongyi Wei. Fisher-rao gradient flow: Geodesic convexity and functional inequalities. arXiv:2407. 1569,, 2024.
Shixiang Chen, Shiqian Ma, Anthony Man-Cho So, and Tong Zhang. Proximal gradient method for nonsmooth optimization over the stiefel manifold. SIAM Journal on Optimization, 30(1): 210-239 2020.
Yongxin Chen, Tryphon Georgiou, and Allen Tannenbaum. Optimal transport for Gaussian mix- ture models. IEEE access practical innovations, open solutions, 7:6269 6278, 2019.
Sinho Chewi. Log-concave sampling. 2024.
Sinho Chewi, Tyler Maunu, Philippe Rigollet, and Austin Stromme. Gradient descent algorithms for Bures-Wasserstiin barycenters. Conference on Learning Theory, 2020.
Christopher Criscitiello and Nicolas Boumal. Negative curvature obstructs acceleration for strongly geodesically convex optimization, even with exact first-order oracles. Conference Learning Theory, 2022.
Christopher Criscitiello and Nicolas Boumal. An accelerated first-order method for non-convex optimization on manifolds. Foundations of Computational Mathematics, 2023a.
Christopher Criscitiello and Nicolas Boumal. Curvature and complexity: Better lower bounds for geodesically convex optimization. Conference on Learning Theory, 2023b.
Christopher Criscitiello and Jungbin Kim. Horospherically convex optimization on Hadamard man- ifolds part Analysis and algorithms. vrXiv:2505.16970 2025.
Claudia Czado. Lecture 8: Gamma regression. Lecture slides on GLM, TU München, 2004.
Michael Diao, Krishnakumar Balasubramanian, Sinho Chewi, and Adil Salim. Korward-hackwant Gaussian variational inference via JKO in the Bures-Wassersse space. In International Confer- ence on Machine Learning. 2023.
Daniel Dodd, Louis Sharrock, and Christopher Nemeth. Learning -attinn--att-Freee Stochastic Optimiza- tion over Riemannian Manifolds. International Conference on Machine Learning, 2024.
Carles Domingo-Enrien Enrichh and fram-Alexanee Pooladian. An explicit expansion of the kullback- leibler divergence along its fisher-rao gradient flow. arXiv:2302. /2229, 2023.
Justin Domke. Provable smoothness guarantees for black-box variational inference. International
Conference on Machine Learning. 2020.
O. P. Ferreira and P. R. Oliveira, Proximal point algorithm On Riemannian manifolds. Optimization, 51(2):257-770, 2002.
Pierre Fillard, Vincent Arsigny, Xavier Pennec, Paul M. Thompson, and Nicholas Ayache. Extrap- olation of sparse tensor fields: Application to the modeling of brain variability. Lecture Notes in Computer Science, 3565, 2005.
A. A. Goldstein. Cauchy's method of minimization. Numerische Mathematik, 4(1):146-150, 1962.
12

---

Published as a conference paper at ICLR 2026
Linus Hamilton and Ankur Moitra. A no-go theorem for robust acceleration in the hyperbolic plane. Neural Information Processing Systems, 2021.
Andi Han, Bamdev Mishra, Pratik Kumar Jawanpuria, and Junbin Gao. On Riemannian optimiza- tion over positive definite matrices with the warrs-WWssersttinn geometry. Neural Information Processing Systems, 2021.
Minyeon Han and F.C. Park. Dti segmentation and fiber tracking using metrics on multivariate normal distributions. Journal of Mathematical Imaging and Vision, 2014.
Wei Han and Yun Yang. Statistical inference in mean-field variational bayes. arXiv:1911.01525. 2019.
Reshad Hosseini and Suvrit Sra. Matrix manifold optimization for gaussian mixtures. Neural Infor- mation Processing Systems, 2015.
Reshad Hosseini and Suvrit Sra. An alternative to EM for Gaussian mixture models: Batch and stochastic Riemannian optimization. Mathematical Programming, 181:187-223, 2020.
Wen Huang and Ke Wei. Riemannian proximal gradient methods. Mathematical Programming. 194 (2):371-413, 2022.
Yiheng Jiang, Sinho Chewi, and Alexande Pooladian. Algorithms for mean-field variational inference via polyhedral optimization in the Wasserstein space. Foundations of Compulational Mathematics, 2025.
Michael Jordan, Zoubin Ghahramani, Tommi S. Jaakkola, and Lawrence K. Saul, An introduction to variational methods for graphical models. Machine Learning, 37(2):183-233, 1999.
Richard Jordan, David Kinderlchrer, and Felix Otto. The variational formulation of the Fokker- Planck equation. Siam Journal On Mathematical Analysis, 29(1):1117, 1998.
Anya Katsevich and Philippe Rigollet. On the approximation accuracy of Gaussian variational inference. Annals of Statistics, 2024.
Jakwang Kim, Jiyoung Park, and Anirban Bhattacharya. Robust estimation in metric spaces: Achieving exponential concentration with a Frèchet median. Artificial Intelligence and Statis- fics, 2025a.
Jungbin Kim and Insoon Yang. Accelerated gradient methods for geodesically convex optimization: Tractable algorithms and convergence analysis. International Conference on Machine Learning, 2022.
Kyurae Kim, Jisu Oh, Kaiwen Wu, Yian Ma, and Jacob R. Gardner. On the convergence of black- box variational inference. Neural Information Processing Systems, 2023.
Kyurae Kim, Yi-An Ma. Trevor Campbell, and Jacob R. Gardner. Nearly dimension-indepeeee convergence of miel-field black-box variational inference. arXiv:2505.21721, 2025b.
Mare Lambert, Sinho Chewi, Francis Bach, Silvère Bonnabel, and Philippe Rigollet. Variational inference via Wasserstein gradient flows. Neural Information Processing Systems, 2022.
Larry Armijo. Minimization of functions having Lipschitz continuous first partial derivatives. Pacific Journal of Mathematics, 16(1):1-3, 1966.
Puya Latafat, Andreas Themelis, and Panagiotis Patrinos. On the convergence of adaptive first order methods: Proximal gradient and alternating minimization algorithms. Learning for Dynamics & Control Conference, 2024a.
Puya Latafat, Andreas Themelis, Lorenzo Stella, and Panagiotis Patrinos. Adaptive proximal algo- rithms for convex optimization under local Lipschitz continuity of the gradient. Mathematical Programming, 2024b.
Hugo Lavenant and Giacomo Zanella. Convergence rate random scan coordinate ascent varia- tional inference under log-concavity. SIAM Journal Optimization, 34(4):3550-371,, 2024.
13

---

Published as a conference paper at ICLR 2026
Brodie A. J. Lawson, Kevin Burrage, Kerrie Mengersen, and Rodrigo Weber dos Santos. The fisher geometry and geodesics of the multivariate normals, without differential geometry. arXiv:2306.01278.222 2023.
John M. Lee. Introduction to Smooth Manifolds. Springer, edition, 2012.
John M. Lee. Introduction to Riemannian Manifolds. Springer International Publishing, 2 edition, 2018.
Tianjiao Li and Guanghui Lan. A simple uniformly optimal method without line search for convex optimization. Mathematical Programming, 2025.
Yulong Lu, Jianfeng Lu, and James Nolen. Accelerating langevin sampling with birth-death. arXiv:1905.09880 09863, 2019.
Yulong Lu, Dejan Slepěev, and Lihan Wang. Birth-death dynamics for sampling: global conver- gence, approximations and their asymptotics. Nonlinearity, 36, 2023.
Hoang Phuc Hau Luu, Hanlin Yu, Bernardo Williams, Marcelo Hartmann, and Arto Klami. Stochas- tic variance-reduced Gaussian variational inference the Bures- Wassersesnntt manifold. Interna- tional Conference on Learning Representations, 2025.
Yura Malitsky and Konstantin Mishchenko. Adaptive gradient descent without descent. Interna- tional Conference on Machine Learning, 2020.
Yura Malitsky and Konstantin Mishchenko. Adaptive proximal gradient method for convex opti- mization. Neural Information Processing Systems, 2024.
Charles C. Margossian and Lawrence K. Saul. Variational inference in location-scale families: Exact recovery of the mean and correlation matrix. Artificial Intelligence and Statistics, 2025.
David Martínez-Rubio and Sebastian Pokutta. Accelerated Riemannian optmization: Handling constraints with prox to bound geometric penalties. Conference on Learning Theory, 2023.
David Martinez-Rubio, Christophe Roux, and Sebastian Pokutta. Convergence and trade-offs in Rie- mannian gradient descent and Riemannian proximal point. International Conference Machine Learning, 2024.
Peter McCullagh and John A. Nelder. Generalized Linear Models, volume 37 of Chapman & Hall/CRC Monographs on Statistics and Applied Probability. Chapman & Hall CRC Press, London, 2nd edition, 1989.
Klas Modin. Geometry of matrix decompositions seen through optimal transport and information geometry. Journal of Geometric Mechanics, 933):335-390 2017.
Xuan Son Nguyen. A gyrovector space approach for symmetric positive semi-definite matrix learn- ing. Computer Vision ECCV, 2022.
Maximillian Nickel and Douwe Kiela. Poincaré embeddings for learning hierarchical representa- tions. Neural Information Processing Systems, 30, 2017.
Frank Nielsen. Fisher-rao distance and pullback spd cone distances between multivariate normal distributions, 2024.
Konstantinos Oikonomidis, Emanuel Laude, Puya Latafat, Andreas Themelis, and Panagiotis Pa- trinos. Adaptive proximal gradient methods are universal without approximation. International Conference on Machine Learning. 2024.
Felix Ott. The geometry dissipative evolution equations: The porous medium equation. Com- munications Partial Differential Equations, 26, 2001.
John Paisley, David M. Blei, and Michael I. Jordan. Variational Bayesian inference with stochastic search. International Conference Machine Learning, 2012.
14

---

Published as a conference paper at ICLR 2026
Jiyoung Park, Abhishek Roy, Jonathan W. Siegel, and Anirban Bhattacharya. Acceleration via silver step-size on Riemannian manifolds with applications to Wasserstein space. The Thirty- ninth Annual Conference Neural Information Processing Systems, 2025.
Xavier Pennec, Pierre Fillard, and Nicholas Ayache. A Riemannian framework for tensor computing. International Journal of Computer Vision, 66, 2005.
Carsten Peterson and Eric Hartman. Explorations of the mean field theory learning algorithm. Neu- ral Networks, 2(6):475 49444,444444 1989.
Matias Quiroz, David Nott, and Robert Kohn. Gaussian variational approximations for high- dimensional state space models. Bayesian Analysis, 18(3):989-0016, 2023.
Rajesh Ranganath, Sean Gerrish, and David Blei. Black box variational inference. Artificial Intelli- gence and Statistics, 2014.
Adil Salim, Anna Korba, and Giulia Luise. The Wasserstein proximal gradient algorithm. Neural Information Processing Systems, 2020.
Filippo Santambrogio. Optimal Transport for Applied Mathematicians: Calculus of Variations, PDEs, and Modeling, volume 87 of Progress in Nonlinear Differential Equations and Their Ap- plications. Birkhäuser Cham, 2015.
Hiroyuki Sato, Hiroyuki Kasai, and Bamdev Mishra. Riemannian stochastic variance reduced gra- dient algorithm with retraction and vector transport. SIAM Journal on Optimization, 29(2): 1444 1472, 2019.
Xu Shi, Rufeng Xiao, and Rujun Jiang. An adaptive algorithm for bilevel optimization on rieman- nian manifolds. Neural Information Processing Systems, 2025.
Michael Spivak. A Comprehensive Introduction Differential Geometry. Publish Perish, Inc., Berkeley, CA, 3rd edition, 1979.
Karl-Theodor Sturm. Probability measures on metric spaces of nonpositive curvature. Heat Kernels and Analysis on Manifolds, 2(1):217-240, 2000.
Jaewook Suh and Shiqian Ma. An adaptive and parameter-free Nesterov's accelerated gradient method for convex optimization. arXiv:2505.11670 2025.
Asuka Takatsu. On Wasserstein geometry of the space of Gaussian measures, 2009.
Yann Thanwerdas. Riemannian and Stratified Geometries on Covariance and Correlation Matrices. Theses, Université Côte d' Azur, 2022.
Yann Thanwerdas and Xavier Pennec, Bures-wasserstein minimizing geodesics between covariance matrices of different ranks. SIAM Journal on Matrix Analysis and Applications, 44: 44:1447-1476. 2022a.
Yann Thhanwerdas and Xavier Pennec. O(n)-invariant Riemannian metrics on SPD matrices. Linear Algebra and its Applications, 661, 2022b.
Constantin Udrişte. Convex Functions and Optimization Methody Riemannian Manifolds. Springer, 1994.
A. W. van der Vaart. Asymptotic Statistics, volume 3 of Cambridge Series in Statistical and Proba- bilistic Mathematics. Cambridge University Press, Cambridge, 1998.
Cédric Villani. Optimal Transport: Old and New. Springer, 2008
Martin J. Wainwright and Michael Jordan. Graphical models, exponential families, and variational inference. Foundations and Trends@ in Machine Learning, 1(1-2):1-305. January 2008.
Bokun Wang. Shigian Ma, and Lingzhou Xue. Riemannian stochastic proximal gradient methods for nonsmooth optimization over the stiefel manifold. Journal of Machine Learning Research, 23:1-33, 2022.
15

---

Published as a conference paper at ICLR 2026
Chengjing Wang, Defeng Sun, and by Newton-CG primal proximal point algorithm. SIAM Journal on Optimization, 20(6):2994- 3013, 2010.
Toh. Solving log log-detemmnnant optimization problems
Haoshu Xu and Hongzhe Li. Wasserstein f-tests for Frechet regression on Bures-Wasserstein mani- folds. Journal of Machine Learming Research, 26(77):1-123, 2025.
Shotaro Yagishita and Masaru Ito. Simple linesearch-free first-order methods for nonconvex opti- mization. arXiv:2509. 111500994670,.. J4670, 2025.
Ho Yun and Byeong U. Park. Exponential concentration for geometric-median-ff in non- positive curvature spaces. Bernoulli, 29(4), November 2023.
Hongyi Zhang and Suvrit Sra. First-order methods for geodesically convex optimization. Conference on Learning Theory, 2016.
Hongyi Zhang and Suvrit Sra. An estimate sequence for geodesically convex optimization. Confer- ence on Learning Theory, 2018
Hongyi Zhang, Sashank Reddi, and Suvrit Sra. Riemannian SVRG: Fast stochastic optimization Riemannian manifolds. Neural Information Processing Systems, 2016.
Danqing Zhou, Shiqian Ma, and Junfeng Yang. AdaBB: Adaptive Barzilai-Bone Borwein method for convex optimization. Mathematics of Operations Research, 2025.
16

---

Published as a conference paper at ICLR 2026
APPENDIX
A PRIOR WORK AND MOTIVATION
A.1 RIEMANNIAN OPTIMIZATION
Riemannian optimization has become an essential framework for solving constrained problems where the solution space admits smooth manifold structure. Foundational works by Udrişte (1994), Absil et al. (2008), and Bačák (2014a) established the theoretical basis for this field by introducing key geometric concepts for optimization purpose. These geometric foundations enable the gen- eralization of classical optimization methods such as gradient descent, conjugate gradients, and Newton's s method to the Riemannian setting.
Early research in this area primarily established asymptotic convergence guarantees. More re- cently, attention has turned toward non-asymptotic convergence analyses. A particular milestone was achieved by Zhang & Sra (2016), who derived non-asymptotic convergence rate for Rieman- nian gradient descent (RGD), obtaining f(xx)= f(x,) ) S O(1/k) under the assumptions of geodesic L-smoothness and convexity, using constant step size 5 = L.L. Further studies of RGD and its variants under a generalized geodesic convex setup are presented in Chewi et al. (2020); Altschuler et al. (2021).
Subsequent efforts extended these results in several directions. Several authors developed acceler- ated Riemannian optimization methods for geodesically (strongly) convex functions (Zhang & Sra, 2018; Ahn & Sra, 2020; Alimisis et al., 2020; 2021; Kim & Yang, 2022). However, full acceleration in the sense of Nesterov has been shown be impossible on Riemannian manifolds with negative curvature (Hamilton & Moitra, 2021; Criscitiello & Boumal, 2022; 2023b).
Beyond the smooth setting, nonsmooth Riemannian optimization has also been explored. Ferreira & Oliveira (2002) proved asymptotic convergence of Riemannian proximal point methods, with non- asymptotic (L/kk) rates later provided by Bento et al. (2017). Chen al. (2020) first developed proximal gradient methods specifically for optimization over the Stiefel manifold with composite objectives.
There are also some results on stochastic optimization over Riemannian manifold. In particular, asymptotic convergence of Riemannian stochastic gradient descent was first shown by Bonnabel (2013), with non-asymptotic rates for geodesically convex functions later established by Zhang & Sra (2016). Subsequent work introduced variance-reduced variants such as Riemannian SVRG and SAGA (Zhang et al., 2016; Sato et al., 2019).
A persistent challenge across these works the reliance on strong assumptions that limit practical applicability. Most convergence analyses that base on geodesic L-smoothness and convexity require: 1) global L-smoothness on a known set; 2) prior knowledge of the smoothness constant L; and 3) boundedness of iterates, either assumed a priori or enforced through projections. Such restriction is crucial on Riemannian manifolds; as noted in Alimisis et al. (2020, 84), there is some sort of duality between convexity and smoothness with respect to the curvature of the manifold. For instance, on non-positively curved manifolds. even a simple function x + d² p) is not geodesically L- smooth. This simple instance implies that the condition I, global L-smoothness, may become very restrictive condition on certain manifold, conjecturally manifolds with non-positive curvature.
Recent years have witnessed efforts to relax these restrictive assumptions. For instance, Martínez- Rubio et al. (2024) and Martinez-Rubio & Pokutta (2023) removed the iterate boundedness assump- tion for RGD and an accelerated variant, respectively. Nevertheless, these methods still require knowledge of the geodesic smoothness parameter I, over predetermined set, limiting their appli- cability in settings where such global information is unavailable.
Our work represents a significant advance beyond the current state of the art by addressing these persistent limitations simultaneously. Unlike most existing methods, our approach eliminates the need for prior knowledge of the smoothness parameter L while maintaining non-vanishing stepsize and achieving ((1kk) non-ergodic convergence rates. This combination of features has not been attained by any prior stochastic or deterministic method in the literature. Theoretically, our method requires only local geodesic smoothness assumptions without needing to assume iterate bounded- ness. While we employ a stronger generalized geodesic convexity condition, this trade-off enables
17

---

Published as a conference paper at ICLR 2026
substantially weaker smoothness requirements that are satisfied by C2 functions and thus widely ap- plicable to practical problems like Gaussian variational inference where global L-smoothness does not hold. We discuss this topic in more detail in the next subsection.
A.2 ADAPTIVE ALGORITHSS
The textbook step size selection of gradient descent
1k++ Zk S% (xx)
s/ 1 for L-smooth convex functions, achieving convergence rate of f(xx) O(1/k).
There are several limitations of this step size selection, which are:
The assumption that the function class requires global smoothness L for V can be restrictive.
(ii) Even when f is T--mmoohh prior knowledge of I. may be difficult to obtain.
(iii) L is global parameter, so + can be too conservative a choice.
One classical approach that addressing this issue is line search (Goldstein, 1962; Larry Armijo, 1966), which has been widely studied in the literature, but its computational cost can be considered a limitation. For a more detailed discussion, please refer to (Malitsky & Mishchenko, 2024, $ 1).
The development of an adaptive step size method that is
(i) line-search-free,
(ii) parameter-free ffre;; adapts the approximated local smoothness parameter to choose the step size, (iii) guarantees theoretical convergence, and
(iv) for locally smooth convex functions
was first proposed by Malitsky & Mishchenko (2020). To clarify the terminology. note that from now on when we refer to "adaptive" we will mean methods that satisfy these four properties." The method AdGI), introduced in this paper, selects the step size based an approximation of the local smoothness parameter using the last two iterates, defined as La = | -V//ts.)-V/(a--ill Even though it achieves an ergodic convergence rate f((x)= f(x.) I O(1/k), where ig is weighted average of æ æi,..., IA, AdGD was the first adaptive method that satisfies four properties mentioned above for convex function minimization.
||a 16-1l
Initiated by Malitsky Mishchenko (2020), various follow-up works studying adaptive methods have since emerged; here, introduce few. Oikonomidis al. (2024) introduced an adaptive method in the prox-grad setup, named AdaPDM. They leveraged two pieces of local information, and C2: = [V V-f(()))
l = -Vffee-)-
((=v-))-V((er ((2)))1)11111-112 while choosing the step
----- size. They also extended their method to primal-dua setting and introduced AdaPGM. Latafat al. (2024a) generalized AdaPGM to a family of algorithms AdaPGer, Malitsky & Mishchenko (2024) also introduced an adaptive method for the prox-gradient setup, called AdProxGD, by sharpening and extending their previous method, AdGD. Zhou et al. (2025) introduced an adaptive method, named AdaBB, based on the Barzilai -Brbeein method Barzilai & Borwein (1988). They chose the
step size based on the approximation L3 (V ")2²* -==t-==-1) and achieved a larger lower
[V((*) 1
k-2+ Ω (k)), closing an open problem raised by Malitsky & Mishchenko (2024). Furthermore, AdaBB also considers extensions to the prox- prox-gradient setup and locally strongly convex functions. Nevertheless, achieving non-ergodic convergence rate remains a challenge for all of these adaptive gradient descent methods, as noted in Malitsky & Mishchenko (2024).
bound on the step size: Σ
The first adaptive method with non-ergodic (or last-iterate) convergence was introduced by Li & Lan (2025) and is named AC-FGM. They choose the step size based on novel approximation of the smoothness parameter,
2 [Vf((*)) Vf(++=1)11
LA = f(*x) f (x+-1) (Vf(xx), Fk-1 1T)
(18)
This clarification is hased (Malitsky & Mishchenko, 2024, 1).
18

---

Published as a conference paper at ICLR 2026
Surprisingly, they achieved not only non-ergodic rate but also an accelerated rate, f(x)))((x)) _ O (1/k2). Leveraging the same approximation of the smoothness parameter, Suh & Ma (2025) sharpened the constant with their new method, AdaNAG, by providing a Lyapunov analysis. In- spired by the Lyapunov analysis, Suh & Ma (2025) also introduced AdaGD, which was the first adaptive gradient descent method with a non-ergodic convergence rate, f(xx) - f(x,) S (11/k), without momentum. As side note, while momentum may not pose significant difficulties in the Euclidean setup, it can become a bottleneck when extending to Riemannian setup (Criscitiello & Boumal, 2022; 2023b).
A.2.1 ADAPTIVE ALGORITHMS FOR RIEMANNIAN OPTIMIZATION AND LOCAL SMOOTHNESS
Similar to Euclidean optimization, RGD J6+1 = exp,, (-sa Grad f(zx)) with s5 = 7 is a pop- ular method for geodesically convex and geodesically L- smoohh functions. While the limitations discussed at the beginning of this subsection also apply to Riemannian optimization, developing adaptive methods can also beneficial.
However, compared to Euclidean optimization, extending to locally smooth functions can be even more critical, since the conventional geodesic T-smoothness assumptions of RGD exclude a number of important applications in Riemannian optimization. For instance, as noted throughout the paper, even the simple function x + (f, p) fails to be geodesically L-smooth on a Hadamard manifold (Criscitiello & Kim, 2025). Consequently, many applications have been considered without fully established theoretical guarantees, even when certain problem parameters are known.
Despite this demand, studies of adaptive methods that satisfy the four properties discussed carlier are limited, particularly regarding the extension to locally smooth functions. This motivated us to develop Riemannian adaptive gradient descent method inspired by AdaGD (Suh & Ma, 2025). Recalling that AdaGD leverages the approximation (18) to achieve a non-ergodic convergence rate, our method RAdaGD employs similar approximation (5). Since the well-definedness of (5) heav- ily relies on the inequality (4), we assume generalized geodesic convexity (Definition 3.2), which often assumed in the study of RGD, as mentioned (Altschuler et al., 2021, A.2). In the meantime, RAdaGD benefits from the extension to local smoothness, as C2 functions satisfy this property (Theorem 3.3), enabling wider range of interesting applications. In particular, RAdaGD enables applications to the Gaussian variance inference problem with a target log-density that lacks T-smoothness, for which we provide more context the next subsection.
A.3 VARIATIONAL INFERENCE
Variational inference (VI) (Peterson & Hartman, 1989; Jordan et al., 1999; Wainwright & Jordan, 2008; Blei et al., 2017) is a task to approximate a target distribution π to some family of distributions M, called variational family. Formally, one can write the variational inference task as follows:
Find μ. := arg min D(ul(T)
(19)
PEM
for some discrepancy metric D between probability measures. The need of such approximation roots from the sampling, often, π is a distribution whose sample is hard to obtain. Natural case is when n is a posterior distribution from Bayesian inference. In principle, one obtain the exact samples of e from Markov Chain Monte Carlo (MCMC) methods. However, MCMC method is often require heavy computational costs, particularly in high-dimensional setting. VI detours this problem by allowing reduced emphasis on sample accuracy in return for computational gains. Specifically, one chooses M in (19) be casy-to-sample distributions. If one obtains μh, by solving (19), one can use samples from μ., as approximated samples from л. In this regard, VI provides a fast and scalable alternative to MCMC.
'While (Cai et al., 2025; Shi et al., 2025) also consider adaptive methods for Riemannian optimization, their problem setups differ from ours (bilevel and strongly monotone settings, respectively), and they only partially satisfies the discussed properties. We also regard methods with nonincreasing step sizes, such as Dodd et al. (2024), as a different line of work. This is because the step size may not capture local smoothness information, as discussed in (Malitsky & Mishchenko, 2024, $ 1), and thus violates (ii). (Yagishita & Ito, 2025) studies nonconvex adaptive methods and also introduces Riemannian variant, but their stepsize is nonincreasing and does not incorporate curvature-aware scaling.
19

---

Published as a conference paper at ICLR 2026
In Equation (19), there two choices that the user have make: variational family M and discrepancy metric D. The choices of these components determine the nature of the optimization problem (19). For the discrepancy metric D, the most canonical choice is KL-divergence. In this setting, if x has a density of the form πt (I:) о (xp(-V(x)) for some measurable function V R > R. one can show that DxL (1 ||)) = F(u) const where F is in (15). Thus, solving the optimization problem (15) is equivalent to performing variational inference with respect to KL KL-divergence. For variational family M, two notable cases are M = 3W(R", the family of Gaussian measures, and M = P(R) 00m the family of product measures. These correspond to Gaussian variational inference (GVI) and mean-fiele variational inference (MFVI), respectively.
The literature on variational inference (VI) can be broadly divided into two strands: analyzing its statistical approximation properties, and developing algorithms with provable guarantees. In this work, we focus on its algorithmic aspects. For statistical approximation, interested readers may refer to Han & Yang (2019) for MFVI and Katsevich & Rigollet (2024); Margossian & Saul (2025) for GVI.
For the algorithmic aspects, the study of algorithms for MFVI are mainly based on coordinate as- cent method, called CAVI (Blei et al., 2017, $2.4). The algorithmic convergence guarantees for MFVI is an active area of research (Bhattacharya et al., 2025; Amnse & Lacker, 2024; Lavenant & Zanella, 2024; Jiang al., 2025; Kim et al., 2025b). For GVI, one of the most direct and stan- dard approaches to solving (15) is to view it as an optimization problem over the parameter space, treating the parameterization 2020; Quiroz et al., 2023; Kim et al., 2023). In contrast, recent developments emphasize respecting the natural geometric structure of the (WR)) space. In particular, beginning with Lambert et al. (2022), works such as Diao et al. (2023); Luu et al. (2025) have studied solving (15) with respect to the Bures Wasserstein geometry. Our work aligns with this direction.
subset of Rd (Paisley et al., 2012; Ranganath al., 2014; Domke,
B PRELIMINARIES
B. RIEMANNIAN GEOMETRY
In this section, review the basics of Riemannian geometry, and introduce the important functional properties. For the full detail, we refer to the classical textbooks, such as Lee (2012; 2018); Boumal (2023).
A n-dimensional differentiable manifold M is topological space that is locally diffeomorphic to an open set in R. At a point x E M, a tangent vector velocity vector of a smooth curve on M that passes through x. A tangent space T,M is the vector space consisting all such tangent vectors. The space TM- IL T, called the tangent bundle.
Often, this tangent space T..M is conveniently expressed in the form of the vector field, which takes a point in a manifold as an input and returns a tangent space vector that point. map X ( - C (M) is called smooth vector field if it is a derivation, i.e., X satisfies
(fg) = (-) + (g).
Here - E M is the input of the function. A set of smooth vector fields on M is denoted by X(M). As the name derivation indicates, one can think of the vector field directional derivative along the direction of the vector field.
A Riemannian manifold (M, (.))))) a differentiable manifold equipped with Riemannian metric (,, which is an inner product on each tangent space T.M and smooth map on x E M. A Riemannian manifold a manifold equipped with an inner product for each tangent space, called a Riemannian metric. Then, the notion of norm naturally extends as ||||| = V(u. v) For xr y E M, the distance d(x y) is the infimum of the length of all piecewise continuously differentiable curves from x to y.
For (Ir, v) E TM, a smooth curve 7. [0. + M with (0) = x and (0) v is called geodesic if it is critical point of the energy functional defined by
/ In (0) dt.
E(7)
20

---

Published as a conference paper at ICLR 2026
For instance, straight line (t) = - to in Euclidean space is geodesic. If geodesic such that (0) and 7(1 = y satisfies d(x, v) = So [t) dt, then we call this geodesic the minimizing geodesic.
Geodesics induce the notion of exponential map and logarithmic map. The exponential map expx M + M is the map defined by (x) (1)... The exponential map transports the point x in the direction of the tangent vector v following the geodesic Ye. is well-known that the exponential map is local diffeomorphism in some neighborhood of 0 E T,M, denoted by U. Thus, the exponential map allows the inverse map on U, which is called the logarithmic map log, cxp,, (U) M. Note that the exponential map and logarithmic map are locally well- defined, but may not be globally. The manifold is called complete if exp,, is defined on entire M.M.
Geodesics also provide a way to transport vectors from one tangent space to another, called parallel transport. A parallel transport (()) Ty(to) M Tr(t,) M is way transport tangent vector along the curve parallely. Given geodesic, the parallel transport need to satisfy the following property:
(s)) >
It is known that the linear operator P with the above property is unique. Hence, the notion of geodesic uniquely determines the parallel transport. 7 is a geodesic curve such that r(0) I, (1) y. then simply denote I())) rÿ, (geodesic) parallel transport from T. to MM
Geodesics and the corresponding parallel transports define an operator called Levi-Civita connection or covariant derivative, which moves tangent vector another tangent vector along the certain direction. Specificaly, a covariant derivative is an operator V. :((M x x(M) - x(M) defined by
V V=(()(((())) = lim a ( X h-0
h
Under Levi-Civita connection, geodesics have zero-acceleration, Meaning, the following equation holds if and only if is a geodesic:
(4)() = 0.
In this regard, often geodesics are denoted by constant speed geodesic.
The last geometric concept induced from the covariant derivative curvature. The Riemannian curvature tensor R(.,)) .). ((M) x ((M) x ((M) - X(M) is defined by the following formula: R(X, Y)Z V v V V
where [. denotes Lie bracket.
The key geometric quantity in our analysis is sectional curvature, which generalizes Gaussian cur- vature 2-dimensional surface, For p E M, denote Ep a set of two-dimensional subspaces in Tp M. For og Mp, the sectional curvature K Σp + IR is defined by the following formula: K K(a,)) (Ruu, u)
where {u, u} is an orthonormal basis of dp-
Lastly, we introduce the notion of Riemannian gradient and geodesic convexity. Given differen- tiable function f .M > R. the Riemannian gradient of f at I, denoted by Grad f(x) e T.M is the tangent vector satisfying the following property:
(Grad f(x), v), = =ff(v) E.M.M
Here, dfz a differential of at x. We say N M is geodesically convex subset of M if for all x, E N there exists a geodesic such that 7(0) z, y(1) y, and ytt) E N for all t E [0, 1]. We say differentiable function N R is geodesically convex if for all y,y N.
fw)) > f(x) + (Grad (x), w)..
We note that the choice of geodesic, parallel transport, and covariant derivative are mutually equiv- alent: once one of these structures is specified, the other two are uniquely determined. In this work, we adopt the perspective of starting with geodesics, defined as critical points of the energy functional, which in turn uniquely determine the associated parallel transport and the Levi-Civite connection. This approach follows Spivak (1979), whereas other texts, such as Lee (2018), take the alternative route of beginning with connection and subsequently developing the other structures.
21

---

Published as a conference paper at ICLR 2026
B.2 BURES- WASSERSTEIN GEOMETRY
In this appendix, we briefly introduce waesssaassrrstninn space (W(R") space of Gaussian mea- sures equipped with W2 metric. This space is known to be a Riemannian manifold (Modin, 2017; Bhatia et al., 2019; Han et al., 2021; Thanwerdas, 2022).
One of the most intriguing aspects of this space is that it admits two equivalent geometric formu- lations: the optimal transport & Otto calculus viewpoint and the Riemannian geometric viewpoint. While both lead to the same underlying geometry, they offer distinct perspectives. Depending on the problem at hand, one may choose the Otto calculus formulation or the Riemannian formulation, as each has its own advantages in different contexts.
In terms of the analyses of gradient methods (including RAdaGD precisely BWAdaGVI), it is more common to take Otto calculus based approaches (Jordan al., 1998; Salim al., 2020; Lam- bert et al., 2022; Diao et al., 2023; Park et al., 2025). In this regard, we mainly provide the Otto calculus formulation of BW (W"),, and briefly mention the equivalence between Riemannian for- mulation. Main takeaways of this appendix is to provide the detail how can implement ingredients in RAdaGD (exponential map. parallel transport, (W-gradients) under (W(R") geom- etry.
We list the results mainly without providing the full derivation. For detail, we refer to Ambrosio et al. (2008); Villani (2008); Santambrogio (2015); Chewi (2024) for general Otto calculus, and Takatsu (2009); Altschuler et al. (2021); Lambert et al. (2022); Diao et al. (2023); Park et al. (2025) for Otto calculus BW("")
B.2.1 ОТTО CALCULUS
We first begin with the general Otto calculus on the space of probability measures, and then focus on (R").. which allows the finer analyses. We denote (R to be the space of probability measures on R") which have the finite pth moment, (IR)
be the space of probability mea- sures on R", which are absolutely continuous with respect to Lebesgue measure and have the finite second moment, and (²(x)= {1 I2 +R" I [I][(=))] dp(x) oo). For u. E Py (BB", let ((,, v) be set of couplings of μ and Wasserstein distance between μ and v is defined as follows. Definition B.1 (Wasserstein metric). Let μ.v E (RR)). Denote ((u, v) to be a set of coupling measures of and v. p- Wasserstein distance between y and is defined as follows:
(H( = inf B(aa))y |||e = JEr(e.)
This is known to be well-defined metric, and the metric space (P,(""), is called p- Wasserstein space.
Brenier (1991); Jordan et al. (1998); Otto (2001) showed (P2.ac (R.) W2) endows a richer geomet- ric properties than mere Pp (TR"). Specifically, while (P2.ae (RR"), ), W2) is not precisely Riemannian manifold, its geometry is almost same to the non-negatively curved Riemannian manifold. The following theorem serves as the key ingredient for such observation (Brenier, 1991):
Theorem B.2 (Brenier Theorem). j,, E P2 )TR",, then
min W (p, v) TEC*(H) st. Tau=u
Evnn [IT(()- - =0|²]
T id id] (=()))
min
TECP Tap
Denote the minimizer Tp... Then Tu is unique p-a.e. and is gradient of some convex function o on R" p-a.e. Furthermore, Ty, Tv. = id. The minima Ty, is called the optimal transport map from μ to v.
Theorem B.2 gives notion of tangent direction at p (Chewi, 2024, $1.3).
Definition B.3 (Riemannian metric in 2- 2-Wasserstein space). For μ (P"),, tangent space
T.P2.0 (IR.) C L² (2) is defined by
T.Po. 2.00 (R") V (R") (P(M) = [X[ - a 0,T is an optimal transport map)" (M)
Here, O (R") is set of compactly supported smooth functions on R". The Riemannian metric here is L (u))--inn product, i.e., (v., w) (xzzu((((x), ((()))
22

---

Published as a conference paper at ICLR 2026
One naturally define geodesic curve in (P2,ac W2), by pushforwarding the interpolation between particles to the measure space.
Definition B.4 (Geodesic in Wasserstein space), geodesic curve T [0. 1] - P2,ac (R*) such that 7(0) and ((1) = y can be defined as follows:
7(t) ((1 t)id + pT,,,) #a
The exponential map and logarithmic map are then defined accordingly.
Definition B.5 (Exponential map and Logarithmic map in Wasserstein space). For jr E (R") and u E (p), exponential map and logarithmic map of (P2.ac (R"), W2) are defined as follows:
(xp = (v + id) #y:
1og = TT.. - id.
Ambrosio Gigli (2008) defined geodesic parallel transport on 2- 2-Wasserstein space. Definition B.6 (Parallel transport). For μ., E P2 (R") and v E T (Rc
F = 1I (U Tw.)). Here, IT. a projection operator (2 + (.P2,ac(R"))
While this definition of parallel transport theoretically sound, in terms of tractability this notion of parallel transport is hard to deal with; it involves the projection operator 1I which lacks an explicit form, In this regard, one often considers dropping out IT and using the un-projected parallel transport Te = v Te as the parallel transport from (v() to L2 (μ)). Such construction is widely adopted Wasssrssn geometry literature; see Ambrosio et (2008); Salim et (2020); Park et al. (2025) for instances. Perhaps surprisingly, this un-projected parallel transport still imposes favorable properties parallel transports.
Proposition B.7 (Un-projected parallel transport). For u E (R,)((R") and v (2),,,, define T : oo Tv, Then,
1. T is linear operator (4)..
2. I . D id.
3. (v. (, = (T T
Proof. The property is direct: for v, w E (1)) and a, E R. [a + bw) = av Tu. + bw Ta, al v bl"
The property 2 from Theorem B.2.
The property 3 is direct consequence of the change of the measure formula:
(u, w),,),, = (u(x), w(x)) (()) = f (v Tzp(z), (m)( du(x) (TTe, Ibw),
By Proposition B.7, we can consider using the un-projected parallel transport Ty. as a parallel transport.
Next, we introduce the notion of Wasserstein gradient. Wasserstein gradient detined analogously to the Riemannian gradient formula df..()) = (Grad f(x),
Definition B.8 (Wasserstein gradient). For a functional F (P2 + [o)), the Wasserstein gradient of F at μ the mapping Gradw, F(u) (²(2) satisfying the following equation:
F.F()))) (Gradw, F(H), vo)
Here, pe is arbitrary sufficiently regular curve of measures on P2 (R") with po μ and ve L2 (pt) is the vector field satisfying the continuity equation d.p -dvv(wee).
23

---

Published as a conference paper at ICLR 2026
One has the following formula::
Gradw, F(1) 6F((1)) 81
Here, is Euclidean gradient and 8F(m) the first variation. 8G
Here, the role of 70) is changed to vo- More generally, one can think of 12 as a velocity at pr playing the similar role to (t) Riemannian manifolds. For the derivation, we refer to Chewi (2024, Theorem 1.4.1).
Consequently, the notions of (generalized) geodesic convexity and local geodesic smoothness can be formulated in terms of this un-projected parallel transport and Wasserstein gradient.
Definition B.9 (Generalized geodesic convexity and local geodesic smoothness in Wasserstein space), Let F (2,,ac(R)) R U {oo} be differentiable functional.
1. F is called generalized geodesically convex with base (2." if for all p.v (Ik))(()))
F(v) 2 F(n) + (Gradw, F(p) TT, T
If for given a and F is generalized geodesically conver with base μ. it is called geodesically convex, and if F is generalized geodesically convex with base a for all a e P2 (R"), then it is called generalized geodesically convex.
2. F is locally geodesically smooth if for all compact set K C (22,a0(8") there exists constant Lx 2 such that
||Gradw, F(v) o Ta. Gradw, F(u)l,, S LKW2(4,
Vu,eeK.
B.2.2 ОттО CALCULUS ON (R))))
Now, we turn attention to specifically to (W(R"), where much finer analyses doable.
The first advantage to consider BW space is that the optimal transport map has the explicit form, particularly the affine map.
Definition B.10 (Optimal transport map between Gaussian). The optimal transport map between
140 N(mo, En) and μ11 = (m1,, E1) is defined as follows:
T Tun.us (x) m2 E 1/2
1/2 m, x
//2
1/2
M 1//2 -m))-
Definition 1.10 provides two favorable results.
First, since afline transform of the Gaussian is also Gaussian, from Definition B.4 every geodesic interpolation between two Gaussians is also Gaussian. This shows (W(R") is a geodesically convex subset of Wasserstein space. In addition it implies (W(R") is totally geodesic submanifold of 2- Wasserstein space (Lee, 2018, Exercise 8.4).
Second, we can identify μ N(m, Σ) = (m, E) E E R" SPD(n) and T,BW( (a, Rn Sym(n). Here, SPD(nn is the space of the symmetric positive definite Raxn matrices, i.e., positive semi-definite matrices, and Sym(n) ) is the space of Raad symmetric matrices. By writing + at u (m,, )) can be parameterized by (a, S). Under this identification, we can view BW 3WR) space product Riemannian manifold of Rd SPD(nn. Then one can parameterize every quantity in Appendix B.2.1 by this product manifold sense. For instance, the vector corresponding to the ",,/",/,///// /2//////////2111111/
an affine map as (xx) S(r- for fixed m (which is the mean of p), any affine map starting
optimal transport map Tho is (m, (2:) 112
We can define Riemannian metric, exponential map, logarithmic map, and Bures-Waas Wasserstein gra- dient in terms of parameters as well. These are basically just parameter-wise representations of the formula in Appendix 1.2.1. We want to emphasize that the formulation we will introduce may be different with some other Bures-Was Wasserstein geometry literature (Thanwerdas & Pennec, 2022b). This is again because we are taking into account the Otto calculus formulation, and these two for- mulations equivalent.
24

---

Published as a conference paper at ICLR 2026
Definition B.11 (Riemannian metric of Bures- Saies-Wassersteiin space), Let (mm, E). The Rieman- nian metric of ( is define by
((ao, So), (ai Si)), (ao, a1)zz + tr(SoES)).
The 2- 2-Wasserstein distance between Gaussian distributions, μ4i N(mis 2ie) for 0, has the explicit form
W (po, I| -m1l tr(Eo + EΙ 2 E2
(Eon
Definition B.12. (Lambert et al., 2022. Appendix B.3) Let bui = Nm,, E.). The exponential map and a logarithm map (R") are defined by
exPHo (a, S)) N (a + mo, + 1) +)),
log (M1) = (m1 - mg. 20 (20 /- Σ1200 /2/22///1//22000
for a - Sym(n).
Definition B.13. (Lambert et al., 2022, Appendix B.3) Bures- Wasserstein metric of the functional F can be written a function IR" x PSD(n), the space of the mean and covariance. Then, for m E and E PSD(n).
GradBW F(m, Σ) Fmm, E), 2V 2 E))
See Lambert et al. (2022); Diao et al. (2023) for further discussion.
Remark B.14 (B (BW-gradients for KL-divergence). In particular for the functional = V H. we considering in our main paper (Section 5.1), gets the following form of the BW-gradients (Lambert et al., 2022):
Gradew V(m, )) (Ex=N(m,E) XX B (VV())
Gradew H(m, 2) 1)
Gradgw F(m, Σ) (Ex~N(m,)) VV(X)) VVV
V is not C or C2, VV and V can be regarded as the weak derivatives.
Using the isometry between the function representation and the vector-matrix representation of T,BW(R" (I2"), can define the following operation, which can be used construct the (un- projected) parallel transport.
Definition B.15. For S) E TH BR) and (b, R) E Tv. BWR"",, we have the following operation.
(a, o (b. R) = (a + Sb - SR).
In particular,
(m) (mi Enj (a, 9) SE 1/2 (E/ ΣイΣリ / /3/
Note the above formulations are just parameterization of the general Otto calculus in Ap- pendix B.2.1, and one can write in either way to proceed. Throughout the paper, we will inter- changably use the transport map formulation and vector-matrix formulations of BW(R") geometry. Precisely, we will mainly use the transport map based formulation for the proofs, and the vector- matrix formulation for the algorithmic aspects.
known that BW(R") space geodesically incomplete, if one restricts to non-singular Gaus- sian distributions (Thanwerdas & Pennec, 2022a). On the other hand, if one includes the singluar distributions, then the space become geodesically complete (Takatsu, 2009), while the space is no longer precisely manifold, but cone. In addition, it is also known that BW(R" is neggavvvyy curved manifold without the curvature upper bound. See Park al. (2025, Lemma A.33, A.42) for the detail.
25

---

Published as a conference paper at ICLR 2026
C DEFERRED PROOFS FOR SECTION 3
C.I PROOF THEOREM 3.3
For any compact set we have diam(K) == D co. We denote B(K, (p d(p, K) s D).
Claim: The below quantity Lx is well-defined and finite:
LK = sup PEB(K,D)
|Hess f(p)llop
To prove the claim, note (KK, D) is compact from Hopf-Rinow Theorem. In addition, from f E C2 (M), the map p + |[Hess f(p)llop continuous. As a result, the function p + |[Hess f(p) op
admits the minimum and maximum compact set B(K, D). This proves the claim.
Take v, E K, let y be the minimizing geodesic with 7(0) and 7(1) = 2. Define
G(() 1 7(t) Grad f(((()).
Observe
T ((+h) Grad f(n(t + h)) rV 7(2)
Grad f(7(t))
(t) = lim
h-+0
h
( r
7(0) 7((+N) Grad f(r(t + h)) TV Grad f(r(t))
7(0)
lim
6-+0
h
pi(t)
7((+)) Grad [()(1 h)) Grad s(>(())
r
lim
200 h
h
n (() 900
V
Grad f()()).
Here, (i) is from the property of the parallel transport and the fact that 7(t), (t+ /)) and y lie on the same geodesic y, (ii) is from the fact that the parallel transport is the continuous linear operator, and (iii) comes from the definition the covariant derivative and the fact that f E (M).
Next, show
V 400 Grad f()())) < |[Hess (()( [[(()II.
(20)
Note this inequality trivially holds if V 4(t) Grad f()(t)) 0. Otherwise, observe the following calculations:
V Grad f(r(t))) <V (()) Grad f()(t)), Va(0) Grad f(((t)))
[Hess (t) V V Grad (((()))
SVV
Hess (()(()) Vace) Grad f (>
W Hess f(r(())ll [|((t)] 5 Grad f(((t))|
Here, (i) is from the definition of Riemannian Hessian, and (ii) is from Cauchy-Schward inequality.
Furthermore, we have r(00, 1]) B(K,D), for all [0,
d()(), s d(v((), < id(x, v) S D.
Together with the construction of Lx. we get
Hess [(T(I)) S LK. eE
sup
(21)
26

---

Published as a conference paper at ICLR 2026
With these ingredients, apply Fundamental Theorem of Calculus to get
1./ G" (t())tt
= $f [[
[T Grad f(z) Grad f())1l
V/ s00) Grad ((7())]
|[Hss f(r(t)) /. [1(())] de
||Hss f(x(tt)llo)) (( (I)( dt - sup te[0,1]
Lgd(x,y). v)-
Here, (i) is from the isometry of the parallel transport, (ii) is from (20), and (iii) is from (21) and the fact that y is the minimizing geodesic.
C.2 AN AUXILIARY LEMMA USED IN THE PROOF OF PROPOSITION 3.4
In this subsection, we state and prove a lemma that is used in the proof of Proposition 3.4. This can be considered local version of inequality (11.11) in Boumal (2023).
Lemma C.1. Let f M R is locally geodesically smooth function. For compact set K C M, there exists LK such that
f(v) S fxx) (Grad fx)), vi= + LK (z 2
Vx, E K.
The proof follows (Park et al., 2025, Lemma A.18). By the definition of the Riemannian logarithmic map, get 70(0) log, y. By Fundamental Theorem of Calculus and properties of the parallel transport,
f(v) = f(((1)) = f(+(0)))))))))))))))))))) + S.'' d (f + 1'' (Grad ((y(t)), dt f()(t)), (0)) di = f(x) f(((rr Grad f(r(t)), v) dt.
= f((x) + (((1110) Grad f(+(t),,
Then, by subtracting f(x) + (Grad fxx), log, u) from the both hand sides,
/''
(()) - (xx - (Grad ((x), log, y) =
(rr 700 Grad ((x()) - Grad f(x), ((xx,,lgz, u)
| Grad ((((t)) - Grad f f(x) |loga ulll de
Lxc (xx y) f idl
(ied()((), y)dl (i
L dl (2.z))
2
For (i) we used Cauchy-Schwarde inequality. For (ii), we used the local geodesic L-smoothness with LK be the local geodesic smoothness parameter for the compact set B(K, D) where D diam (K). This holds as 7(00, 1J) C B(K, D), which was shown in the proof of Theorem 3.3. For (iii) we used the fact that the geodesic curve satisfies d(z, (f)) = td(x, y) due to the constant speed property.
C.3 PROOF OF PROPOSITION 3.4
Proof of Proposition 3.4. We combine the proof strategy of Park et al. (2025, Lemma B.2) and Suh & Ma (2025, Lemma 1). Let == exp, (- (Grad f(y) - ry Grad [(x))). Note x, y E
27

---

Published as a conference paper at ICLR 2026
B(c, R. C 3R), We first show E B(c,3r).
1
d(z,c) S +(c, s ++ + d(y,
-Grad (w)) r Grad / (x))
1
- R +
|Grad f(y) - 1 Grad f(x) S R RR+((z,) v)
LK
2 R + d(I, c) + d(c, v) 3R.
Therefore, E Bee, 3R). Note K is compact by Hopf-Rinow theorem. Now, since E N by the assumption, observe
[(x) - f()) f(x) - f(z) + f(z) f(y)
@
< TT Grad (xx, logy = -OBg 2),, + (Grad [(y), z),
L K
|logy 20 =1(
1
n Grad fx), log, x)
|Grad f(y) - 1 Grad f(x)11?
2
Grad (z), log y)
Grad f(z) - I Grad sw))l
2
We used generalized geodesic convexity and Lemma C.I for (i), the construction of for (ii), and Lemma F.2 for (iii). Rearranging the final formula and swapping x and yield the desired result. C
We additionally provide a useful corollary of Proposition 3.4 in the convergence analysis proofs presented later.
Corollary C.2. (Suh & Ma, 2025, Corollary 2) Let f N - be locally geodesically smooth and generalized geodesically convex function. Assume æ, y E B(c, R) and K == B(c, 3R) S N for some N and R > Then, the following statements hold:
1. If f(x) - f(y) + (Grad (3) logg x)) 0, then
[| Grad f(y) Grad f(x) O.
2. ff)()) - f(v) + (Grad [v). log, xt) + 0, then
r Grad f(v) Grad f(x) f(y) - f(x) + (Grad f(y), log x),,
0
where Tog is the local smoothness parameter of K. Proof. For the first argument, the result directly follows from Proposition 3.4. For the second argument, geodesic convexity and non-zeroness yield
f(x) - f(v) + (Grad f log,, x> > 0.
Since ||.11 > 0, this proves the non-negativeness. For the upper bound, since f is generalized geodesically convex and locally geodesically smooth, Proposition holds. Dividing the both hands sides of (4) by Tw)- [2)),og, leads to the desired claim.
D DEFERRED PROOFS FOR SECTION 4
D. PROOF OF PROPOSITION 4.2
We show the following lemma to show Proposition 4.2.
28

---

Published as a conference paper at ICLR 2026
Lemma D.1. The following equality holds for any k 2 -1:
Vk+1 V&
(|18SRn Nz 1| - 10B Fz -|| S6+2 |0gmm matzl| + 2 + tk+2. log m²))
2
+ (9k+1 Ax + 8k+1
+ sk+1 (f(2++1) - f.+ f4 (Grad f(1++1), (or Z.)
+ 8k++ Av [((R+1)) ((xx) + Grad (fx+1) 10, log Ib
mk11
Do 2641
Grad ((*++1) Grad ((2)))
2L
+ +Q=:
(22)
where if
AR Qk
Bk+1 Sk+1 S441L8+1 ||Grad f(xx+1))l
2
As
Grad ((xx+1), Grad ((xx) + (1 (1+
BR oh Lo+t l| |Grad f(xx)]]
-2 8x bxLk++)
pese 2841
Av. 86+1
and Q-1 O.
Proof. Reorganizing difference of
1og F.
term:
1
log
1ogzzt2
lognets F.
1
[ |Iogrees
|10 5x
2
S44182+1 ||Grad (**++)) 25k+1 (Grad (IT log 22))
2
1
+
S4+1
Grad f(xx+11)
s2+1
Grad f(Ix+1); log Z.
2
1
Nz -| -1
Se+1 logrens 0kt2 122 + 2 (10800e 22, loBeens v..) ms.)
2
1
+
S4+1sh+1 Grad f( + sk+1 2
Grad f(Ix+1); loga.. X.)
Reorganizing difference of 8k+1 1+ (x(( - f4) term:
3642A441 (f - fa) 84 - AA - f.)
= (88441 As 8k+1 8k+2 ((*411) + 8k+1 8k+1 (f(x+))) fa)
8841 Ak [(0241) - ((xx) (Grad I (xx+11 1), logoezs 5hh
1x11
1
Dre Grad f(x41)) Grad ff ((ww))
2L
85+1 Ag
As (Grad f(x++1) logmen
Grad f(xx+1) Grad f(xx)
2
1e+1
29

---

Published as a conference paper at ICLR 2026
Combining two, we have
log logm,_, 2211 || + (f f.)) --21 111 2, a1l" + su+ 14 f())
fi)
1 |108800 2.1l - |108er. Zx 11 Sa+1 logeens 16t2/ | 2 + Ik+2, Ikt2, logxe+1 ItA) 2
+ (s + Sk+t1 Sk42As411 [f+ f(z2+1))
+ Ak f((x+1)) - f(xx) + Grad ((x+11), logax (00441)),gg1111 Fk) 12+1
1
Iete 2641
Grad f(xx+)) Grad ((xw))
2Lx+1
+ Qu.
where
1
Qx Sk+1Sk+1
||Grad f(x++1)l S Grad f(0++1), #k
8k+1 As
T #k11
Grad ((x2+1) Grad of (xxx)
2Lx+1
First, consider the houndary case k = -1. By the definition that A_ 0, B_1 Bo + Co. and T-1 = Jo, we obtain
1 =" B- B- |/raa f(1-1)1? =0
i- +
2 Bo |[rrad (=x))) -
which proves the result.
For 0, we conclude the desired equation by verifying
Qx sk1 Ba+1 [[Grad f(*x+1)12 2 $ Br ||rrad f(x))1?
8k+1 Ak
[1- B++1 + Ca+1
Lx+1 [|rrad fxxx+1)11
8k41
2L
A
(ras FK- Grad [(x++1), Grad ((xa) + (1 BR 22 Lk+1 |Grad ( )!(
2 (1- - 8h Lx411
Ag 8k-1
where right hand side is the definition of Qk. This can be verified by using Lemma F.2 and compar- ing the coefficients carefully.
Leveraging Lemma D.1, we can show Proposition 4.2.
Proof of Proposition 4.2. First, for the fixed (as stated in the Proposition, which satisfies l: 2 -1), we can verify the following inequality by the stepsize rule (7):
A +1
Sk+2 S
A++1 Sk+1-
(23a)
Furthermore, if this A satisfies & 2 0, i.e. k + -1, then the following inequalities also hold based on stepsize rule (7) and the condition (12):
A
8k41 I
(23b)
B&+1 + Cktz
8k+1 5
if sh S 1,
(23c)
B +
GKt1 Liwi
Ak
Bett + Sh+1
8k+1
if if 2 1.
(23d)
Bg
As
L
30

---

Published as a conference paper at ICLR 2026
The desired inequality (13) then follows from applying Lemma D.1 with Ck+1- We now show how each term is upper bounded in Lemma D.1 with Ç = (k+1-
From the definition of Sk+1, obtain ((d(#x+1,x.)) Sk+1- Applying Lemma F.1 with Ik+1. = Tk+2, P = 2+ yields:
[oB 2. =1I2s | |oBBets ax
+ C++1 [oo mk+2 ||
(1080e4s 2k+2 lob Fx)
which implies the first term in Lemma D.I is nonpositive.
The second term appears exactly in the upper bound of (13), and is nonpositive by (23a) and f. S
f(x++1).
For the third term, geodesic convexity (which is guaranteed by the generalized geodesic convexity) between Tk+1 and x, leads to
f(xx+1) - -(xx)) (Grad f f(*++1), logmat σx) 1x-1 S 0.
The fourth term vanishes, by the construction of Lx+1- It only remains to show the bound on Qk, i.e., it remains to show:
{
Ax 8k+1
Be
Qk S min
si
|[Grad f(xx)1?
La+1
For 11, both sides equal zero by definition, verifying the inequality.
For 2 0, it suffices to show:
Bs 5k ||rrad (xxx
Qk +
<
if sk La+11 S
As 8k+1
Qx
||Grad (x <
if if 2 1.
L
First, we show the coefficient of the ||Grad f*x+11ll term is nonnegative. If st SaL++1 S 1, from (23d), we have
Ba+1 1 + CK+1 884 La+1 2 1
A B₁
Bx+1 SA+1)
+
8841 2 0.
A
A
Sy. Tx441 2 1, inequality (23c) directly ensures the nonnegativity of
Bet1 sk+ Tat11
An Then we use Lemma F.3 with Grad (*x) and w Grad f (I*++1). Considering the discrimi- nant of the quadratic form, our goal reduces to showing:
BR+1 + GR+1
Sk+ Lx+1
50
if sE 1,
Ak
(1
BR+1 GK+1
Bk o
(1 SxL411)
sk+1 BettIR+1/
Lk+1 S
if sf 2 1.
Au
As SK+1
Each case can be proved as follows:
stLk+1 S
From (23b) we have
Bk+1 + C++1
8k4 Lx41 SaLai1
A
Note that For arbitrary 61,62 E [0, 1]. if 82 - 8n we know
61- (1 62) s (1 61) (1 (1)))) = 82 50
We obtain the desired result by substituting &n 8 Ja++1 and 8,2 B2+1+
sct TaAt-
An
31

---

Published as a conference paper at ICLR 2026
S4Lk+1 2 1. From (23d) we know
BR+1 + Ck+1
+
Sk+1 L++1
A&
holds. Therefore, recalling that 28 S 0, we obtain:
(1
Bu+1 + Sk+1
SktLLAtI
Bs $ Lk+1
-k+11
Ak
As sk+1
Bg
Be+1 + Ck+1
sf.2241 2 12L22 (1 5xLR+1)
Ag
814
Ak
Now we move on to the case L6+1 0. Note that this implies 8; La+1 0 S 1. From Corol- lary C.2 we have Grad f(xx) L
Grad f(*++) 1). Then from the same calculation done to (22)
ZE+1
but without 1
Ezyi Grad f (zx) - T ad+1 Grad f (xx++1) I|[ and leveraging Lemma F.I as discussed
before for the first term we have:
V4+1
V& S
(s++IAx + 88+1 8k+2 (x++1))
+ Sft1(f - fo (Grad f ). Ze41 - -xa))
Ak (f(x*+1) f(xx) (Grad (Ix++1)); Tk41 xx))
(24)
Qk
S ( 8841 - 8142 ((x++1)) + Qk:
where
Qx = Ab ( (Ite,
Grad f(I++1), Grad f(Ix))
2 sh+ (Be+1 + S6+1) |Grad [( [(*++1)]]
2 5 Br ||rrad [(x))]
Note that following incquality follows from (23b):
8k-1 (B&+1 + Ck+1) - Ak S
Therefore, by leveraging the fact Grad f (xxx = Tok Øk41
Grad f (x++1))
1 Qu 2 sk+1 ($K+1 (Be+1 Sk+1) SxAA)) ||Grad f(x))1
1 T ||rrad fax))l ShSStt1A4
< 2 s Bu Grad 2
Substituting into (24), we obtain the desired result (13).
Proof of Corollary 4.3. From (13), we know that V4+1 Vx holds for all A: 2 ko- 6o-1. Consequently, we obtain
|log.r, _.Il 11 2Vx-1 S 2V R²,
Wk ko
which implies our desired statement. This concludes the proof.
D.2 PROOF OF LEMMA 4.4
Define S₃ min Lo
1
We prove the following by induction:
8 2
rS
Vk 2
(25)
From (7) and the lower bound r
B
from (6), we have
Bg
A
B
An
81 2 min
max
Ao:
Ao
min
Ao
B
Lo
L
min
To
= rSi
32

---

Published as a conference paper at ICLR 2026
Thus (25) holds for k 1. Assume (25) holds for k. show that it also holds for & + 1. Applying the induction hypothesis 5k 2 rSp and (10), we obtain the following inequality, which completes the proof of (25):
AR + As
AB
Bu+1
1
sk+1 min
min
TSk.
Ak
Bx+1
Br
An
La+1
min
rk,
rSS+1.
Tk++1
To prove (14), recall that we have I4 E BR(x.)) from Proposition 4.2 and To E Ba ( (x.).
Therefore, from Corollary C.2, we know that La S L for k 2 0, where L is smoothness parameter of f on (x. Fsd() Again, the compactness of this set guaranteed by Hopf- Rinow theorem. Therefore, we have S₃ 2 +. From (25), we conclude (14).
D.3 PROOF OF THBOREM 4.5
The statements in Theorem 4.5 follow from Proposition 4.2 and Lemma 4.4.
For notation simplicity, denote C2 min
As s6-1 Le+1?
Ba 2 5
First, summing up (13) from hi - to
have
V&+
Σ
Al 10:1A4-84420 8142 + +:11 f4) C, [|rad C f(:))) Vho-1 R.?
1
-
Combining with (x - f+) $ Vs, we have
C. Grad f(;))1 I 1 R². (26)
(441Ax(((IN) - f.) Σ (s,A,-1
1, + (f)((()))) - + f-)
Σ
2
k=
s=
Function value, last iterate.
Recalling (23a), we know that the summations (26) are nonnegative. Therefore, from (26), we have 1 Ax (fxx)) - f-) Y/2 R. Dividing both sides by 8k+1 A and applying Lemma 4.4, we obtain:
f(as) - S Ar
R² O
A4
Gradient norm square, last iterate.
As mentioned in Theorem 4.5, by applying Proposition 3.4 with y = Ty and r = x, to the above inequality, we have:
I =
T.2
|| Grad f(xx)]1 < R²
A
An
Function value, minimum norm selection.
Next, considering only the function value terms on the left-hand side of (26), we obtain
1 n22
8k+ Ac([(za) f.)
884 A; + 8;) (z) f.) S (MAA-
(27)
Ž=
Observing the sum of the coefficients of the ((xx) fx terms, we see that
r
Sk+1
Ag
- S + s;) = $1
(s
8
2
- k
L
No
T Aoo
where the second inequality is from Lemma 4.4. From (23a), we know 8; A for all 2 kn- Gathering these observations, from (27) we conclude
8141
A
L
L
min (f(x;)) f-) < 2rk
R² =
ie (ko.....))
k
33

---

Published as a conference paper at ICLR 2026
Gradient norm square, minimuam norm selection. From Lemma 4.4, we know Sk,, 8k+1 2 4. and from Corollary C.2, we know T 2 J. Focus- ing on the squared gradient norm terms in (26), enlarging the upper bound of the summation to stj. we obtain:
k, and recalling the definition C min
Au 8E42
Ta+1
1 [[]rad 2 R².
min (A, B,) Grrad f(x))l $
Σ
(28)
2L
Ko
k=
Note that, for simplicity, we have changed the lower bound of the summation to ky- Therefore, conclude
1
min ||Grad f(x;)]]
R². < co. Therefore, we
VE (k)
r Σi min {As: B.}
Lastly, from the second inequality (10), follows that lim sup
Ae
Br.
obtain A = O (min {Ax, B2 we conclude El= [i,
Σi A.
Point convergence to minimizer:
We follow the strategy of Suh & Ma (2025, Theorem 5). Take X' = (xr or is minimizer of f)- From the proof of Proposition 4.2, we have Ią Bn(x.), which is (relatively) compact by Hopf- Rinow theorem. addition, from (28), we have lim ||radd fxx)) 0. Hence, cluster points of {ax} } belong to X. Then, setting 4 = 2- - f.) + s BB ||Grad (xx)))
and using Lemma F.4 leads to the desired result.
D.4 PROOF OF THEOREM
Note that Theorem 4.1 is special case of Corollary 4.6. The condition BN B₁ + is clear from the definition, it remains to check that the parameter choice of Theorem 4.1 satisfies the conditions in Theorem 4.5.
is sufficient to check the conditions in (10). Since AR+1 - A₁ = a S 1, the first inequality of (10) holds. Also, A Bx+1 = + O, the second inequality if (10) holds. Now, observe that
As B&
BE+1
BR +1+ Be + a = 2 (a Be + 7)+ BB + +1+
=2 a a+5+ (1+5)
+
+
A
BB
B& B
B +1+1
decreasing for 2 0, since Bg is increasing. Hence, the third inequality (10) holds with
Aa
B Bi+2) Finally, we obtain:
+
Ag
1 min (A, B,) = = aaii+1) a(k + 1)(k + 2). 2
Σ
0
1=0
We obtain the desired result by Corollary 4.6.
D.5 PROOF OF COROLLARY 4.7
Proof of Corollary 4.7. We first establish that Vo < RR. Observe that f = V-1 B_1 is defined by Bn + Co- and by the definition, Ço = ((d(zo, x.)). Applying Proposition 4.2 with k 1 yields the result Vo < V_1 RR, required.
Next, we prove by induction that for all & Jb+1 BR(X.), G4+1 S To- and the inequalities in (13) hold.
First, consider the base case when
Since ((r), defined in (2), is nondecreasing in we have Ç(d(zo, x.)) < ((do). Moreover, we have d(a),æ,)² < 21% R²² < dg + (Bo + {( [Gradd (x0))] Hence, In E Br(xx), and, by the monotonicity of (()), (1
34

---

Published as a conference paper at ICLR 2026
x()( sC (V (Bo + ((do))& |Grad f(x0))] Su- With B₁ = Bi+ 2 +i+
the conditions of Proposition 4.2 satisfied with & 0. Applying Proposition 4.2 with l: 0 implies (13) holds with k = 0.
We now proceed to the inductive step. Assume that for some is 2 0 that Tk+1 E Bu(x,), S6+1 S Co and (13) hold for all 0 v k S kt. Then we have
2V. S R S d + (Bo + ((do))s || |rad fzo))l
I.) x - 2V
which implies 22+2 E BR(x.), and, again by monotonicity of <(-). Cā42 S So- Noticing B = R + So B 2k+2' and applying Proposition 4.2 with k = k shows that (13) holds with & = k + 1. This completes the induction.
Finally, following the same progression in the proof of Theorem 4.5, the desired result follows directly.
D.6 DETAILS RELATED TO REMARK 4.9
Proof of Corollary 4.8. Since limg->oe (Bx- B&) oc, there exists ho such that B₁ B₁ 2 holds for all k 2 ko. The statement follows immediately from Theorem 4.5.
Details related to Remark 4.9. Recall that the choices under consideration are A& = a(k: + +) +1, B₁ a(k + 1), and B₁ = a:(k + 1), with a E (0, 1] and ā E (0, We first check that
lim (Ba Ba) = (a a [k + 1) oc k
holds since á E (0,a), Therefore, (12) holds by Corollary 4.8.
The rest of the argument can be carried out in a similar manner to Appendix D.4, and we provide here for completeness. First, from a E (0, 1], it follows that:
A A& a
A Ba+1 = 1 2 0,
vk 2 0.
Next, since
A
B
B
1
1-a
A (1 BB
a +1
1 (a ā) + añ
+1
+1+
B&
A
Bb
As
decreasing for k 2 since both A₁ and Ba increasing. Therefore, (10) holds. Lastly, we know
min {A, B.) = Σ a(i
1 a ((k + 1)(k 2) (ko + 1) + 2))
h=
Therefore, by Corollary 4.8, we obtain
L
f(xx) - fo s 2r alk + 1) +1
R² =
k
and
2L
1
L 12
min ie (Ka.......k)) |[Grad fxx;)]] ar (k + 1) + 2) - (Ai + 1) + 2)
R²2
0
D.7 ADDITIONAL DISCUSSION ON AE. BR. AND Br
In this subsection, we discuss the respective roles and the design mechanism behind Ag. Be. and Br- As reminder, the sequences A₁, Bg, and B, are fundamentally designed to satisfy the inequalities in (23), thereby ensuring the validity of the Lyapunov analysis.
The convergence rates in Theorem 4.5 are written in terms of Aç. and a larger A guarantees faster theoretical convergence rate. However, the first inequality of (10) implies A I A S
35

---

Published as a conference paper at ICLR 2026
Au ++ 1, thus Aş cannot grow arbitrarily large, and its largest possible order is (kk). While the tightest choice A&+1 = A& + yields the theoretically fastest guarantee, the resulting step size is nonincreasing, since 5k+2 I Aut1
A2+1 Sk by (6) and (7). We therefore sacrifice part of the theoretical guarantee and consider a practically faster variant in the experiments in Section 5.
B& is introduced to account for the curvature of the manifold. It can be viewed as a variant of B& that upper-bounds Bg while incorporating an additional curvature term Ck. as seen in (12). The quantity B& appears in (6) and ensures that (23b), (23c), and (23d) hold, which governing 8 and validating the I Lyapunov analysis. A larger Bg leads to smaller step size, which makes these incqualities easier to satisfy, but smaller step size often results in slower convergence in practice.
the core inequalities
As a natural follow-up, one might wonder whether B₁ can simply be taken arbitrarily small so as to obtain a smaller BR, and hence a larger step size. However, if lim inf
B, co, then rk
As
may converge to zero due to the term Therefore, B₁ must have the same asymptotic order As, and Bg cannot be made arbitrarily small.
As
B which implies that the step size vanishes. A
E DEFERRED PROOFS FOR SECTION 5
E. DEFERRED PROOFS FOR SECTION 5.1
E.1.1 PROOF OF LEMMA 5.1
is sufficient to check A₁ 2(k + 2)0.1 + 2, Bg = 2(k + 2) and B₁ = 2(k + 2)0 + satisfy
the conditions in (10). First, by mean value theorem, for some & E [k 2, + 3] we have
A++1 - A₁ = 2 ((k + 3) (k+ 22) = 2 x 0.1 x -09
Since 1 + l 0.9 is nonincreasing function and i: 2 2, we have
A A S 2 x 0.1 x 2 <
Therefore, the first inequality of (10) holds. Also.
A Bs41 = 2(k 2) 2(k + 3)0 0.1 + (As+1 - A₂) +120,
the second inequality if (10) holds. Finally, the existence of E (0, 1) in the third inequality of (10)
follows from the facts
A
B
A BR
Bk+1
=1+11) - 0.
1
O. Vk 1 0,
lim
A
As
1.1 PROOF OF LEMMA 5.2
We show the generalized geodesic convexity of V and H. and then how combining these leads to the generalized geodesic convexity of F. For the definition of generalized geodesic convexity in BW (IR) space, see Definition B.9.
First, we show V is convex then V is generalized geodesically convex. To this end, fix any base (R"). Then, for all μ, v E BW(R") there exist optimal transport maps Tx and T
Let Z ~ From the convexity of V, one gets
VTT+,
(VV(T,,,,Z)) Tx, -T T (Z))
Take the expectation over ~ т. Then, using the fact that Wasserstein gradient of V = VV (from Definition B.8 and Santambrogio (2015, Section 8.2)), one gets the generalized geodesic convexity of with base x. Since the above argument holds for any π E (W(R"), we get the generalized geodesic convexity. For non-differentiable V, one consider VV as the weak gradient.
36

---

Published as a conference paper at ICLR 2026
The generalized geodesic convexity of H is established in Diao et al. (2023, Lemma 3.2).
Above two facts show that is convex then V and H. are generalized geodesically convex. Then, for any 7, U, E (W)),,
( 2 V(n) + (Gradaw V(n) o T T, TT, - T H(v) H(j)) + (Gradew H(p) T T TTT T T.
Summing these two inequalities lead to the generalized geodesic convexity of F.
E.1.3 PROOF OF PROPOSITION 5.3
Proof (i). We first show any C² (PBW(R))) function is always locally geodesically smooth, as stated in Theorem 3.3. However, since (W))) is not geodesically complete manifold, Hopf- Rinow theorem is not direct, and as a result one cannot directly apply Theorem 3.3. That said, note the only part the completeness of the manifold is used is when show the compactness of K := B(, D), which is the set that contains all the minimizing geodesic segment. For BW (IR"), we use an alternative set for B(K, D) which still guarantees the same properties (compactness, containing geodesics). Precisely, show the lemma below:
Lemma E.1. For given po. BWRR"), denote the minimizing geodesic μi to Hl as Vao.s Them, the set
K = U
Nu (10, 1])
is compact set in (R").
Proof of Lemma E.I. Define G: BW(R") x BW(R") x [0. BW(R") as:
G(uo: p1. 71 (t).
On the other hand, from Chen et al. (2019, Equation (9)), we know that any Bures wrres-Wasserstein geodesic between μo == (mis Mo) and μn == (m,, Σn} is (uniquely) defined by
G((mo, Mo), (m,, E)),t) (1 t)mo + tm1, 2 1/2 ((1 (1-t)22 1(2//
2 21 23122) 2
Leveraging the fact that the mappings > + Σー1 and + 1/2 are continuous on SPD(n), and that the binary operations used above (addition, scalar multiplication, matrix multiplication, etc.) are continuous, we know that G is continuous function. Since K K x [0, 1] a compact set,
K g (K K 10,
an image of a compact under continuous function G, and therefore is compact set. This concludes the proof.
Now, by taking Lx == P ||Hsss flp one can proceed in Theorem 3.3 and conclude that any C2 function on BW (WRR) is locally geodesically smooth.
Proof of (ii). Since we showed every C2 function on BW(R)) is locally geodesically smooth, it is sufficient to show that potential energy functional and entropy functional are C². Note the differentiable structure of BW(R") is determined by the differentiable structure IR" x SPD(n). Hence, by Theorem 3.3 it is sufficient to show these functions are C² SPD(n)) (and in fact they are C x SPD(n))).
First, for V. write the Gaussian density of N (m, Σ) as qm, ie..,
1
--(2 T - m)
1
x2)) (2 m/2 exp
Since smoothness is local property, consider (m, E), and some compact neighborhood K. We show V is smooth on U (K).
37

---

Published as a conference paper at ICLR 2026
As V is written by the expectation, we invoke the dominated convergence theorem, which constitutes the most technical part of the proof. First, since K is compact, there exist fixed constants > and 0 < mn S MM such that
mll - B, mK S Amiu (E)= S Amax (Σ) S Mr
V(m, Σ) E K.
Now, for any multi-index a = (am: ar) and any (m, Σ) E U,
Dm = (2) x Pa , m m))m.(())
(29)
with some Pa which is polynomial of the degree Pa S |aml + 22///- Now, recall for any (m, Σ) E K,
| S 1/mx, (del E)- S Amin (E)-n/2 S (mx) n/2
and Amax (-Σ-1) = X ( Mx
Writing y = 2 - m, for any (m, 2) EU C C.
(2) (mm,mm=m)))))))) S C])(()) (1+ [m,,)()) S C (=) (1+ [ axie llvll for some constant C > 0 that depends on mx. Mr, and a, but not and E. Next, we plug-in (16). This leads to
V (2)) smm xx)] C(1 + (1=11²) + )e 6 -znz llxll
1ell
Then, from ||m||| < B, we have (1 + ||=|²) < C(1 + |lyll)) for some C > 0. In addition, from [0, 2), there exists another constant C such that
<|= S Ce ae llul
Therefore, we get
Ve2) (.,mm(m)))))))) & C(1 |lvll
uxt wll =: G(vy)
where C' is the constant that does not depend on (m, 2).
Clearly, G(y) is integrable over C* from the fact that f Gly)dy is the Gaussian moment. Hence, G(y) is the integrable majorant of (m))aa for any (m, Σ) U.
Therefore, by dominated convergence theorem, for any |a| 2 0 and any (m, 2) EU
Do Vi,, Σ)
V(z) m. Imm sm))ddd
(30)
Rm
Once dominated convergence theorem is established, the rest is straightforward. Since xm,x(2) is smooth and Σ. Do (m) is well-defined. In addition, from (29) and the fact that all moments of the Gaussian are well-defined, the integral at RHS of (30) is well-defined. Hence, Da D (m, Σ) is well-defined for any multi-index This proves V E C (RR SPD(n)).
Proof of (iii). For the entropy, one observes the Gaussian entropy has the form of
H(m, Σ) = - -" (1 + log 2T) log det Σ. This function does not involve m and is clearly smooth on Σ E SPD(n). a result, H E C (R" SPD(n)).
n
(31)
2
Proof of (iv). Lastly. V. H being locally geodesically smooth trivially implies the local smoothness of F. For == V + H. and for any compact set K, let the local geodesic smoothness parameters of V and H are LK. L respectively. Then, for all 2,, y E K,
[[T GradF(r - Grad F(V)lI S I Grad V(4)= - Grad V(v) + [[ Grad H(n)) Grad H(v) L ) (2 v).
S k
One concludes by taking 1 I. I.HL
38

---

Published as a conference paper at ICLR 2026
E.1.4 PROOF OF COROLLARY 5.4
Algorithm 2 Bures-Wasserster Adaptive Gradient Variational Inference (BWAdaGVI)
Input: Log-potential V(X), mo R", Eo = el. so > 0. {Ak]s20 {Bk] A-1 0.
8 (0, 1). 2: for k 0, 0.1....d do
3:
m++1 my X [V(()),
NR NM : 23k)
S = Ex~mm (X)) 2x
Σ4+1 = (T - SR
) ) S1.
0k+1 = Ex-MEI1 [V()), MR+1 N(m Ex+1)
Sk = Ex-ARI1 [VV(X)] -
( - (ak, Sk) Sk)
Lx+1 = 2 F((U+1) F(uu) + (T (aR+1, Sk+1), (aR, S&))
A
As
AA
Bk+1
8
=
min
min
Av
Bk+1
8kk
B₁
+
As
L841
max; I..... (5
(32)
4: end for
We follow the proofs of Section 4, while incorporating the details needed to handle the geodesic incompleteness and specific geometry. We fact prove more general statement here: any that C²(BW(R"" and generalized geodesically convex satisfies the claimed convergence rate. Over the proof, we will write μ = (m, >))
First, from the transport map representation of Bures- Wasserstein space and Definition B.13, one can write Gradew (x) = V F(m, Σ) + (m,, [) - m)) Hence, Gradew F(u) = 2V 2vnF(m, Σ), and |/ Grad F(u)l,, 2 max;=1...... | (V Em))
Now, we check the new choices of stepsizes ensure the well-defined gradient iterates. First, observe
1-6
1
8/ S max I...., [A:(V Gradnw F(H))) | Gradnw F(uw)llop
Hence, 7 8% GradEw F(ux) > O, meaning T T (id 8% GradBw F(ux)) convex. By the Brenier theorem, this means the gradient iterate is the minimizing geodesic, and result the logarithmic map is properly defined as Tus..044 id. In addition, / - sy Gradpw F(us) implies it is invertible, that the exponential map well-defined.
Now, we set kn 0 and B₁ B& 1 for all & in the proof of Proposition 4.2. In addition, one sets Lyapunov functional by
1 T ThetiiN. id|l?
Ve 8 Ax F(P.)) 2 Be ||rradew F(ux)11 +
M411
This exactly coincides with (11) under (W( geometry. Then, Lemma D. with ( still valid, as it is just the rearrangement of the terms. Next, one checks whether Proposition 4.2 is still valid under the changed stepsizes. As in Proposition 4.2, we check the non-positivity of the terms appearing in Lemma D.1.
First term: One gets the non-positivity of the first term by the use of Lemma F.5.
Second term: As in the proof of Proposition 4.2, the stepsize rule (7) and F(u.) S
F(1x+1) gives the non-positivity.
The third term is guaranteed to be non posiiivie from the geodesic convexity, which is guar- anteed from the generalized geodesic convexity of F.
39

---

Published as a conference paper at ICLR 2026
To check the fourth term, notice
Gradow F(us+1) o TpgliiPE - Gradow F(uss) 1
1 2 F(P++1) F(uss) + Ss. ( (Gradøw F(U++1) T 1 Gradpw F(Hu))
LK-1
The denominator is non-positive, due to the generalized geodesic convexity. Rearranging the above inequality, the fourth term becomes non-positive.
One remains to check Qx < O. For Lk+1 * O, the sign discrimination test does not involve the explicit value Lx+1, so we have Qx S 0 for Lx+1 # O. For LK+1 = 0, since L++1 2 0 it implies Gradøw F(ux) Grad F(uk+1) Tun. by the use of Corollary F.11. Hence, one can proceed the same as in Proposition 4.2.
Above calculations prove that the conclusion of Proposition 4.2, ie., non-increasing Lyapunov func- tional, still valid for this problem.
Then, one is left with providing the lower bound of the stepsize, Le., analogous result to Lemma 4.4. To this end, we first observe that non-increasing Lyapunov functional gives
(2 p.) 2V_1 (B + 1) ||Gradww F(uo)ll + W.2 (No, N-) =: R??
Then, from Definition B.11.
[ - m R and dow (Ek, =tr==tEE + Ex 1 S RR. (33) Thus, for all mo. mo...., my are contained in Km = B (mm,, R) where the ball is with respect to Euclidean metric.
Next, we claim the condition on zw (2k,, E.) gives the uniform control over the largest eigenvalue of Σg. To check this, from Bhatia et al. (2019, Theorem 1),
tr(Ex + Ex - 2(x 2. E. /1/))///---
= min UEO(n)
|
--1///1
2 S
RR²
(34)
By triangular inequality, for any U E O(n).
Amaxx (Ex)
|E llop
11/2 - -
/
1 $111//
-1//2-_1!! E1/
+ |///
op
op
1/2 - =
+ V Amaxx (E,).
Since the above upper bound holds for any choice of U (n)), by taking the minimum over U E (n), one gets
Amax (Ex) S R+ v Amax (E.).
(35)
In addition, since we assumed the uniform lower bound the eigenvalues of Σx, we know for all i == ,....,k
Σ E Ky = I (R R V (I.) r].
Note Ky is compact (e.g., use Heine-Boue theorem with respect to Frobenius norm). Now, we construct the compact set
C (R"(R"):
K = (Km x Kx) U (mo, En)}
Then, K the compact set that contains all μu as well as Too Then, one gets K 5 K from Lemma F.8 and Lx from Proposition F.9.
Then, one repeats Lemma 4.4, with substituting I. by
T.x == max [Txx suppEK |/ Graduw F(1))
1
40

---

Published as a conference paper at ICLR 2026
and
1- 8
S min
Li
V
Gradaw F(ue)) (pp))lop
Note LK < oc from the compactness of K. The induction part is direct from the condition 1-ò 2 r. The uniform lower bound of Ss can be established by L& S T.K S Ix for k E N, which is from Corollary F. and suPkeN GradEw F(p SuppeR K || GradEw F(H))lop (1-8)KK-
Since we established the non-increasing Lyapunov functional well as the lower bound of the stepsizes, one can follow directly the proofs in Theorem 4.5 (i) and (ii) obtain the desired result.
For the point convergence result (Theorem 4.5 (iii)). one can take the above K as the (relatively) compact set containing (ur Since Lemma F,4 holds for general metric space, the argument is valid regardless of geodesic incompleteness.
We provide three theoretically ideal cases when ensure the uniform eigenvalue lower bound assumption.
Comments Remark 5.5
Case 1. Uniform upper bound the function value gives the uniform eigenvalue lower bound:
We show if sup;=1..... F(u)) MM for some M 0, then inf inf.=t,..... Amin (E:) 2 À for some
Ao > 0.
First, recall my B(m,, R) from (33). Since V is continuous, one has
Vini == min ((m) I ((m))
mc (m.,
uniformly for all N.
Since V is convex, from the Jensen's inequality and the explicit form of H in (31), one has
1
M 2 F(ux) E Ex) [V()) マ (1 log 2xx log det Ex
> V(mx)
7 + log22)) log det Ex
2 Vinin
z + log 2T) log det Σ
Then, one has
Amax (Ex)²-1 Amin ( (Ex) 2 det Σn 2 exp (2Vnin n11 log 2n) 2MM.
Using (35), one gets
Amin (Ex) exp (2Vmin - n1( + log 2nT) -2M) Aa
(R + v Ammax
which is the uniform bound independent of k.
Case 2. A good initialization gives the uniform eigenvalue lower bound.
Next, we show a good initialization gives also gives the uniform lower bound of the smallest eigen- values. Precisely, show if yn is initialized satisfy
R² + Bg + 1) [[Gradew F(wo))l W3 (No, u.) < Amin (Ew),
(36)
then there exists A which is the uniform eigenvalue lower bound.
From (34) and the standard inequality between the operator norm and Frobenius norm, there exists U E O(n) such that
|>//2 2//2 op
<
<|>
Σ!
1/2 F S R.
41

---

Published as a conference paper at ICLR 2026
Then, from Weyl's inequality singular value,
Mminn ( (=)²)) min ()) - | -
- R.
Since orthogonal matrix preserves the singular values, we have
Amin((2.) R S v Amin (2i).
Therefore, if R < V Amin (E.), we have the uniform eigenvalue lower bound by A = Vmmiu(E., (E+) R.)
Case 3. A priori upper bounds Amax (E.) and W2 (No, (4,) give the uniform eigenvalue lower bound:
Lastly, show if has the information about (max((.))) and (Ho, pa), then one sure the eigenvalue lower bound, though requiring much smaller stepsize. To this end, write R 2 R + Ammax(2), in (35). Note this quantity requires the priori upper bounds of Amax(E.) and W (Po: (+=). We precisely show if En ž + then Σx+1 also has the same bound.
To this end, we claim that by choosing 8 to be
2
E
(37)
117R
ensures the uniform eigenvalue lower bound. satisfies
Under this choice of 6, observe the stepsize rule
8k+1
(38)
117H
V Gradaw F(us++1llop
Unlike Case and 2, Case directly modifies the stepsize rule, so we need to check whether the proof of Corollary 5.4 is still valid. Fortunately, this is direct; since the proof of Corollary 5.4 works regardless of the choice of (for the condition - 6 2 r, since r is just the lower bound (10), one make the condition valid for given 6, with the cost of constant factor in the convergence rate), all the proof arguments will be valid. The only thing we are left is to show the above choice of ó ensures the uniform lower bound the eigenvalues of the iterates.
can choose smaller
We will follow the proof of Lambert et al. (2022, Lemma 6). One can consider their = in our setting. Then, by the same arguments, one can check Σx+1 is the generalized WasserWassrnniiin barycenter of
sk/ En [VVV and their (2v(()) as the normalized E& [X)]
P
1 1
Sk
8I. +2
(52.1+55) + 653
E
V (X)
IR
V
(X)
where
Es (VV(X)
E.
VVV))
Σ
21
Σ
27
E, V
VVV
|E
Next, suppose eT < Σg. Then, (35) is still valid (note this bound holds regardless of the existence of eigenvalue lower bound), we have
1 _ s A(==
E A(Ex)) RR
e es x < 4R
where the last upper bound from the same argument as in Lambert et (2022, Lemma 6). Next, as in Lambert et al. (2022, Lemma 6), we apply Altschuler et al. (2021, Theorem 1). Using the above eigenvalue bounds, one can check their a, B in Altschuler al. (2021, Theorem 1) for our case becomes
Ve
2
1
a =
3 =
+ 2R.
2
2
2c
42

---

Published as a conference paper at ICLR 2026
Then, by Altschuler et al. (2021, Theorem 1) if 2 Bun. VVX(x]] Amin (Ei++) 2 a/4. Since we would like to guarantee Amin (2i+1) 2 6, set a/4 = c. Then solving e, a, and with respect to gives
op I a/(23), get
4
13
a=
B
=
R.
9R
9R
2
Summarizing the above arguments, if Amin (Ek) 2 (9R)
and sy satisfies
2
S6 S
(39)
117R || (X)((
then Amin(2an (E++1) 2 (9R)-1. Since (9R) does not depend on k, we have the uniform lower bound. Then, since
V/V Gradew Flp || [VV(( 2 1lap 2 Bp [VV())]
the stepsize (38) satisfies (39). This means if we initialize Eo (9R) I. the new stepsize rule (38) ensures the eigenvalue lower bound E (9R) ( for all k 0.
In conclusion, all the proofs of Corollary 5.4 are valid under the choice in (37), and the uniform eigenvalue lower bound condition is guaranteed without any additional assumption. Hence, with 8 being (37), Corollary 5.4 valid without any additional assumption.
Special case when one can drop the uniform eigenvalue lower bound assumption Particularly, Case 3 gives the fully explicit algorithm which does not require any cigenvalue assumptions, one considers fixed mean Gaussian variational families with strongly convex potential V. This is because from Lambert et al. (2022, Equation (7)), if V is a-strongly convex, then one has
Amax (E.) S 1/a, W2 (()))))) S tr(Eo) + tr(E..) s (r((2)) + d/a.
Hence,
Te.² - (V + 1)2 [[Grad F(un)) tr(Ee) + d/a v1/a
which is fully explicit and known for any fixed Moo- One can plug-in this value to (37), which guarantees the convergence result of Corollary 5.4 without any additional assumption.
E.1.5 EXPLICIT FORM OF POTENTIAL AND GRADIENT
Lemma E.2 (Closed form of potential and SW-gradient). For any potential V R" of the form
1 - Q0 + -0+ Σ Gexp(aT0)+
V())
2
i=1
where Q > 0. E RR, C 2 O, and constant C. V is convex and satisfies (16). Moreover, the corresponding objective F(m, E) and its Bures- Wasserstein gradient Gradgw F(m, E) admit explicit closed-form expressions given by:
1 F(m, Σ) = = (1r(Q))) + T (m)) + q m +
Σ
ceep
X. x, [xTEX.) - log det E + C,
1=1
GradEw F(m, Σ)
Qm (Qm+ +m+ Σ cx exp x7 m xxTEX,
{
x
Q + Σ XXXx exp
x. + 2
1 ExTEX; - -
Proof. Convexity V follows directly from
V²V(0)= - - Σ ax,x exp(X 2
1=1
43

---

Published as a conference paper at ICLR 2026
In addition, V trivially satisfies (16).
For the closed form expression, the formula of the moment generating function of multivariate nor- mal distribution gives for any X; IR",
1 In addition, the expectation of quadratic form for any 0 with mean and the covariance Σ is explicitly known:
Nm.,) [exp(X? 0)] exp X x xxx2xx)
(40)
Eg Q0] tr(Qz) - m Qm. Using (40), (41), linearity of the expectation, and the closed form H. (31) yield the explicit form of F.
(41)
For WWggradient F, we have
(m, [Ea + a
ci Eo Σ
1=1
Ee-N(m, (V0
Σ
xxXx Eg (exp(X 0)].
1=1
Plug-in (40) to the above formula. Then, using Remark B.14 14 yields the claimed result.
F AUXILIARY LEMMAS
This section includes the auxiliary lemmas needed for Appendix C, and E.
F.I AUXILIARY LEMMAS FOR COMPLETE RIEMANNIAN MANIFOLDS
We first introduce an existing inequality concerning metric distortion on manifolds, whose proof can be found in (Zhang & Sra, 2016, Lemma 5).
Lemma F.1. Let Ç [0, x - [1, oc) be the function defined by (2). Then, for the vertices y, E M of a uniquely geodesic triangle, we have
log, pl
< ||log, pl1 C(d(p, x)) |log, vll² - (log, log, p)
(42)
We also introduce some basic property of the Riemannian logarithmic map for the intermediate calculation. For the proof, see Park et al. (2025, Lemma C.7).
Lemma F.2. For all x, y E M, let r be parallel transport from x to y induced from the minimiz- ing geodesic connecting and y. Then,
ry log,
log,
Next, we show the sign determination of the quadratic form by discriminant test is still valid in Riemannian manifold.
Lemma F.3. Let M be a Riemannian manifold, and p.q E M. Then, if a 2 O. c 2 0, and
b² 4ac < O, then
a ||el| + bT w),, + +cllw|l? 20, 0, Vo E T,M, w E T.M.
Proof. If one writes x := ||ulp = |1 then the standard discriminant test on real numbers gives
ax² - |b| ay + cy² 0. From the Cauchy-Schwarl inequality and the isometry of the parallel transport, for all IR
b 6T w), 2 |6| ||113011,, ||l|| |b] xy.
Therefore,
a Illl (Tzr, + cc|ww|l 2 ax - |b| ay cy²² 2 0.
44

---

Published as a conference paper at ICLR 2026
Next, provide the generalization of the result of Malitsky & Mishchenko (2024, Lemma 2) to metric spaces, which will be used to prove the asymptotic convergence to a minimizer for our algorithm.
Lemma F.4. Let (M,d) be metric space. Let A S M be the relatively compuct set, and X C M be some fixed set. Suppose {I (oR)KeN C A be a sequence whose cluster points belong to X and {ax} tER be sequence Rt. If (TR)KEN and {ax]ket satisfy
d (( + ax+1 S (T x)) + Ve E X,
then {aRIKeN converges to some element in X.
Proof. Let Let x² be any cluster points of {aihven. Then, since is relatively compact, there exist Jk. and Fky, two subsequences of (=*)keN= such that Fı. + x and #k, - x For fixed E X, the sequence d' (xx,, x) + as is non-negative and non-increasing, hence by the monotone convergence theorem the limit lim (d (Ix 2) + ax) exists. Then,
lim (d (xx, aa) lim (d (IR, x²) ax.) lim ak, 10000
10000
10000
lim (I x²) + ax2) = (I) lim Oky 1-30
1000
Hence, one gets lim;- 0k. = d + lim lim aky- By repeating the above argument sym- metrically using lim (d + on) lim (dr x²) + +ax,). one gets lim Ok₃ de (x x²) + lim ak. Thus, gets d² (x x²) and as result =x.
F.2 AUXILIARY LEMMAS FOR BW (IE)) SPACE
To prove Corollary 5.4, we need the similar results we obtained in Appendix F.I on BW(R") space. Despite the absence of the completeness, turns out that Bures-Wassen Wasserstein space have the similar properties. We provide the analogous statements for 2- 2-Wasserstein space and BW BW(R") discussed in Section 3 and Appendix F.1 that not require the completeness.
Lemma F.5. For any J., V, и E (T. and Ta,b being the optimal transport map from a the following inequality holds:
1[T2,x idl? S ||
idl; + 11 id] 2 (Trp id, Td,Tyn Tp,x
id)
The result trivially extends to BW(R") space.
Proof. Note T Te,p is also (possibly non-optimal) transport map from v to π. By the optimality of the optimal transport map T..
|Tu, - id]3 $ [|p Tuua - -idl @ |T Tp
[T, i. + 1|Tar id] (Tr. id, Tyix - id)
For (i), we used Proposition B.7 and Theorem B.2.
Note while we never used the curvature information to prove this, this result is indicating Ç = 1 for P2 (R²") space. In fact, P2 space is non-negatively curved space, so this is natural result. Lemma F.6. For all μ., V E (2aac(R))
have
(T u. - id) Tu., = - (T id).
Again, the result trivially extends to BW(R") space.
Proof. This comes directly from Theorem B.2 and Proposition B.7.
Next, we show the sign determination of the quadratic form by discriminant test is still valid in 2- Wasserstein space.
45

---

Published as a conference paper at ICLR 2026
Lemma E.7. For all u.v V E (R))) ifa 0. c 0, and b² 4ac S 0. then
a [|lll + b (v o nu w), + + /cllolly 2 0, Ve E P2,0 w E ToP2ac(R")
Proof. The same proof of Lemma F.3 works, once one simply substitutes the isometry of the parallel transport by Proposition B.7.
We now extend Proposition 3.4 and Corollary C.2 to BW (P*) space. Due to the fact that BW space is not geodesically complete, one cannot directly apply Proposition 3,4 to this setting. How- ever, turns out that one can still obtain the same conclusion under this specific geometry. Specifi- cally, observe that the geodesic completeness is used when verifying the compactness of the ball K Bee 3R), which is chosen to contain the point & the proof of Proposition 3,4 for all x, y K. We show one can still construct the similar compact set K for this problem, by using the explicit formula of tangent space, exponential map, and parallel transport.
Lemma F.8. Let K C BW(R") be compuct set, and F BW(R") + R. be a continuously differentiable function. Then, there exists Lo > O that only depends on K and F, and K 3 K that only depends on La, F, and such that for all L 2 Lo and all (mo, Ea); (m1, E K,
(m, Σ) = cm En) (-z Graddw F(m,, En) -
pmmm, GradBw F(mo: Eu)) K. (mo.Fo)
Proof. We first consider the space of O-mean Gaussians, so that can restrict the analyses to SPD(n).
Since K is compact in SPD(n) and the eigenvalue map is continuous, there exists Mx 2 0 such that for all K
< mK S Amin (2) S Amax (E) S MR oo.
Denote H == Gradaw F(E,) - - Gradew F(Eo). Since Gradaw F is continuous, there exists GK [0, oo) such that
sup
|Graddww F(E)]ls: = Gk-
yeK
Then, by triangular inequality and the isometry of the parallel transport (Proposition B.7)
H s
< |[Gradew F(Eo)) + [[Gradgw F(E,)) s 2G
Since || = HHHHH, one has
|| Mmk
2G Ck. /mk
H op S H
<
Now, we claim any Lo > Cg is the desired Lo. Define = 1 - -CK/LL Lo. which lies on (0, 1] from
the condition Ly > Ck 2 We then set
K = e 3E 4MK1}
We show K is the claimed compact set. This set is compact from Heine-Bort theorem on R with matrix norm (e.g. Frobenius). Then, from 1, K C K. Thus, we left show E E K. This is sufficient to show the eigenvalue bounds of Σ.
Now, for any L 2 Lu, since the tangent space is Sym(n), A == I H/L is symmetric. Using the explicit form of the exponential map (Definition B.12).
Σ = A
(2///2 A1) (>i/² AL)
Thus,
Amin (E)
amin (E//2
(Omax (2//2 Amax (5) 0(0maxx
46

---

Published as a conference paper at ICLR 2026
where a is the singular value. Now, using the inequalities between the singular values,
Omin 19/2AL) 2 miin )0 (AL) v Amin ( (At) 2 Mmxomin (AL), max 1// AA)== S maxx (2/// m = v Vmmax (amax(11)nn ( (At) & M (At).
Now, we bound min ( (Ar.) and (Ar))Ar)). Since Ar. is symmetric, the singular value o; [A(AAr) = X T./. Hence,
X
max; Ap
|| ||
La
minn (AL) 2 min
(1
26
L
L
L
L
A
max; A
|||||
Lo L
mmax (AL) E max
1
S 1 +
11
< ++
22
L
L
L
Aggregating the results, one gets
émp S Xmiul (E) Amax (E) S 4MK.
This means D E K, which completes the proof for the O-mean case.
The extension arbitrary BW is direct, since product of compact sets are compact. We iden- tified the compact set for matrix component, and can simply take the same BR), did in Proposition 3.4 for Euclidean mean part.
Proposition F9.. Suppose F BW(R") - R is generalized geodesically convex and locally geodesically smooth. For any compact set K, there exists LK >> such for all p., v K
1
F(u), T(.),TT.V id)
||rradaw F(v) Tp. Gradaw F(N)) 2 0.
2LK
Proof. First, we take Tug the same in the proof of Lemma F.8. Then, by Lemma F88, there exists a compact set K such that it contains all (m, Σ) (which is defined in Lemma F.8) for all L 2 Lo- Now, define
Lx == max (Lo, Lx}
where LR is the local smoothness parameter of the compact set K. Then, Lx is also lo- cal smoothness parameter of K, which is the compact set containing both v as well as æ exp, (- [Graddw F(v) ry Gradew F(u)) by Lemma F.8.
Now, one can proceed with the exactly same proof as Proposition 3.4, with the following substi- tutions:
Use μt instead of a. k,
Use Lx instead of Lig-
Use Definition B.9 for the definition of generalized geodesic convexity and local geodesic smoothness.
Use Lemma F.6 instead Lemma F.2.
Use Definition B.5, Proposition B.7, and Definition B.8 for the exponential map, logarith- mic map, parallel transport, and Riemannian gradient.
Remark F.10 (From P2.cc(R)) to (R"))) Note unlike the previous lemmas dealt with the gen- eral (R. from this proposition narrowed down focus to BW(R" specifically. This is because the compactness argument Lemma F.8 required the explicit form of the BW(R") geom- etry.
47

---

Published as a conference paper at ICLR 2026
Corollary F.11. Let F:
(W)R") + R. be generalized geodesically convex and locally geodesi- cally smooth function. Then, the following statements hold:
1. For any p.v E (W(RR) with F(v) F()) + (Gradew F(v),Tu., -id) 0, [Graddww F(v) o Tp Gradow F(1)) =0. =0.
2. Let K c. (R)) be compact Then there exists LK 0 such that for all μr, with F(v) F(p) + (Gradøw F(v), - id) + O,
K
||rradew F(v) o Tn GradEw F(u) S 2 F(v) - F()) + (GradBw F(v), Tu. id) I Lx.
Proof. One can proceed the same as in the proof of Corollary C.2, but based on Proposition F.9 Note dropped the condition (c, 3R) and N, since set N M (W(R).
G ADDITIONAL EXPERIMENTS
We provide additional experiments that will also validate the strength of RAdaGD over various tasks. All experiments were conducted either in free version Google Colab environment or Apple Macbook Air M4, and experiment took more than minutes
G.1 ADDITIONAL EXPERIMENTS ON GAUSSIAN VARIATIONAL INFERENCE
We conduct GVI experiments on two additional problems. Since both fall within the scope of Lemma E.2, BWAdaGVI accordingly endowed with the theoretical guarantees provided by Corol- lary 5.4.
1. π is a Gaussian target, corresponding quadratic potential V. In this case, V con- vex and T-smooth, so both BWAdaGVI and FBGVI (Diao et al., 2023) enjoy provable guarantees.
2. is the distribution of the response variable in log-link Gamma regression, generalized linear model. Although the canonical link for Gamma regression is the inverse function, the log-link is often preferred due the unfavorable behavior of the inverse (Czado, 2004, Page 21). In this setting, if
Y X; ind T shape = a, scale exp(X,T0) a
then the corresponding potential is
V = a (xxxp(-XTT0 xXT6.
Typically, the shape parameter a is either assumed to be known (e.g., = 1 corresponds exponential regression) or replaced by an estimate. For this experiment, we assume it is known. Notably, due to the exponential dependence on 0, this potential is not L- smooth. Consequently, FBGVI (Diao et al., 2023) lacks theoretical guarantee in this setting, whereas BWAdGVVI continues to provide one, as condition (16) is satisfied.
We chose the same Ag, Bg. and B& as in Section 5.1. For the Gaussian target task, we used V(z) } (I M'nwe 1 " (a- - Mune) for randomly generated Piree and 10-7 < Eine 3 / on R" x PSD(n) with n
This problem has 1-smooth target log density. The experiment results are aggregated in Figure 2. Note our method outperforms FBGVI even when the task is geodesically L-smooth. This is due to the fact that the adaptiveness allows to choose stepsizes at each step which may not be valid globally but valid locally.
For Gamma regression task. we used X; i.d.d
I.) and Y; X, and I(5, xTO for i = 1,...,50 50 and n. 25. Again, Figure 3 shows that our method outperforms FBGVI with various
48

---

Published as a conference paper at ICLR 2026
Staal ulpenwstit trocking
Oo-ooraancc out Faration.
135
2-Mtzperiinn citt Dh sestion
045
033
14
830
690
3003
Figure 2: Comparison between Diao et al. (2023) and our methods with different choice of Ak. Bg Gaussian target E (W).. Left: KL-divergence (KL(24 || Middle: Squared 2-
Wasserstein distance W (pu, T). Right: Minimum eigenvalues Amin (2&) over the iterations.
Smallest eggenveluee tracking
Germme GVI funciion
Facoin 1,50
FBGVI
FBGVI Is 1/12)
Raderrnn njk
12"
0.8
RAda00 0y
0.6
7042
02-
1B0Vils 3.01
TBOUl 125
F8CV/11 m 3/12)
10
M4duOn CId
02-
10
Laplace approsimation
a02
10'
200
400
500
830
Figure 3: Comparison between FBGVI (Diao al., 2023) and our methods Gamma regression target π with l = 50, = 25. Left: F(ux) Fmmn. where Fmin is the minimum value among all experiments. Right: Minimum eigenvalues Amin (Ex) over the iteration.
step size selections 1/3, 1/3,1/6, 1/6, 1/12. In particular, 1/6 was in fact the numerically optimal stepsize for FBGVI among our experiments. In addition, we observe that, in both experiments, the smallest eigenvalues remain uniformly bounded away from zero
G.1.1 GAUSSIAN VARIATIONAL INFERENCE WITH FISHER- RAO GEOMETRY
On the manifold of non-singular multivariate Gaussian distributions, there exists another natural Riemannian metric besides the Bures- wassersesnn metric: the Fisher-Rao metric. This metric arises from information geometry (Amari, 2016), and its geometric properties have been extensively stud- ied; see, for example, Andai (2007); Nielsen (2024).
Accordingly, one may consider RGD and RAdaGD equipped with the Fisher-Rao metric as alterna- tive algorithms for solving (15). At present, however, only RGD is practically implementable. The obstruction lies in the parallel transport operator under the Fisher-Rao geometry for multivariate Gaussians: is characterized by a differential equation whose closed-form solution not tractable (Han & Park, 2014; Lawson et al., 2023). As a result, implementing RAdaGD, which relies on explicit parallel transport, currently infeasible in this setting.
For the Gaussian variational inference (GVI) problem, the Bures waees-Wassrrttein geometry is theo- retically preferable to the Fisher-Rao geometry. The distinction is structural. The objective func- tional F(u) : DkL. (x..((l)) [T) exhibits more favorable properties under the Bures-Wasserstin metric- most notably, generalized geodesic convexity for any log-concave target density म, as established in Lemma 5.2. In contrast, under the Fisher-Rar Rao metric, strong negative result holds: for any smooth target distribution n, the functional F(u) = DRI. (μ ||*) is never geodesically convex when considered the Fisher-Rao manifold of all absolutely continuous probability measures on Rd with
49

---

Published as a conference paper at ICLR 2026
Poisson regression GVI furiction volue
10"4
20.5
70x
10-7
20'1
FaCV (x- 322)
PABaGD Citf
10-13
Laplace Approsinaaioo
Itemton
Figure 4: Comparison between Bures- Wasserstein gradient methods and the Fisher-Rao gradient method on Poisson regression target.
smooth densities (Carrillo et al., 2024, Theorem 3.3). Admittedly, this result does not exclude the possibility that geodesic convexity might hold after restricting the domain of F to Gaussian mea- sures. Nevertheless, it highlights the intrinsic difficulty of establishing geodesic convexity of the KL-divergence under the Fisher-Rao metric. In fact, where both μ and π are univariate Gaussian distributions, geodesic convexity already fails even in this highly structured setting.
be seen from the case in their proof
Despite the absence of convexity guarantees, our numerical experiments demonstrate strong empiri- cal performance when employing the Fisher-Rao metric for GVI. this appendix, we report results on the same Poisson regression task studied in Section 5. comparing the performance under the Mares-Wasserstehh and Fisher-Rao geometries. The results are presented Figure 4.
The strong empirical performance of the Fisher-Rao metric suggests that developing adaptive meth- ods under this geometry poyy--iillllllyy via approximating the parallel transport- may yield practi- cally effective algorithms. A systematic study of such methods, however, lies beyond the scope of this work. We therefore limit ourselves to reporting this empirical observation and defer a further investigation to future research.
We reiterate, however, that the Bures-Wassersiee Wasserstein metric provides substantially stronger theoreti- cal guarantees, most notably the (generalized) geodesic convexity of the KL-divergence, which is unlikely to hold in most cases under the Fisher-Rao metric. particular, among the methods com- pared in Figure 4, WWAdGGVI is the only algorithm that is supported by corresponding theoretical guarantee at this moment.
G.2 ADDITIONAL EXPERIMENTS ON THE SPHERE
We additionally demonstrate the numerical strength of our algorithm on Rayleigh quotient maxi- mization problem, standard benchmark task non-negatively curved manifold. While this prob- lem is not geodesically convex, RGD methods have shown strong empirical performance (Alimisis al., 2021; Kim & Yang, 2022). In this regard, we provide the numerical experiments on this problem. Given H E Sym(n), Rayleigh quotient maximization problem is written as follows:
1 - Hx. 2
min f(x) = #EBn
The above function f is geodesically Amax ()) Amin H -mmothh (Kim Yang, 2022, Proposi- tion 7.1), while not geodesically convex. Hence, we note neither RGD RAdaGI provides the theoretical guarantee. For numerical experiment, We considered n 1000 and H = ] (B + BT)
'This behavior is not entirely unexpected. Under a favorable initialization, the Fisher- Rao gradient flow of the -l-deieeeee exhibits some favorable properties than the corresponding Wasserstein gradient flow; see for example Lu et al. (2019; 2023); Domingo-Fnrich & Pooladian (2023).
50

---

Published as a conference paper at ICLR 2026
Rayleigh quotient -0.5x the iteration
10
10-1
RSD (s 7867977711966655577
10
10
Figure 5: Rayleigh quotient maximization. The stepsize was chosen by 1/L with L >max Amin(H)
where the entries of E Ruxa were randomly generated by the Gaussian distribution N (O, 1/n). This setup is the same setup as in Kim & Yang (2022, $7). For RAdaGD, we considered the same instance in Lemma 5.1. The result is displayed Figure 5, We observed numerical strength of our method.
G.3 ADDITIONAL EXPERIMENTS ON HADAMARD MANIFOLDS
While our primary motivation is the GVI problem, we also demonstrate the applicability of our method on Hadamard manifolds i.e., non-positively curved manifolds. Our algorithm is particularly suitable here: as discussed in Section 3, geodesic smoothness is highly restrictive on Hadamard manifolds, whereas the C condition remains general.
For the experiments on Hadamard manifolds, we consider two parameter choices: the benchmark parameters from Lemma 5.1, and an alternative set derived from Corollary 4.8 which can cope with unknown Ç.
Lemma G.1. Assume that č co. Then Av 2(k 2)0.1 B₁ 2(k 2) and Be aBr with a = 0.999 satisfy the conditions in Corollary 4.8.
Proof. The verification of (10) can be done in almost the same way as in Lemma 5.1, whose proof we provide in Appendix 1.1.1. The only difference that lim ( instead of }.-
Ax Bz
+
B
A = (2 +
-1
What remains is show that limg- (Br B&) oo, which is clear since
Be BB (1 - a) = 0.001 2(k + 2)0.1
k
00.
G.3.1 OPTIMIZATION ON MINNAATT FUNCTION
In this section, we focus the optimization problem on SPD(n). While Bures- Mures-Wasserstein (BW) metric provides one specific geometry on SPD(n), many algorithms lose its theoretical guarantee under BW metric due to the geodesic incompleteness of the space (Han et al.. 2021; Thanwerdas & Pennec, 2022a; Thanwerdas, 2022). On the other hand, there is another geometry on SPD(n) called Affine invariant metric (Fillard al., 2005; Pennec al., 2005). This metric is defined by (A, B) := |log A A1/22 -1/2| and (S, R2),
SA---SAA-1 R) for S, R Sym(n). This SPD(n) with geodesic completeness and the curvature bound |11/2, 0] (Crisciticllo Boumal, 2023a). There is a family of functions under this geometry where the assumptions of our methods are met.
metric induces non-posity curved geometry
51

---

Published as a conference paper at ICLR 2026
Proposition G.2. For any C² convex function & - R. the functional fx(X) (log det X) is
locally geodesically smooth and generalized geodesically convex on (SPD(n), daI).
Proof. First, we consider g(X) log det X. We first show for any arbitrary base and X, Y SPD(n),
g(V - x Grad ((X), logg Y logg x)
(43)
To verify this, first observe that Grad G()) X. This can be verified by the definition of Rieman- nian gradient. Riemannain gradient is defined by the operator satisfying for all H E Sym(n)
tr(x-1 Grad (XX)x H) = (Grad ( H)x = dgx = tr(x-1 H).
For the last inequality used the well-known formula for the derivative of log log-determinant func- tion. Since X E SPD(n). this implies Grad g = X.
Next, we show x B. To check this, from the definition of the parallel transport in (SPD(n), dai) (Nguyen, 2022, Supplement 1.1),
r = -BX- ( =B.
Hence,
; L,X,log log Y log x) = (B, Y) (B.logn X).
Now, for any matrix M, from the definition of the Riemannian metric on daI and logarithmic map (Nguyen, 2022, Supplement 1.1),
(B, log M)
tr(B log M) = tr(BB 1/2 1BB "B
log (B 1/33 log det B log det M.
Therefore,
( X, logg Y -logg x) = llog det Y log X = gg(Y) x(X)
which shows (43).
Now we are left with the & part. By the convexity of Φb,
( - gBB) 2 (o(g()))(((A) g(B)) ('(( (Th Gradg(B), loge A -log B) (rÉ Grad( - (B), logo A - logc B)
where the last equality is from the chain rule and linearity of the parallel transport and gradient. This completes the proof.
For the local geodesic smoothness, since (SPD(n), dai) is the smooth and geodesically complete Riemannian manifold (Thanwerdas & Pennec, 2022b), Theorem 3.3 is applicable. As discussed in the proof of Proposition 5.3, the differential structure of the manifold is invariant with the Rieman- nian metric, and one gets the desired result by the smoothness of log (-). Since á is C2 and log and det are smooth on (pD(n), fo is C².
log det X coincides (up to an affine transformation) with the Gaussian entropy H. discussed in Sec- tion 5.1. Hence (X) corresponds to the entropy matching or regularization. Such functions frequently appear in composite objectives the penalty term. While the full composite may not satisfy our thcoretical assumptions, the fx( part does, and our method has shown strong empiri- cal performance in these cases.
Accordingly, report three experiments: two directly optimizing f(X) = |log det X - cc and fXX |log det c| for some c e I, which fully meet our assumptions in Theorem 4.5, and the other optimizing the composite function f(X) = tr(XC) log det X for some C' E SPD(n) which, despite weaker theoretical guarantees, has important practical implications as generaliza- tion of linear semidefinite programming as well as the log- lig-likelihood of Wishart distribution (Boyd & Vandenberghe, 2004; Wang al., 2010; Han et al., 2021). Note llog det X cl is geodesically 2n-smooth, while the other two functions are not geodesically L-smooth. The experiments are sum- marized in Figure 6. Our method achieves either comparable superior performances and shows robustness to step-size choices.
52

---

Published as a conference paper at ICLR 2026
13:11
an ud
19"1
1370
10:3
Figure 6: Comparison between RAdaGD and various stepsizes RGD on SPD(n) for different tasks. We set n = 50. We chose the RGD stepsize for non L-smooth tasks by the quantities that make the RGD iterates stable. Left: f(X) = [log det X - 1|2 Middle: f(X) = [log det X - 11² Right: f(X) = (r(XC) log det X.
SPD space..ccc.. dim- 10, sample 1000
Hyperbbllc space, dim 50, sampre 1000
10
40
10
10
10
10
10'
10
10
10
60 (5 1/2200)
GB 8- 1/1100)
10.4
RGD 1.1000)
RGD (s 1,20000
10
5 1/550.00
31-16
RCD 1/1000 01
C-3
BAdaGO 03-11 (- Hone
12
RAdeGO ov (-None
Iasation
Figure 7: Comparison between RGD with different stepsizes and our method on Fréchet mean estimation problem under different geometries. We chose L that makes the RGD iterates to well- defined. Left: SPD(10) with 10000 points. Right: H with 1000 points.
G.3.2 FRECHET MEAN ESTIMATION
We show the validity of our method on Fréchet mean estimation (or barycenter) problem, which is standard benchmark in Riemannian optimization literature (Alimisis et al., 2021; Kim & Yang. 2022). Fréchet mean estimation problem is widely studied problem in non-Euclidedd statistics. In particular, Fréchet mean imposes favorable statistical properties on non-positively curved spaces (Sturm, 2000; Bačák, 2014b; Yun & Park, 2023; Brunel & Serres, 2024; Kim et al., 2025a). For hadamard manifold M and for given {pr..... .P.] M, the empirical Fréchet mean over p;
p. := arg min f(x) = ZEM
(l,
This problem is known to be geodesically n-strongly convex (Kim & Yang, 2022, Proposition 7.2), while not generalized geodesically convex. On the other hand, this function is C² and thus locally geodesically smooth, while not geodesically L-smooth. Hence, neither RGD nor RAdaGD guaran- tees the convergence for this problem.
In this regard, we conduct experiments on Fréchet mean estimation problem on Hadamard mani- folds. We analyze two experiments: one is SPD(n), which was what we used in Section 5, and the other is the hyperbolic space. Hyperbolic space is widely used in machine learning for hierarchical representation learning (Nickel & Kiela, 2017), and as a result Fréchet mean on this space is also one of the standard benchmarks (Alimisis et al., 2021; Kim & Yang, 2022). The results are presented in Figure 7. While not outperforming, our method still comparable to the standard RGD on both tasks.
53

---

Published as a conference paper at ICLR 2026
H
LLM USAGE
This paper made use of Large Language Models (LLMs) primarily for editing grammar, spelling, word choice, and sentence clarity. In addition, LLMs were used occasionally to help with minor coding tasks during the experiments.
54


Reviews:
[{"review_id": "TUhXTgkG1v", "reviewer": "Reviewer_jiEp", "summary": "The authors propose adaptive gradient descent methods on Riemannian manifolds. Following the approach introduced in Suh and Ma (2025), they show that the iterates of their algorithm converge in objective value at a rate of $\\mathcal{O}(1/k)$ when the potential function $V$ is locally geodesically smooth and generalized geodesically convex. They apply this framework to Gaussian Variational Inference (GVI), establishing provable computational guarantees while relaxing the standard $L$-smoothness assumption.", "questions": "1. The citation to \\emph{Ansari-Onnestam \\& Malitsky (2025)} should be removed, as that paper has been withdrawn.  \n\n2. Some discussion should be included for non-convex objectives, clarifying whether the proposed methods or convergence results extend in any way beyond the geodesically convex case.  \n\n3. Since only geodesic convexity is assumed, there may exist multiple minimizers. It would be helpful to clarify whether the choice of $x^*$ in the statements of results (e.g., Theorem~4.1) affects the conclusions or constants.  \n\n4. In Section 4.2, the authors address the case of manifolds with negative sectional curvature and unbounded $\\bar{\\zeta}$ by assuming that $d(x_0, x_*)$ is bounded above (in Corollary 4.7). This condition appears vacuously true, since otherwise the constant $R$ defined in Theorem 4.1 would be infinite, implying that the smoothness parameter $L$ should be infinite as well. It is unclear what substantive difference is between Corollary 4.7 here and Corollary 4.6 in the previous section. In particular, as suggested by Corollary 4.7, what happens if we set $\\tilde B_k = B_k + \\bar{\\zeta}_0$?  \n\n5. Although the proposed method is a direct extension of Suh and Ma (2025)  to the manifold case, more intuition and explanation should be provided regarding Algorithm~1. For instance, what are the respective roles of $A_k$, $B_k$, and $\\tilde B_k$ in the algorithm? Is there an intuitive interpretation of these quantities or of the update mechanism?  \n\n6. In line 326, the sentence  \n    ``The following corollary states that, as a trade-off, if $\\zeta$ is known to be bounded, we can still achieve the same asymptotic convergence rate without any prior knowledge of $\\zeta$.''\nis unclear. What does $\\zeta$ refer to here? Is it the function defined in equation~(2), or is this a typo?  \n\n7. In Section E.1.3, the proofs of parts (iii) and (iv) are provided, although these statements do not appear in Proposition 5.3. It would be clearer to reorganize the proof into Steps 1--4 instead. In Lemma E.1 (line 1832), the phrase should read:  \n    ``is a compact set in ${\\rm BW}(\\mathbb{R}^n)$.''", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "Overall, the paper is carefully written and mathematically rigorous. It extends the work of Suh and Ma (2025) to the setting of optimization on Riemannian manifolds and establishes convergence guarantees for the proposed algorithm. In particular, this work provides a algorithm with convergence results for Gaussian Variational Inference (GVI) even when the target distribution is not log-smooth.", "weaknesses": "The main idea appears to be a rather direct extension of Suh and Ma (2025) to the Riemannian setting, with additional technical complications arising from the non-Euclidean geometry. In particular, Algorithm 1 in the paper is largely a manifold-based translation of Algorithm 2 in Suh and Ma (2025). The comparison with the reference paper is somewhat limited and does not sufficiently emphasize the novel challenges or technical difficulties encountered in extending the analysis to the manifold case.\nIn Section 4.1, the authors mention that the constants $\\tilde B_k$ are introduced to handle the terms $\\zeta_k$, which are specific to the manifold setting (since $\\zeta_k \\equiv 1$ in Euclidean geometry). A more detailed explanation and comparison in this direction would be appreciated—especially to clarify how these constants affect the algorithm’s behavior and the convergence analysis.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "C3z1IU0z3g", "reviewer": "Reviewer_LVsh", "summary": "The paper proposes RAdaGD, a line-search-free family of adaptive gradient descent methods on general Riemannian manifolds. The method adapts the stepsize by estimating a local smoothness quantity from two consecutive iterates and gradients, and it is analyzed under local geodesic smoothness and generalized geodesic convexity, which are weaker than the usual global geodesic $L$-smoothness. The authors prove a non-ergodic convergence rate\n$$\nf(x_k)-f(x^\\star) \\le O(1/k)\n$$\nmatching classical Riemannian gradient descent with known $L$, and also give bounds on $|\\operatorname{Grad} f(x_k)|^2$. A technically involved part of the analysis is the treatment of curvature via the parameter $\\zeta(\\cdot)$ and the triad $(A_k,B_k,\\tilde B_k)$ in the Lyapunov function. As a main application, the paper shows how to instantiate RAdaGD for Gaussian variational inference on the Bures–Wasserstein manifold, obtaining (to the authors’ knowledge) the first convergence guarantee for GVI without assuming global $L$-smoothness of the log-density, under an additional eigenvalue-boundedness condition. Experiments on Poisson regression GVI and negatively curved spaces demonstrate competitiveness with prior Riemannian and BW-based optimizers.", "questions": "1. Can the authors clarify whether the adaptive rule for $s_{k+1}$ can be made compatible with retractions/vector transports (as opposed to exact $\\exp$ and parallel transport) while preserving the rates, possibly with higher-order error terms?\n2. In the BW-GVI application, the extra step-size clipping in equation (17) depends on the largest eigenvalue of an expectation involving $\\nabla^2 V(X)$. How is this implemented in practice when only Monte Carlo estimates of expectations are available, and how sensitive is convergence to this clipping?\n3. The convergence proof in Section 4 relies on the compactness/boundedness of the iterates (via the Lyapunov argument). Is there an example of a nonnegatively curved manifold and locally geodesically smooth, generalized geodesically convex function where RAdaGD would fail without this boundedness step?\n4. For the Poisson regression GVI example, can the authors provide a small ablation showing RAdaGD vs. a Riemannian method with backtracking line search, to separate the benefit of adaptivity from the benefit of avoiding line search?\n5. The paper mentions possible stochastic variants. Are there obstacles in controlling the estimate of $L_{k+1}$ in (5) under stochastic gradients, or is it mainly a matter of concentration bounds?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* A concrete, implementable Riemannian adaptive gradient method (RAdaGD) that removes line search and still attains $O(1/k)$ in function value.\n* Relaxation from global geodesic $L$-smoothness to local geodesic smoothness, supported by the statement that every $C^2$ function on a complete Riemannian manifold is locally geodesically smooth; this enlarges the class of admissible objectives.\n* Careful curvature-aware analysis via the $\\zeta(\\cdot)$ parameter and the triplet $(A_k,B_k,\\tilde B_k)$, including extensions to unbounded or unknown curvature.\n* Application to Gaussian variational inference on the Bures–Wasserstein manifold, yielding (under an eigenvalue assumption) what appears to be the first convergence guarantee for non-$L$-smooth GVI.\n* Empirical results on Poisson regression GVI that illustrate the claimed setting where prior methods lack guarantees but RAdaGD still converges.", "weaknesses": "* The core convergence theorems require generalized geodesic convexity, which is stronger than plain geodesic convexity and may limit applicability on general manifolds or nonconvex statistical objectives.\n* The method still assumes access to exponential maps and parallel transport on the manifold; this is reasonable for BW$(\\mathbb R^n)$ but can be a practical limitation on more general manifolds where only retractions and vector transports are available.\n* The step-size rule in Algorithm 1 depends on quantities that must be computed from two consecutive gradients and parallel transports; while explicit, it is more involved than standard RGD and the paper could more clearly quantify this overhead.\n* Experiments are focused on a few GVI scenarios and do not yet explore larger-scale or stochastic-gradient settings, which are natural for VI in practice.\n* The GVI convergence result needs an additional uniform lower bound on the covariance eigenvalues, which is plausible but not guaranteed by the algorithm itself.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "jhTzLI9oJa", "reviewer": "Reviewer_j1ST", "summary": "This paper introduces RAdaGD, a novel family of adaptive gradient descent methods on Riemannian manifolds. The method adapts the step size automatically—without explicit line search—by estimating a local smoothness parameter, achieving the first known non-ergodic convergence rate \\mathcal{O}(1/k) under local geodesic smoothness and generalized geodesic convexity. The authors extend adaptive optimization concepts from the Euclidean setting (e.g., Suh & Ma, 2025) to Riemannian manifolds via nontrivial generalizations of Lyapunov analyses, parallel transport arguments, and curvature-dependent scaling factors.\nA notable application is given to Gaussian Variational Inference (GVI), where the method provides the first convergence guarantee without assuming global L-smoothness of the target log-density. Empirical experiments confirm the theoretical results and show competitive or superior performance compared to existing algorithms.", "questions": "Please comment about the weaknesses", "limitations": "", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The paper is rigorous, technically elegant, and addresses a clear theoretical gap in Riemannian optimization—namely, the lack of adaptive step-size methods with formal convergence guarantees. The introduction of local geodesic smoothness as a relaxation of traditional L-smoothness is both natural and powerful, broadening the class of admissible objective functions while preserving tractability.\nThe convergence proofs are carefully executed and extend nontrivial elements from Euclidean adaptive methods to manifold geometry (e.g., dealing with exponential/log maps, curvature terms, and parallel transport). The theoretical guarantees (non-ergodic \\mathcal{O}(1/k) rate) are significant, and the paper connects the framework elegantly to practical applications in Bayesian inference. The exposition is also exemplary—clear definitions, consistent notation, and intuitive motivation for the introduced assumptions.", "weaknesses": "While the paper fills a critical theoretical gap, I have a few questions regarding its positioning and broader implications:\n\t1.\tRelation to recent Riemannian adaptive methods.\nThe paper positions itself as the first of its kind for adaptive algorithms on manifolds. However, a recent arXiv preprint (link https://arxiv.org/abs/2306.16617) also studies adaptive Riemannian optimization with curvature-aware preconditioning and adaptive metric scaling. Could the authors clarify whether RAdaGD fundamentally differs from these approaches (e.g., in the use of local smoothness estimation versus curvature-adapted preconditioning)? Are the convergence guarantees directly comparable, or do they apply to distinct settings (e.g., deterministic vs. stochastic, local vs. global smoothness)?\n\t2.\tTightness of assumptions.\nThe local geodesic smoothness assumption, while weaker than L-smoothness, still presupposes control of curvature-dependent Lipschitz constants L_K. Is it possible to relax this further, for example to Hölder-smooth functions, or to obtain adaptive control of curvature-dependent constants during optimization?\n\t3.\tEmpirical scope.\nThe experiments are clean but limited to Gaussian Variational Inference tasks. Given the generality of the method, could the framework extend naturally to Riemannian neural optimization or hyperbolic embeddings (where curvature varies strongly)? A short experiment in such a context could enhance the impact.\n\t4.\tAdaptivity mechanism.\nThe paper introduces adaptive sequences A_k, B_k, \\tilde{B}_k, but it would be helpful to clarify whether the adaptivity primarily tracks local smoothness, curvature, or both. In other words, is the method “geometry-aware” in practice, or purely driven by gradient variance?\n\nOverall, these are minor concerns—the theoretical contribution is strong and original.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "Oljrpcj2oP", "reviewer": "Reviewer_Jjvv", "summary": "This paper proposes a adaptive optimization method, RAdaGD, for smooth optimization problems defined on Riemannian manifolds. The authors extend adaptive gradient methods to the non-Euclidean setting under relaxed smoothness and convexity assumptions—specifically, local geodesic smoothness and generalized geodesic convexity—and provide a non-ergodic $\\mathcal{O}(1/k)$ convergence rate. The approach is further applied to Gaussian Variational Inference (GVI) problems formulated on the Bures–Wasserstein manifold of Gaussian distributions. However, this paper heavily relies on the algorithms and theories presented in Suh & Ma (2025). Furthermore, several theoretical assumptions and proof steps require clarification or stronger justification.", "questions": "1.\tThe main theoretical results (e.g., Theorem 4.1) rely on an upper bound $L \\geq L_k$ for $k \\geq 0$. While the authors address the existence of such a bound, this assumption obscures the distinction between global and local smoothness, thereby reducing the theoretical appeal of the proposed framework. \n2.\tThe constants in the $\\mathcal O{(1/k)}$ rate (e.g., dependence on $r, L, R, \\bar \\zeta$) should be discussed. If these constants can be large in curved manifolds, the theoretical rate may not be meaningful in practice.\n3.\tTheorem 3.3 claims that all $C^2$ functions are locally geodesically smooth. While intuitively true, the proof relies on the boundedness of the Riemannian Hessian. The manuscript should more explicitly state the curvature-dependence of this claim.\n4.\tThe authors mentioned another adaptive gradient descent algorithm (Ansari-O¨ nnestam & Malitsky (2025)) on a Riemannian manifold in the previous text , but did not conduct a comparison in the experiments. The authors claim that the paper was withdrawn from arXiv on September 15, 2025, yet its V1 remains accessible for viewing.\n5.\tIn Corollary 5.4, the author claims that $\\bar \\zeta =1$, based on the assumption in Corollary 4.6 that $\\bar \\zeta =sup_{x \\in N} \\zeta(d(x,x_*))$, but I have not found a proof for $\\zeta = 1$.", "limitations": "", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1.\tThe authors propose a family of Riemannian adaptive gradient descent algorithms, RAdaGD, which automatically adjusts the step size by approximating the local smoothness parameters. \n2.\tThe introduction of local geodesic smoothness and generalized geodesic convexity provides a new theoretical framework that generalizes existing smoothness/convexity assumptions.\n3.\tThe application to Gaussian Variational Inference under the Bures–Wasserstein geometry is well motivated and could be useful in practice.", "weaknesses": "1.\tThe main theoretical results (e.g., Theorem 4.1) rely on an upper bound $L \\geq L_k$ for $k \\geq 0$. While the authors address the existence of such a bound, this assumption obscures the distinction between global and local smoothness, thereby reducing the theoretical appeal of the proposed framework. \n2.\tThe constants in the $\\mathcal O{(1/k)}$ rate (e.g., dependence on $r, L, R, \\bar \\zeta$) should be discussed. If these constants can be large in curved manifolds, the theoretical rate may not be meaningful in practice.\n3.\tTheorem 3.3 claims that all $C^2$ functions are locally geodesically smooth. While intuitively true, the proof relies on the boundedness of the Riemannian Hessian. The manuscript should more explicitly state the curvature-dependence of this claim.\n4.\tThe authors mentioned another adaptive gradient descent algorithm (Ansari-O¨ nnestam & Malitsky (2025)) on a Riemannian manifold in the previous text , but did not conduct a comparison in the experiments. The authors claim that the paper was withdrawn from arXiv on September 15, 2025, yet its V1 remains accessible for viewing.\n5.\tIn Corollary 5.4, the author claims that $\\bar \\zeta =1$, based on the assumption in Corollary 4.6 that $\\bar \\zeta =sup_{x \\in N} \\zeta(d(x,x_*))$, but I have not found a proof for $\\zeta = 1$.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "dirNbTv11s", "reviewer": "Reviewer_eAPu", "summary": "The paper proposes RAdaGD, a family of adaptive gradient descent algorithms on Riemannian manifolds. The method generalizes Euclidean adaptive step-size rules to the manifold setting, avoiding the need for a known smoothness constant or line search. The authors prove non-ergodic O(1/k) convergence under local geodesic smoothness and generalized geodesic convexity, relaxing the assumption of global smoothness. They apply RAdaGD to Gaussian Variational Inference (GVI) and provide convergence guarantees even when the target log-density is not globally smooth.", "questions": "1. Can the authors clarify how restrictive the local geodesic smoothness assumption is in practice? For example, are there realistic manifolds or objective functions where it might fail to hold, and how would RAdaGD behave in such cases?\n\n2. As I've understood it, the experiments are limited to Poisson regression on the Bures–Wasserstein manifold. Could the authors provide results or discussion on how RAdaGD performs on other manifolds (e.g., SPD or hyperbolic spaces) or compare it directly to Riemannian Adam or stochastic methods?\n\n3. How sensitive is RAdaGD to the choice of initial step size or curvature estimates in practice, and what is the computational overhead of computing the adaptive step-size parameters compared to standard Riemannian gradient descent?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The paper makes a novel theoretical contribution by defining local geodesic smoothness and proving that it holds for all twice-differentiable functions on complete manifolds, thereby relaxing the need for global smoothness assumptions. It introduces generalized geodesic convexity as a principled framework for analyzing adaptive Riemannian optimization methods and provides the first adaptive Riemannian gradient algorithm with non-ergodic O(1/k) convergence rates. The application to Gaussian Variational Inference is well motivated and theoretically sound, extending results from the Euclidean to the non-Euclidean setting. Overall, the theoretical analysis is rigorous and carefully adapts adaptive step-size techniques to curved spaces.", "weaknesses": "The paper’s main weaknesses lie in its presentation and empirical scope. The exposition is mathematically dense and at times difficult to follow, with heavy notation and definitions that could be streamlined for clarity. Empirical validation is limited to a single Poisson regression experiment, which does not fully demonstrate the algorithm’s practical benefits or robustness across different manifolds. Comparisons with other Riemannian adaptive or stochastic methods, such as Riemannian Adam, are missing, leaving open questions about relative performance. Some theoretical assumptions, including bounded curvature and known constants, may be hard to verify or implement in practice, and the computational cost of adaptive step-size evaluations is not thoroughly discussed. Overall, while the theory is strong, the work would benefit from clearer exposition and more extensive empirical evidence.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)