A Framework for Bilevel Optimization on Riemannian Manifolds
Authors: Andi Han
Bamdev Mishra
Pratik Jawanpuria
Akiko Takeda
Venue: NeurIPS 2024 poster (2024)
Keywords: Bilevel optimization
Riemannian optimization
Riemannian hypergradient
Primary area: optimization_for_deep_networks

Abstract: Bilevel optimization has gained prominence in various applications. In this study, we introduce a framework for solving bilevel optimization problems, where the variables in both the lower and upper levels are constrained on Riemannian manifolds. We present several hypergradient estimation strategies on manifolds and analyze their estimation errors. Furthermore, we provide comprehensive convergence and complexity analyses for the proposed hypergradient descent algorithm on manifolds. We also extend our framework to encompass stochastic bilevel optimization and incorporate the use of general retraction. The efficacy of the proposed framework is demonstrated through several applications.

A Framework for Bilevel Optimization on Riemannian Manifolds
Andi Han' Bamdev Mishra? Pratik Jawanpuria' Akiko Takeda'
RIKENN AIP Microsoft, India UUiversity of Tokyo
andi.han @ rik
(bamdevm, jawanpurrll @ microsoft.com.. takeda mist..a-tokyy..a..ppp
Abstract
Bilevel optimization has gained prominence in various applications. In this study, introduce framework for solving bilevel optimization problems, where the variables in both the lower and upper levels are constrained on Riemannian man- ifolds. We present several hypergradient estimation strategies manifolds and analyze their estimation errors. Furthermore, we provide comprehensive conver- gence and complexity analyses for the proposed hypergradient descent algorithm on manifolds. We also extend our framework to encompass stochastic bilevel optimization and incorporate the use of general retraction. The efficacy of the proposed framework demonstrated through several applications.
Introduction
Bilevel optimization is hierarchical optimization problem where the upper-level problem depends on the solution of the lower-level, i.c.,
min F(x) f(x,y ZER
s.t. yd (z) arg gin y).
yERdy
Applications involving bilevel optimization include meta learning [16], hyperparameter optimization [8]. and neural architecture search (NAS) 5 name few. The lower-level problem is usually assumed to be strongly convex.
Common strategies for solving such problem can be classified into two categories: single-level refor- mulation [0] and approximate hypergradient descent [99 40) The former aims to reformulate the bilevel optimization problem into single-level
using the optimality conditions of the lower-level problem as constraints. However, this may impose large number of constraints for machine learning applications. The latter scheme directly solves the bilevel problem through iteratively updating the lower and upper-level parameters and, hence, is usually more efficient. Nevertheless, existing works
have mostly focused on unconstrained bilevel optimization [19 32
52 45 [4].
In this work, we study bilevel optimization problems where x and y are on Riemannian manifolds My and My. respectively. We focus on the setup where the lower-level function g 1/) is geodesic strongly convex (a generalized notion of convexity on manifolds, defined in Section in y. This ensures the lower-level problem has unique solution y (x) given The upper-level function can be nonconvex on M2 X My- Because the unconstrained bilevel optimization is special manifolds, such formulation includes wider class of applications. case of our formulation Examples of Riemannian bilevel optimization include Riemannian meta learning (64 and NAS over SPD networks (62). Moreover, there has been a surge of interest of min-max optimization over Riemannian manifolds (37 43 27 25 67 35 which also gets subsumed in the framework of bilevel optimization with g ff
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

Contributions. (i) We derive intrinsic Riemannian hypergradient via the implicit function theorem and propose four strategies for estimating the hypergradient, i.e., through Hessian inverse, conjugate gradient, truncated Neumann series, and automatic differentiation. We then provide hypergradient estimation error bounds for all the proposed strategies. (ii) We introduce the Riemannian hypergradi- ent descent algorithm to solve bilevel optimization problems on manifolds and provide convergence guarantees. We also generalize the framework to the stochastic setting and to allow the use of retraction. (iit) The efficacy of the proposed modeling is shown on several problem instances includ- ing hypper-reressentttinn over SPD matrices, Riemannian meta learning, and unsupervised domain adaptation. The proofs, extensions, and experimental details are deferred to the appendix sections.
Related works in unconstrained setting. Unconstrained bilevel optimization where the lower- level problem is strongly convex has been widely studied [19 [32] 40
[T 52
45 14. A crucial
ingredient is the notion of hypergradient bilevel optimization problems and its computation. There exist strategies for approximating the hypergradient, e.g., using conjugate gradient [4), Neumann series [1.. iterative differentiation [1], and Nyström method [31. While bilevel optimization with
constraints is relatively unexplored, few works exists that impose constraints only for the upper level problem [2 10] Recently, linearly lower-level constrained bilevel optimization has been explored
in 665 68
where a projected gradient method is employed for the lower-level problem.
Related works on manifolds. There has been limited work on bilevel optimization problems on manifolds. [7] studies semivectorial bilevel optimization on Riemannian manifolds where the upper- level is a scalar optimization problem while the lower-level a multiobjectivee problem under greatest coalition. (50 44] reformulate bilevel problems on manifolds into a single-level problem based on the KKT conditions on manifolds. However, for all those works, is unclear whether there exists an algorithm that efficiently solves the problem in large-scale settings. In contrast, we aim to provide a general framework for solving bilevel optimization on Riemannian manifolds. [47 is a contemporary work that also proposes gradient-based algorithms for bilevel optimization on Riemannian manifolds. The main differences of our work with respect to 47 are as follows: (1) We provide an analysis for various hypergradient estimators while [47] focuses on conjugate gradient for deterministic setting and Neumann series for stochastic setting; (2) We provide an analysis for retraction which is more computationally efficient than exponential map and parallel transport employed in [47], and (3) We explore the utility of Riemannian bilevel optimization in various machine learning applications, which is not the case with 47).
2 Preliminaries and notations
A Riemannian manifold M is a smooth manifold equipped with a smooth inner product structure (a Riemannian metric) (p T.M x T.M- + IR for any = E M and its tangent space T.M. The induced norm is thus ||u]|= = (u, u)æ for any u E T.MM A geodesic e: [0, 1] + M generalizes the line segment in the Euclidean space as the locally shortest path on manifolds. The exponential map on a manifold is defined as Exp.(u) = c(1) for a geodesic that satisfies c(0) = z, 2, c0)) = u. In totally normal neighbourhood 2 where exponential map has smooth inverse, the Riemannian distance d(xx, v) = 1Exx_"(u) (v) I|= [|Exp,, (ax) The parallel transport operation r M - T.,M is linear map which preserves the inner product, i.e., (u, u)2, = <T r) Vu, E M. The (Cartesian) product of Riemannian manifolds M2 x M, is also Riemannian manifold.
For a differentiable function ʃ .M - IR, the Riemannian gradient G5(=) E T.M is the tangent vector that satisfies (Gf(=), u): [=) for all u E T.M. Here D is the differential operator and [)( represents the directional derivative of f at z along u. For a twice differentiable function S, Riemannian Hessian fz.(2) defined as the covariant derivative of Riemannian gradient.
Geodesic convexity extends the convexity notion in the Euclidean space to Riemannian manifolds. A geodesic convex set Z C .Min is where any two points can be joined by a geodesic. A function f .M- + R is said to be geodesic (strongly) convex if for all geodesics e : [o. 1] - Z, f(((t) is (strongly) convex in f E [0, 1]. If the function is smooth, then f is called p-geodesie strongly convex if and only if f (Exp, (tu)) 2 fz) + t(ff(z), u)= + t2 [|ul]]. vt E [o, 1]. An equivalent second-order characterization is H(2)) Z pid, where we denote id as the identity operator.
For bifunction o: M2 x Mg + IR, we denote G. y): Gyb v) as the Riemannian (partial) gradient and Hz v). Hgo(z, (z. y) the Riemannian Hessian. The Riemannian cross-derivatives

---

are linear operators (z TyMy - T2M2. Gyd((z,)) T,M2 + TyMy defined as Pzy y)[u] D u for any v My (with D representing the differential operator) and similarly for (xr y). For linear operator T T,M2 + TyMy, the adjoint operator, denoted as T is defined with respect to the Riemannian metric, i.e. (T(u), = (2' u).- for any E T2Ma, My. The operator norm of is defined as |['lly = MuPueT,,M.:lvll..- [|
3 Proposed Riemannian hypergradient algorithm
In this work, we consider the constrained bilevel optimization problem min aEM, F()) f(z, y" (x)).
s.t. yx)) arg min g(r, y),
(1)
yeM,
where My are two Riemannian manifolds and J.g MM a x My + R are real-valued jointly smooth functions. We focus on the setting where the lower-level function g(r, y) geodesic strongly convex. This ensures the lower-level problem has unique solution y" (a) for given The upper-level function f can be nonconvex on MM x My.
We propose to minimize F(x) directly within the Riemannian optimization framework. To this end, we need the notion of the Riemannian gradient of F(x) f(x (x:), which we call the Riemannian hypergradient.
Proposition The differential of j (x) and the Riemannian hypergradient of F(x) are given by
Dy (x) = -H, g(x, y (z)) Gyzgzz (z, y² (x)) -H,,'g(a,,
(2)
(F(x) = G..f(a: y" (x)) Pzyg(z, y² (2)) g(z, (x))), y²,,)
The above proposition crucially relies on the implicit function theorem on manifolds and requires the invertibility of the Hessian of the lower level function f with respect to y. This is guaranteed in our setup as f geodesic strongly convex in y. Hence, there exists unique differentiable function y" (x: that maps æ to the lower-level solution. We show the Riemannian hypergradient descent (RHGD) algorithm for (T in Algorithm
Steps 3 to solve the lower-level problem using the Riemannian gradi- ent descent method. Since computing the optimal solution y (x:) is compu- tationally challenging, we obtain an approximate solution yk+1-
Algorithm 1 Riemannian hypergradient descent (RHGD)
1: Initialize E MM, E My
2: for O.... K do
3: v yk.
4: for O...., S- do
5:
Vi ++1 Expy, (- Gyg(Es,
Step involves computing the Riemannian hypergradient F(x)) of F(z). For computational effi-
6:
end for
7: Set yk++1 = yr.
ciency, compute approximation
8: Compute approximated hypergradient (F))))
F(()).
9: Update xx+1 Expay. (Fr)))
Step 9 is the usual exponential map find the updated point rx+1-
10:
end for
We highlight that Step 8 of Algorithm I |approximates the Riemannian hypergradient. In the rest of the section, we discuss various computationally efficient ways to estimate the Riemannian hypergra- dient and discuss the corresponding theoretical guarantees for RHGD. The error of hypergradient approximation comes from the inaccuracies of yk+1 to (xx) and also from the Hessian inverse.
3.1 Hypergradient estimation
When the inverse Hessian of the lower-level problem can be computed efficiently, estimate the hypergradient directly by evaluating the Hessian inverse (HINV) at yk++1 i.c., Fniny F'(x)) MI+1))9,ff(a (xF,, [k+1)]] However, computing the inverse Hessian is computationally expensive in many scenarios. We now discuss three practical strategies for estimating the Riemannian hypergradient when y/++1 is given.
fzf, ( G2 Ik, 3k++1) (H,, ( '9,
Conjugate gradient approach (CG). When evaluating the Hessian inverse difficult, we can solve the linear system Hy [u]] f 2 for some E Ty To this
3

---

end, we employ the tangent space conjugate gradient algorithm (Appendix F Algorithm 3 that solves the linear system on the langent space Ty My with only access to Hessian-vector products, i.e., (an) = Gz yh+1) 1)- Gayg yk+ [) where oig is computed as solution to Hyg(as, yh+1 [=f]] gy 1k+1), where T is the number of iterations of the tangent space conjugate gradient algorithm.
Truncated Neumann series approach (NC). The Neumann series states for an invertible operator H (id - HY', where id is the identity operator. An alterna- tive approach estimate the Hessian inverse is to use truncated Neumann series, which leads to the following approximated hypergradient, FusF(arx)
such that ||]]| || S I. its inverse H- = (z
- (Th, 1h++1)l Ei (id- y Haggaaa Me++ (Q (ars, 1)]l where 7 is chosen such that (id Hyg((ss Mk++)) >
G2 f (z 3)-11 g
can be set as Y Z. where the gradient operator is I.-L.ischitzz (discussed later in Definition Empirically, we observe that this approach is faster than the conjugate gradient approach. However, it
requires estimating 7' and L beforehand.
Automatic differentiation approach (AD). Another hypergradient estimation strategy follows the idea of iterative differentiation by backpropagation. After running several iterations of gradient update to obtain Yk++1 (which is function of zy), we can use automatic differentiation to compute directly the Riemannian gradient of f 4/++1 (x%)) with respect to Iik- We can compute the Riemannian hypergradient from the differential in the direction of arbitrary E My using basic chain rules.
3.2 Theoretical analysis
This section provides theoretical analysis for the proposed hypergradient estimators as well as the Riemannian hypergradient descent. First, we require the notion of Lipschitzness of functions and operators defined on Riemannian manifolds. Below, we introduce the definition in terms of bi-functions and bi-operators and state the assumptions that are required for the analysis.
Definition 1 (Lipschitzness). (1) For bifunction f M2 x My - R. say f has L Lipschitz satisfies for any x, 21, x2 E Uz y. y1, U2 E Gyf(2, y2)l S Ld(y1, v2), fl 31) - G.f(2, 12)llx S Ld(y1,92). 2, [ - G.f v2 S LL((11 12) and (1,,f(1,,)) ( w)ll, L2)
Riemannian gradient in Ug x Uy S Me x My if
Uy || (,f v1)
(2) For an operator 9(r, y) : Ty + TaMa, we say G(r, is p-Lipschitz if it satisfies, - ((22, (x),))
22 x2) and [19(x, [(r,, ) llz S pd((),92))
(3) For an operator H(z, v) My + My. we say H(r y) is p-Lipschitz if it satisfies, 81 (x,,))))))))))) y2 - H(x, 12) Ilva S pdly) y2) and [H((..,)) H(r2 1)ll, S 2p 22)
is worth mentioning that Definition D implies the joint Lipschitzness over the product manifold M2 x My, which is verified in Appendix C.2 Due to the possible nonconvexity for the upper level problem, the optimality is measured in terms of the Riemannian gradient norm of F(x).
Definition 2 (c-stationary point). We call MM c-stationary point of bilevel optimization 1 if it satisfies [())( <.
Assumption 1. All the iterates in the lower level problem are bounded in a compact subset that contains the optimal solution, i.e., there exists a constants D2 >> 0, for all k such that d(y y (#x)) s Dj for all s. Such neighbourhood has unique geodesic. We take D = max ( (Di,.., D
Assumption 2. Function f(z, y) has bounded Riemannian gradients. i.c., f y)ll, I M. [x..f(r.) y) & M for all (æ: y) E 2/ and the Riemannian gradients are I.-Lipschitz in U.
Assumption Function g(z, y) is u-geodesic strongly convex in y E U, for any x Udg and has L Lipschitz Riemannian gradient Gag(ı, y), Gyg(1, v) in U. Further, the Riemannian Hessian Hyg(2, v), cross derivatives G J. v), Gyn9(1, y) are p-Lipschitz in U.
Assumption D is standard in Riemannian optimization literature by properly bounding the domain of variables, which allows to express Riemannian distance in terms of (inverse) Exponential map. Also, the boundedness of the domain implies the bound on curvature, as is required for analyzing convergence for geodesic strongly convex lower-level problems 41 71 Assumptions 2 and are common regularity conditions imposed on f and y in the bilevel optimization literature. This translates into the smoothness of the function F and Dy" (x) (discussed in Appendix C.3).

---

Table 1: Comparison of first-order and second-order complexities for reaching c-stationarity. For stochastic algorithms, including HGD-NS, RSHGD-HINY, the complexities are measured with respect to the component functions 9i. 9s. Here, G5, Gg are the gradient complexities of function f. g. respectively, to reach an e-stationary point of (T). Also, we denote JV HV, as the complexity of computing the second-order cross derivative and and Hessian-vector product of function g.
Methods
G5
Gg
JVg
HV
HGD-CG
40
O(KGEt)
"(
O(
Ö(K) O(fic())
-AD 40
((ffc"))
SHGD-NS
40
0
2
RHGD-HINV
O(K(c))
Ö(ffce!)
O(x(c))
NA
-CG
O(kfe=t)
Ö(1)5
-NS
O(K7E-))
Ö(NicE-t)
O(KiEt)) O(K)ge"))
-AD
O
RSHGI--HINV
O
NA
We first bound the estimation error of the proposed schemes of approximated hypergradient as follows. For the hypergradient computed by automatic differentiation, we highlight that due to the presence of exponential map in the chain of differentiation, it is non-trivial to explicitly express D=, v Here, we adopt the property of exponential map (which is locally linear) in the ambient space [1. i.e., Exp,, (w) + u + O( This requires the use of tangent space projection of in the ambient space as Pa (E), which is solved for the e such that (u, C)a = (u, {) for any & E T.M.
For notation simplicity, we denote ią f. and Kip == f.- For analysis, we consider Kpp = (mi).
Lemma 1 (Hypergradient approximation error bound). Under Assumptions 2 error for approximated hypergradient as
2
3
we can bound the
1. HINV: Iinn F(x)) GF S S( M + K2 + KIR Ad (xx). 1k+1).
2. CG: 11 ((r)) - (F(x))) S (L. M + L( + 2//1i) (Kt + Maa ) d(y (mx)), UK+1)
VKi-1 T
2L 2L//
110g peets ((x) villbx where of = H ( 1/" (Ix)) (Gn))) (fx, y (xx)))-
Vet+1
3. NS: (i - GF( I U MLL + M + KIR M M)d(y (Zx). vh+1) + Ki Kg (1 711
4. AD: Suppose further there exist C1, C2 C3 > such that [D live I C₁ | (E) (id + DRu) + E where || Bxpp S
olly S C2d(a, v)llll and Ex.
Cal Cullllullel/ for any x, U and v E T Mg. p, E T,Me. Then,
F GF(#x)llz S ( 1 L(1 + ki)) (1 + n nyH) 8-1 d(yk. y (œx)) +
2MC
MR - ny111)444 nyH4s, where Ō = (m + 1)p+ (C2 + ny ((1-ng + tn
From Lemma it is evident that the exact Hessian inverse exhibits the tightest bound, which is followed by conjugate gradient (CG) and truncated Neumann series (NS). Automatic differentiation (AD) presents the worst upper bound on the error due the introduction of curvature constant
resulting in (1 - 0 tx ))$ = -1 O )) for the trailing term, which could be much larger
than (1 - 7/4)" = (1 (T))) for NS and ( (1 - O(
( 1 )) for CG. Further, the Aay
error critically relies on the number of inner iterations S compared with T for CG and NS, and
the constants C1, C2, C3 can be large for manifolds with high curvature. We now present the main convergence result with the four proposed hypergradient estimation strategies.
Theorem Denote Ao F(ro) d y (zo)) and LF (
PLM
PLM 1
(ke)), Under Assumptions
2
3 we have the following bounds on the hypergradient norm obtained
by Algorithm 7
HINV: Let qr 201
and 2 @(x{(). We have mink=o..... K-1 [1FF(x=))2
80L,A0/ K.
5

---

CG: Let A C x2{ (2( 225M000 + 1). where Cv Mig
Mi. + xj + Kj. Choosing
24 2 (({{(), and Tog 2 @( we have o K-1 || (*x)()x)1))))))))))))())))))(( S
1
Th
96A K
(Au
+ 11 (
NS: Choosing f/5 20LF , (nf((), and Tus 2 ((xlog())) for an arbitrary
0, we have
mink=o...... K-1 [[ I S
+ 2. and z
D
k
AD: Choosing aJ = 20LF
Onf (og()))()) for an arbitrary 6 > 0, we have
mitk OK..,K -- F (1x)) S SOLE do
Complexity analysis. Based on the convergence guarantees in Theorem I we have analyzed (in Corollary D the computational complexity of the proposed algorithm with four different hypergradi- ent estimation strategies in reaching the e-stationary point. The results are summarized in Table] T For reference, we also provide the computational cost of Euclidean algorithms which solve bilevel Euclidean optimization problem [40]. We notice that except for CG, the gradient complexity for f (i.e., Gy) matches the Euclidean version. For conjugate gradient, the complexity is higher by (1s). which is due to the additional distortion from the use of vector transport when tracking the error of conjugate gradient each epoch. In terms of gradient complexity for g (i.e.. Gg). all deterministic methods require higher complexity by at least ((r(() compared to the Euclidean baselines. This is because of the curvature distortion when analyzing the convergence for geodesic strongly convex functions. Similar comparisons can be also made with respect to the computations of cross-derivatives and Hessian vector products.
3.3 Extension stochastic bilevel optimization
In this section, we consider the bilevel opti- mization problem [ in the stochastic setting, where fx (xr) = + fi (z)) and g(x, == 12
Algorithm 2 Riemannian stochastic bilevel optimization with Hessian inverse.
1: Initialize zo E MM y/ E My.
Σr (x, The algorithm for solving the stochastic bilevel optimization problem is in Algorithm [ where sample B1, B2, B3, B3 afresh every iteration. The batch index is omitted for clarity. The batches are sampled uniformly at random with replacement such that the mini-batch gradient an unbiased estimate of the full gradient. Here, we denote fø (zz y)
2: for k O.... K 11 do
3: v = Vks
4: for 0, ...,. S. 1 do
5:
Sample a batch Bi.
6:
of Expry (- 9ugs, (m,, vf)).
7:
end for
8:
Set vk++ yh-
9:
Sample batches B2, B3. B1.
[6] EiEs f( v) and similarly for g. We let [n] = nI....,
10: Compute (xxx)
11: Update rk+1 Expaa (-T)))(Ix))
In Step 10 of Algorithm 2 we can employ any hypergradient estimator proposed in Section In this work, only show convergence der the Hessian inverse approximation of hypergradient, i.c., (F((x))
12:
end for
Je (I&,, yk++1)
Pgy.9. (TK, 3/++1) [H 9B, (Kk, Vk++) I (I 1)]l Similar analysis can be followed for other approximation strategies. The theoretical guarantees are in Theorem 2 where we require Assumption
4 which is common in existing works for analyzing stochastic algorithms on Riemannian manifolds
42 23 22
Assumption 4. Under stochastic setting, Assumption I holdd and Assumptions 2 3 are satisfied for component functions (z, y), gy (z, y), for all E [n].] E [m]] Further, stochastic gradient, Hessian, and cross derivatives are unbiased estimates.
Theorem 2. Under Assumption 4 consider Algorithm Suppose we choose 1/s- :
((+{s), and |B,1 [B2), [Bs). |Ba mink=o.... K-1 ElggF 112 < S0L + and the gradient complexity to reach e-stationary
O(Kie-1) for arbitrary > 0. Then have
solution is Gy = On7e -2, Gg = Ö(rice IV, = O
In Table T we compare our attained complexities with that of stocBiO 40, which makes use of truncated Neumann series. With exact Hessian inverse, we can match the Gy and J complexities

---

with stocBio. For the Gg complexity, the additional curvature constant is inevitable from the convergence analysis for geodesic strongly convex functions. Nevertheless we observe the same order dependency on Kt. This is mainly due to the analysis where we choose smaller stepsize fy = = (z compared θ (( L+H
40). The larger stepsize, despite increasing the convergence
rate, also increases the variance under stochastic setting. We believe an order of (Kz) lower can be established for stocBio, following our analysis.
3.4 Extension to retraction
Our analysis till now has been limited to the use the exponential map. However, the retraction mapping is often preferred over the exponential map due to its lower computational cost. Here, we show that use of retraction our algorithms also leads to similar convergence guarantees.
Assumption 5. There exist constants ē 2 1, 2 0 such that (xx,)) y) S a cllll and |[Exp (w)
u|I < cn ||, for any 2, 1 Retr, (1) E E..
in bounding the error between exponential map and retraction given that retraction is a first-order approximation to the exponential map.
Assumption | is standard (e.g. in 42 23
Theorem Suppose Assumptions / and Soll and let LF 4xTeRM 50Lr. Then consider with exponential map replaced with general retraction. We can obtain the following bounds.
Algorithm 7
HINV: Let 72 = 0(1//p)) ((((()) Then mink ([9F I 16L
C2 CG: (1/A),, @(*i(), Teg ((VKT, where λ = n² "
+ c). Then
mink=o,.., K- (x) s _ (Do ||:511 (mo))
NS: Let 0(1/Lp), 2 (k]{(). Then for an arbitrary e > 0, Tus 2 (mg we 162r
have ming=o,... K- [F <
<1x= Du tots 2.
AD: Let The = 0(1/Lp), s 2 0(G) (og(()))))))).. Then for an arbitrary > 0, we have m K K-1 [ I 16L Ao +
Theorem 3 demonstrates that employing a general retraction preserves the same order of convergence and complexity as the exponential map in Theorem 1 This is due to the fact that L ((Lr) and Ā (A), where Lp and A are defined in Theorem In addition, when exponential map is used, Theorem recovers the results in Theorem as CR = Oand č= 1.
Experiments
This section explores various applications of bilevel optimization problems over manifolds. All the experiments are implemented based on Geoopt 44 and the codes are available at hhtpss //githuu.ccom//ne wtt.uub...mm//ayyn////ar/.ee
4.1 Synthetic problem
We consider the following bilevel optimization problem on the Stiefel manifold (t,, r) (W Ødxr W W I,} and SPD manifold s (M Rdddd M > 0} (in Appendix A
tr(M-XTTWW) YWT), M arg min (M,X X) + (M-1, WY YWT vI).
West(d.r)
ME8y
where E Raxd Y E Roxr, with 2 d2 r, are given matrices and parameter. The above is synthetically constructed problem that aims to maximize the similarity between X and Y in different feature dimensions. We align X and to the same dimension via W E rt(d,r) and also learn an appropriate geometric metric M E $ in the lower-level problem [69]. The geodesic convexity of the lower-level problem and the Hessian inverse expression are discussed in Appendix H.I
the regularization
7

---

2010
1350
0930
Pypeerree WITET
2010
0051
0:20
m250
20122
Epocts
Spocta
Spocts
(a)
(b)
(c)
(d)
Figure 1: Figures (a) & (b) show the plot of objective of the upper-level problem (Upper Objective) for different strategies. HINV and CG strategies have fastest convergence, followed by NS and AD. The corresponding estimation errors are shown in (c). Figure (d) specifically shows the robustness of approximation error obtained by NS across different and T values.
Results. We generate random data matrices X, Y with 100, d 50, and 20. We set = 0.01 and fix f = 0.5. We compare the three proposed strategies for approximating the hypergradient where we select y 1.0 and Tus 50 for Neumann series (NS) and maximum iterations Tag for conjugate gradient (CG) to be 50 and break once the residual reaches a tolerance of 10-10. We set the number of outer iterations (epochs) K to be 200. Figure I compares RHGD with different approximation strategies implemented with S 20 50 number of inner iterations.
4.2 Hyper-representation over SPD manifolds
Hyper-reppresentation [54] 61 aims to solve a regressionnclaassffcatiio task while searching for the
best representation of the data. It can be formulated as a bilevel optimization problem. where the lower-level optimizes the regression/Classffccttin parameters while the upper-level searches for the optimal embedding of the inputs. Suppose we are given a set SPD matrices, D {A,} where A; 8q and the task is learn low-dimensional embedding of A while remaining close their semantics labels. In particular, we partition the set into training set Dr, and validation set Deal-
Shallow typer-repressnttaioon for regression. We consider shallow learning paradigm D through the regression task. The representation is parameterized with W A, for W e rt. The requirement of orthogonality on W follows [8 representations are SPD. The learned representation is then transformed to a Euclidean space for performing regression, namely through matrix logarithm (that acts bijective map between the space of SPD matrices and symmetric matrices) and a vectorization operation vec(.) that extract the upper-triangular part of the symmetric matrix. The bilevel optimization problem is
30 33
that ensures the learned
(vec(logmt W A. w)) -0:
min WeSt(d,)) IeDi
2P
(vec(logm) T B + 1113.
s.t. 3 arg
2P
BERrirtE/2 IED
The regularization À > 0 ensures the lower-level problem is strongly convex. The upper-level problem on the validation set while the lower-level problem is on the training set. We generate random W, A, and ( and construct y with y; vec(logm W A, W) + e;, where e; (N 1). We generate 200 A; with [Dvall = 100 and [D.r = 100. In Figure 2a we show the loss on validation set (the upper loss) in terms of number of outer iterations. We compare both the deterministic (RHGD) and stochastic (RSHGD) versions of Riemannian hypergradient descent. We again observe that the best performance is attained by either the ground-truth Hessian inverse the conjugate gradient. NS requires carefully selecting the hyperparameters 7, T, which pose difficulties real applications. For the stochastic versions, all the methods perform similarly.
Deep hyper-representation for classification. We now explore 2-layer SPD network (38] for classifying ETH-80 image set [46). The dataset consists of 8 classes, each with 10 objects. Each object is represented by an image set consisting of images taken from different viewing angles. Here, represent each image resizing them into 10 x 10. This results in 80 SPD matrices A, of size 100 x 100 for classification. Let ((A,) vec(logm W, ReEig( WT W.A. A, WT W2)) be the output of the 2 layer network where ReEig( A) U max(el. el, E)UT is the eigenvalue rectifying activation with the eigenvectors U and
by taking the covariance matrix of the images in the same set after

---

Upper look
Vil Aacuuueee
15 15-e0-ica
BH00-3
AMG2--a
Spocti
Epochs
Time
Buntims
(a) Shallow hyperrep
(b) Deep hyperrep
(c) Deep hyperrep
(d) MinilmageNett
Figure 2: Figures (a), (b), and (c) show the performance of RHGD the hyper-representation problems on SPD networks. Figure (d) shows the good generalization performance our proposed RHGD algorithms over the projected gradient PHGD baselines on the MinilmageNet dataset.
Table 2: Classification accuracy on the Caltech-Office dataset.
Methods
A-C
ADD
A-W
C-AA
C->D
C-
D-
D-C
D-W
W >
W->
W-+D
OT-EMD
66.67
47.77
45.76
67.52
36.31
42.71
62.17
59.71
85.08
55.41
51.16
96.82
OT-SKH
76.83
75.80
69.83
84.35
78.34
68.14
80.92
71.57
93.90
74.17
67.02
87.26
Proposed
78.70
80.25
69.83
88.21
80.25
68.47
82.74
75.69
97.97
83.49
73.62
98.73
eigenvalues E of A. We consider the same bilevel optimization as above except the least-squares loss function becomes the cross-entropy loss. Here we sample 5 samples from each class form the training set and the rest as the validation set. We set dj 20, d2 5, and fix learning rate to be 0.1 for both lower and upper problems. Figures 2b and 2c show the good performance on the validation accuracy (upper-level loss).
4.3 Riemannian meta learning
allows adaptation of models to new tasks with minimal amount of additional data and training, by distilling past learning experiences. A recent work [64] considers meta learning with orthogonality constraint. In particular, the upper-level optimization searches for the base parameters shared by all tasks while the lower level optimizes over the task-specific parameters to ensure generalization ability. Let Py denote the distribution of meta tasks and for each training epoch, we sample m tasks De Pr, 1, ..., m. Each task is composed of support and query set denoted by Dr. Dr and the task is to learn set of base parameters Q such that the model can quickly adapt to the query set from the support set by adjusting only few parameters w. For each task, the task-specific parameter wi is learned from the support set, which is used to update the base parameters by minimizing the loss over the query set. standard settings, we corresponds to the final linear layer of neural network (39 40). Here, we adopt the setup with wy to be the last layer parameters in the Euclidean space while enforcing O on the Stiefel manifold. The problem of Riemannian meta- learning is mineest 1. >>> Z(O, wf: ; DE) s.t. w arg min 2 L(O, We; D + R(wx), where D., Dí are the support and query sets for task and R(-) is regularizer that ensures strong convexity of the lower-level problem.
Meta learning 16 34
55 where the backhone network is a 4-block CNN with the kernel of the first 2 layers constrained to be orthogonal in terms of the output channel (following (48). The kernel size is 3 x and consider 16 output channels with a padding of 1. Each convolutional block consists of a convolutional layer, followed by ReLU activation, a max-pooling and batch normalization layer. O, thus, has the dimension (16 * 3 * 3) x 16 = 144 x 16, which is constrained to the Stiefel manifold.
Results. We consider 5-ways 5-shots meta learning over the MinilmageNet dataset
In Figure 2d plot the test accuracy averaged for over 200 tasks. We compare RHGD with an extrinsic update baseline PHGD, which projects the update from the Euclidean space to the Stiefel manifold at every iteration. We observe the RHGD converges faster compared to the extrinsic update PHGD, thereby showing the benefit of the Riemannian modeling.
4.4 Unsupervised domain adaptation
Given two marginals I E R v E Rm with equal total mass, i.e., p 1I, TTI, 1m where assume unit mass without loss of generality. Let 1I {I Raxm r 0,

---

μ. 1,, = be the set of doubly stochastic matrices with strictly positive entries. From [5 57 is known that the set forms Riemannian manifold with the Fisher metric.
Given a supervised source dataset X E Raxd and an unsupervised target dataset Y E R (n,m 2 d), consider the unsupervised domain adaptation problem to classify target domain instances. Using the optimal transport framework [58] 12), we pose this as a bilevel problem:
minreN((..)) (I, CD.C((MM))) 1/2 VM--1/1
(T)),
(3)
s.t. M' arg min min adist (M, X x + (1 oidist PM, Y TYY,,
where HIT) - (T. Iog is the entropy regularization, C(X,Y) 2XY is the pairwise squared distance matrix, and IT denotes the doubly stochastic manifold 15 Here, dist. is the geodesic distance between SPD matrices and a E [0.
diag(XX (oex(XX)1,, 1
1,diag( YY
The lewer-level problem in (3) finds M which is the weighted geometric mean between X ' and TTTTY [6]. Conceptually, learning of M allows to align the features of the source and target instances. The upper-level problem, on the other hand, minimizes the Mahalanobis distance between the source and target domain instances parameterized by M An interpretation is that the matrix M' leads to whitening of the data (i.e., -M2 the whitened data) S5]
After the transport plan I' is learned, we employ barycentric projection using I to transport the source points to the target domain and and employ the nearest neighbour (I-NN) classifier for the target dataset classification. For barycentric projection, we project the source samples X to the target
Y by solving x; = arg eR Σj= T | | |/ - M*-1/2y, y,12
Pi EEj= Fi))))) Then, a hearest-neighhout (NN) classifier is used to classify the samples in the target given the source labels based on the distance computed with M i.e., C(XM-11/2 -M/2)
Results. We consider the Caltech- dataset [20), which is commonly used for domain adaptation. The dataset contains images from four domains in ten classes, i.e., Amazon (A), the Caltech image dataset (C), DSLR (D), and Webcam (W), each with containing 958, 1123, 157, and 295 samples respectively. Hence, there are 12 domain adaptation tasks, e.g., A-D implies A is the source and D the target. Each domain has the same ten classes. The goal is to classify images from target domain given source domain. For preprocessing, we normalize the samples to have unit norm and reduce the dimensionality 128 by mean pooling every 64 columns.
We compare our proposed bilevel approach (3) with single-level optimal transport baselines, i.e., solv- ing (l))))( (T, C(X, Y)) - (T), followed by the same barycentric projection. Specifically, the baselines are: (1) optimal transport where X = 0 (labelled as OT-EMD) and (2) optimal transport with the Sinkhorn algorithm (labelled as OT-SKH), OT-EMD employs the earth mover distance while OT-SKH employs the Sinkhorn distance [33. We implement the two OT baselines with the POT Python library [17] À is tuned in OT-SKH for each source-target pair. The best validation results are obtained by setting À= = 5 x 10 for all the problem pairs except the W-D pair for which A 10 gives the best result. For our proposed bilevel approach, we set À = 0 and a 0.5.
In Table we observe that the proposed bilevel approach obtains better generalization performance than the baselines across all the tasks. This showcases the utility of learning the whitening metric M-1 in a bilevel setting.
5 Conclusion
In this work, we have proposed framework for tackling bilevel optimization over Riemannian mani- folds. We discussed various hypergradient approximation strategies (conjugate gradients, truncated Neumann series, and automatic differentiation) and provide error bounds. Our proposed algorithms rely only on gradient updates and make use of retraction which scale well across problems. We illustrate the efficacy of the proposal approach in several machine learning applications.
Although in this work, we focus on geodesic strongly convex lower-level problems, our framework can be extended to relax such assumption to (geodesic) convexity with an extra strongly convex regularizer 22 or to (Riemannian) PL condition where global minimizer exists Ø believe there is potential to improve the current results in stochastic bilevel optimization by reducing the strict requirements on batch size. Additionally, the dependency on the curvature constant could also be further optimized.
Furthermore, we
10

---

References
[1] P-A Absil, Robert Mahony, and Rodolphe Sepulchre. Optimization algorithms on matrix manifolds. Princeton University Press, 2008
[2] Jan Harold Alcantara and Akiko Takeda. Theoretical smoothing frameworks for general nonsmooth bilevel problems. Technical report, arXiv preprint arXiv:2401.17852. 2024.
[3] Foivos Alimisis, Antonio Orvieto, Gary Bécigneul, and Aurelien Lucchi. A continuous time perspective for modeling acceleration in Riemannian optimization. International Conference on Artificial Intelligence and Statistics (AISTATS), 2020.
[4] Vincent Arsigny, Pierre Fillard, Xavier Pennec, and Nicholas Ayache. Geometric means in novel vector space structure on symmetric positive-definite matrices. SIAM Journal on Matrix Analysis and Applications, 29(1):328 347, 2007.
[5] Mikel Artetxe, Gorka Labaka, and Eneko Agirre. Generalizing and improving bilingual word embedding mappings with multi- framework of linear transformations. In the AAAI Conference on Artificial Intelligence (AAAI). 2018.
[6] Rajendra Bhatia. Positive definite matrices. In Positive Definite Matrices. Princeton university press, 2009.
[7] Henri Bonnel, Léonard Todjihoundé, and Constantin Udrişte. Semivectorial bilevel optimization on Riemannian manifolds. Journal of Optimization Theory and Applications, 167:464 2015.
[8] Nicolas Boumal. An introduction to optimization on smooth manifolds. Cambridge University Press, 2023.
[9] Lesi Chen, Jing Xu, and Jingzhao Zhang. On finding small hyper-gradients in bilevel opti- mization: Hardness results and improved analysis. In Annual Conference on Learning Theory (COLT), 2024.
[10] Tianyi Chen, Yuejiao Sun, Quan Xiao, and Wotao Yin. A single-timescale method for stochastic bilevel optimization. In International Conference on Artificial Intelligence and Statistics (AISTATS), 2022.
[11] Tianyi Chen, Yuejiao Sun, and Wotao Yin. Closing the gap: Tighter analysis of alternating stochastic gradient methods for bilevel problems. In Advances in Neural Information Processing Systems (NeurIPS), 2021.
[12] Nicolas Courty, Rémi Flamary, Devis Tuia, and Alain Rakotomamoniy. Optimal transport for domain adaptation. IEEE Transactions on Pattern Analysis and Machine Intelligence. 39(9):1853-1888 2017.
[13] Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In Advances in Neural Information Processing Systems (NeurIPS), 2013.
[14] Mathieu Dagréou, Pierre Ablin, Samuel Vaiter, and Thomas Moreau. framework for bilevel optimization that enables stochastic and global variance reduction algorithms. In Advances in Neural Information Processing Systems (NeurIPS), 2022.
[15] Ahmed Douik and Babak Hassibi. Manifold optimization over the set of doubly stochastic matrices: A second-order geometry. IEEE Transactions on Signal Processing, 67(2)):5761- 5774, 2019.
[16] Chelsea Finn, Pieter Abbecl, and Sergey Levine. Model-agnostic meta-learning for fast adapta- tion of deep networks. In International Conference on Machine Learning (ICML), 2017.
[17] Rémi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar Alaya, Aurélie Boisbunon, Stanislas Chambon, Laetitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, et al. POT: Python optimal transport. Journal of Machine Learning Research, 22(1):3571-3578
2021.

---

[18] Luca Franceschi, Paolo Frascomi, Saverio Salzo, Riccardo Grazzi. and Massimiliano Pontil. Bilevel programming for hyperparameter optimization and meta-learning. In International Conference on Machine Learning (ICMI.), 2018.
[19] Saeed Ghadimi and Mengdi Wang. Approximation methods for bilevel programming. Technical report, arXiv preprint arXiv:1802.022
2018.
[20] Boqing Gong, Yuan Shi, Fei Sha, and Kristen Grauman. Geodesic flow kernel for unsupervised domain adaptation. In Conference Computer Vision and Pattern Recognition (CVPR), 2012.
[21] Riccardo Grazzi, Luca Franceschi, Massimiliano Pontil, and Saverio Salzo. On the iteration complexity of hypergradient computation. In International Conference on Machine Learning (ICML), 2020.
[22] Andi Han and Junbin Gao. Riemannian stochastic recursive momentum method for non-convex optimization. In International Joint Conference on Artificial Intelligence (IICAI), 2021.
[23] Andi Han and Junbin Gao. Improved variance reduction methods for Riemannian non-convex optimization. IEEE Transactions on Pattern Analysis and Machine Intelligence, 44(11):7610- 7623, 2022.
[24] Andi Han, Bamdev Mishra, Pratik Jawanpuria, and Junbin Gao. Learning with symmetric posi- tive definite matrices via generalized Bures-Waat Wasserstein geometry. In International Conference on Geometric Science of Information (GSI), 2023.
[25] Andi Han, Bamdev Mishra, Pratik Jawanpuria, and Junbin Gao. Noneonvex-noneconan min- max optimization on Riemannian manifolds. Transactions on Machine Learning Research, 2023.
[26] Andi Han, Bamdev Mishra, Pratik Jawanpuria, and Junbin Gao. Riemannian accelerated gradient methods via extrapolation. In International Conference on Artificial Intelligence and Statistics (AISTATS), 2023.
[27] Andi Han, Bamdev Mishra, Pratik Jawanpuria, Pawan Kumar, and Junbin Gao. Riemannian Hamiltonian methods for min-max optimization on manifolds. SIAM Journal on Optimization, 33(3):1777-1877, 2023.
[28] Andi Han, Bamdev Mishra, Pratik Kumar Jawanpuria, and Junbin Gao. On Riemannian optimization over positive definite matrices with the Musse-Wassrrstein geometry. Advances in Neural Information Processing Systems (NeurIPS), 2021.
[29] Pierre Hansen, Brigitte Jaumard, and Gilles Savard. New branch-and-bound rules for linear bilevel programming. SIAM Journal on Scientific and Statistical Computing, 13(5):1194-12) 1992.
[30] Mehrtash Harandi, Mathieu Salzmann, and Richard Hartley. Dimensionality reduction on SPD manifolds: The emergence of geometry methods. JEEE Transactions Pattern Analysis and Machine Intelligence, 40(1):48 62, 2017.
[31] Ryuichiro Hataya and Makoto Yamada. Nyström method for accurate and scalable implicit differentiation. In International Conference on Artificial Intelligence and Statistics (AISTATS), 2023.
[32] Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang. A two-timescale stochastic algorithm framework for bilevel optimization: Complexity analysis and application to actor- critic. SIAM Journal on Optimization, 33(1):147-180, 2023.
[33] Inbal Horev, Florian Yger, and Masashi Sugiyama. Geometry- principal component analysis for symmetric positive definite matrices. In Asian Conference on Machine Learning (ACML), 2016.
[34] Timothy Hospedales, Antreas Antoniou, Paul Micaelli, and Amos Storkey. Meta-learning in neural networks: A survey. IEEE Transactions on Puttern Analysis and Machine Intelligence, 4449)15149-5169, 2021.
12

---

[35] Zihao Hu, Guanghui Wang, Xi Wang, Andre Wibisono, Jacob D Abernethy, and Molei Tao. Extragradient type methods for Riemannian variational inequality problems. In International Conference Artificial Intelligence and Statistics (AISTATS), 2024.
[36] Feihu Huang and Shangqian Gao. Riemannian gradient methods for stochastic composition problems. Neural Networks, 153:224-234, 2022.
[37] Feihu Huang and Shangqian Gao. Gradient descent ascent for minimax problems on Riemannian manifolds. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2023,
[38] Zhiwu Huang and Luc Van Gool. A Riemannian network for spd matrix learning. In AAAI Conference Artificial Intelligence (AAAI), 2017.
[39] Kaiyi Ji, Jason D Lee, Yingbin Liang, and H Vincent Poor. Convergence of meta-learning with task-specific adaptation over partial parameters. In Advances in Neural Information Processing Systems (NeurIPS), 2020.
[40] Kaiyi Ji, Junjie Yang, and Yingbin Liang, Bilevel optimization: Convergence analysis and enhanced design. In International Conference on Machine Learning (ICML), 2021.
[41] Michael Jordan, Tianyi Lin, and Emmanouil- Vaaiaeouil-Vasileioss Vlatakis- SaaaaiisGGaragkounii. First-order algo- rithms for min-max optimization in geodesic metric spaces. In Advances in Neural Information Processing Systems (NeurIPS), 2022.
[42] Hiroyuki Kasai, Hiroyuki Sato, and Bamdev Mishra. Riemannian stochastic recursive gradient algorithm. In International Conference on Machine Learning (ICML), 2018
[43] Philip A Knight. The Sinkhorn- Knopp algorithm: convergence and applications. SIAM Journal on Matrix Analysis and Applications, 30( (100))11:122))) 2008
[44] Max Kochurov, Rasul Karimov, and Serge Kozlukov. Geoopt: Riemannian optimization in Py PyTorh. Technical report, arXiv preprint arXiv:2005.02819, 2020.
[45] Jeongyeol Kwon, Dohyun Kwon, Stephen Wright, and Robert D Nowak. A fully first-order method for stochastic bilevel optimization. In International Conference on Machine Learning (ICML), 2023.
[46] Bastian Leibe and Bernt Schiele. Analyzing appearance and contour based methods for object categorization. In Conference on Computer Vision and Pattern Recognition (CVPR), 2003.
[47] Jiaxiang Li and Shiqian Ma. Riemannian bilevel optimization. Technical report, arXiv preprint arXiv:2402.020) 2024.
[48] Jun Li, Fuxin Li, and Sinisa Todorovic. Efficient Riemannian Optimization on the Stiefel Manifold via the Cayley Transform. In International Conference on Learning Representations (ICLR), 2019.
[49] Jiagen Liao and Zhongping Wan. Inexact restoration methods for semivectorial bilevel pro- gramming problem on Riemannian manifolds. Axioms, 11(12):696, 2022.
[50] Jiagen Liao and Zhongping Wan. On the Sarush-Kuhn-Tuucee reformulation of the bilevel optimization problems on Riemannian manifolds. Filomat, 36(1 1) 1):3609-3224, 2022.
[51] Zhenhua Lin. Riemannian geometry of symmetric positive definite matrices via Cholesky decomposition. SIAM Journal on Matrix Analysis and Applications, 40(4):1353-1370, 2019.
[52] Bo Liu, Mao Ye, Stephen Wright, Peter Stone, and Qiang Liu. Bome! bilevel optimization made easy: simple first-order approach. In Advances in Neural Information Processing Systems (NeurIPS), 2022.
[53] Hanxiao Li Liu, Karen Simonyan, and Yiming Yang. DARTS: Differentiable architecture search. In International Conference on Learning Representations (ICLR), 2018.
13

---

[54] Jonathan Lorraine, Paul Vicol, and David Duvenaud. Optimizing millions of hyperparameters by implicit differentiation. In International Conference on Artificial Intelligence and Statistics (AISTATS), 2020.
[55] Luigi Malagò, Luigi Montrucchio, and Giovanni Pistone. Wasserstein Riemannian geometry of Gaussian densities. Information Geometry, 1:137-179, 2018.
[56] David Portinez-Rubiio Christophe Roux, Christopher Criscitiello, and Sebastian Pokutta. Ac- celerated methods for Riemannian min- max optimization ensuring bounded geometric penalties. Technical report, arXiv preprint mtXiv:2305.16186 2023.
[57] Bamdev Mishra, NTV Satyadev, Hiroyuki Kasai, and Pratik Jawanpuria. Manifold optimization for non-linear optimal transport problems. Technical report, arXiv preprint wxXiv22003009002 2021.
[58] Gabriel Peyre and Marco Cuturi. Computational optimal transport: With applications to data science. Foundations and Trends in Machine Learning, 11(5-6)3355-66 607, 2019.
[59] Mengye Ren, Sachin Ravi, Eleni Triantafillou, Jake Snell, Kevin Swersky, Josh B. Tenen- baum, Hugo Larochelle, and Richard S. Zemel. Meta-learning for semi-supervised few-shot classification. In International Conference on Learning Representations (ICLR), 2018.
[60] Chenggen Shi, Jie Lu, and Guangquan Zhang. An extended Kuhn-Tucker approach for linear bilevel programming. Applied Mathematics and Computation, 162(1):51 63, 2005.
[61] Daouda Sow, Kaiyi Ji, and Yingbin Liang. On the convergence theory for Hessian-free bilevel algorithms. Advances Neural Information Processing Systems (NeurIPS), 2022.
[62] Rhea Sukthanker, Zhiwu Huang, Suryansh Kumar, Erik Endsjo Goron, Yan Wu, and Luc Van Gool. Neural architecture search of spd manifold networks. In International Joint Confer- ence on Artificial Intelligence (IJCAT), 2021.
[63] Yue Sun, Nicolas Flammarion, and Maryam Fazel. Escaping from saddle points on Riemannian manifolds. Advances in Neural Information Processing Systems (NeurIPS), 2019.
[64] Hadi Tabealhojeh, Peyman Adibi, Hossein Karshenas, Soumava Kumar Roy, and Mehrtash Ha- randi. RMAML: Riemannian meta-learning with orthogonality constraints. Pattern Recognition, 140: 109563, 2023.
[65] Ioannis Tsaknakis, Prashant Khanduri, and Mingyi Hong. An implicit gradient-type method for linearly constrained bilevel problems. In International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2022.
[66] Mengdi Wang, Ethan X Fang, and Han Liu. Stochastic compositional gradient descent: algo- rithms for minimizing compositions of expected-value functions. Mathematical Programming, 161:419-449, 2017.
[67] Xi Wang, Deming Yuan, Yiguang Hong, Zihao Hu, Lei Wang, and Guodong Shi. Riemannian optimistic algorithms. Technical report, arXiv preprint arXiv:2308. 16004, 2023.
[68] Quan Xiao, Han Shen, Wotao Yin, and Tianyi Chen. Alternting projected sgd for equality- constrained bilevel optimization. In International Conference on Artificial Intelligence and Statistics (AISTATS), 2023.
[69] Pourya Zadeh, Reshad Hosseini, and Suvrit Sra. Geometric mean metric learning. In Interna- tional Conference on Machine Learning (ICML), 2016.
[70] Dewei Zhang and Sam Davanloo Tajbakhsh. Riemannian stochastic gradient method for nested composition optimization. Technical report, arXiv preprint arXiv:2207.099: 0095500 2022.
[71] Hongyi Zhang, Sashank J Reddi, and Suvrit Sra. Riemannian SVRG: Fast stochastic opti- mization on Riemannian manifolds. In Advances in Neural Information Processing Systems (NeurIPS), 2016.
14

---

[72] Hongyi Zhang and Suvrit Sra. First-order methods for geodesically convex optimization. In Conference on Learning Theory (COLT), 2016.
[73] Peiyuan Zhang, Jingzhao Zhang, and Suvrit Sra. Sion's minimax theorem geodesic metric spaces and Riemannian extragradient algorithm. SIAM Journal Optimization, 33(4):2885- 2908, 2023.
15

---

Appendix
A Riemannian geometries of considered manifolds
Symmetric positive definite (SPD) manifold. The SPD manifold of size d denoted as $4 {X E Rdxd X = X.X > 0} and the commonly considered Riemannian metric is the affine- [6. for U, V Tx Other Riemannian metrics, invariant metric (U.V), tx
such (generalized) Watrs-Wassersteinn [55 24 Fuglidnn 4 and log-Cholesky (51] metrics
can also be considered. The exponential map is given by Expx (U) Xexpm [X U where expm(-) denotes the principal matrix exponential. The corresponding logarithm map is given by logx() (Y) = Xlogm X YY). Its Riemannian gradient of a real-valued function f derived as grad (X) = X and the Riemannian Hessian is Hessf (UX (X))s X (X)]UUX (UV](X)X Is where
(
:= (A A )/2. Stiefel manifold. The Stiefel manifold is the set of orthonormal matrices, i.e., (t, = (X
(A}s
pdxr, X X = 1]. A common Riemannian metric the Euclidean inner product. We consider the QR-based retraction in the experiment, which is Retrx (U) qf(X + U) where qf(-) extracts the Q-factor from the QR decomposition. Let the orthogonal projection the tangent space be denoted as (U) = U - X Uts. Then, the Riemannian gradient and Riemannian Hessian are given by ( = Px(Vf(()) and (ess = (V [U - UŞX T V (s))ss.
Doubly stochastic manifold. The doubly stochastic manifold (or coupling manifold) between two discrete probability measures μe, v with marginals a E IR"", E I2" is the set 11(,, v) = {T' E Rmxn T > 0, rl, = a, r 1m = b). It can be equipped with the Fisher information metric, defined (U r = E., (1gVe)))
for any U, V E Tr,I(u, The retraction is given by Retrr(U) = Sinkhorn(T's @ exp(U @ r)) where exp, C Ø are elementwise exponential, product, and division operations. Sinkhorn() represents the Sinkhorn- Knop iterations for balancing a matrix
43
B Important Lemmas
Proposition (1)). In a totally normal neighbourhood C M, a function f :: :u + IR is p-geodesic strongly convex, then satisfies for all x, y U
f(v) f(x) (grad ExpEExppX))1) de df (x,g).
2
If a function f has L-Lisschitz Riemannian gradient, then it satisfies for all x, y E u
L d (x. y).
(v) f(x) + (x)) (1))
Lemma (163 26)). There exists constant Co such that for any y1 12. Vs E Uy E My || r 1 < Cod(y1. (2)2 V3) ly.
82 81
Lemma 3 (Trigonometric distance bound 72 28)) Let Fa Tb. Ie E S M and denote
a d(zs, x.), b d(I,, xe) and e = d(za, Lbb as the geodesic side lengths. Then,
2 < ςb² + - (zt), Exp=. (xc))
2(Exp=
VI-ID
where Ç tanh( ID) if K < 0 and C = 1ifx 2 0. Here, D denotes the diameter of U and denotes the lower bound of the sectional curvature u.
C Proofs for Section 3.1
C.1 Proof of Proposition 1
Proof of Proposition[] 7 By the first-order optimality condition, y²(x) satisfies Gyg(z, yz(()) 0 e My. Based on Theorem 5 in (25), taking the (implicit) derivative of the equality with respect to x yields Gz y [x] + H,,(1, y² (D)) [x))[x] 0 for any T,Ms. This gives
16

---

Dy*(x) (=)===-a)) - g( 8" (x))) o Gyngle, y (z)). Notice that Dy T,M2 + (,()) its
adjoint operator (xy" is derived as follows. For any E T,M2, U T (w) ((Dy"()))* u). =
(()) (H, 5 (z)) y H y (()) (1
(oapla y' (z)) Hy gl (x) ())l
= v (x)) H. g(,, y (r),)
where the first equality uses the definition of adjoint operator and the third equality is due to the self-adjointness of Riemannian Hessian (inverse) and the last equality is due to Proposition D.2 in [27 that 22yg and Gyag are adjoint operators. By identification, we have (Dy"(x)) = - (z, y (x))
H,'' g(2, (x)).
Finally by the chain rule, obtain (from the definition of Riemannian gradient), for any M, (GF()), u). (9-f(x, y (z)). u). y,ff(zy (x)) (=)[u]]
(9.f(z, y (z)). + (9,f(z, (x)), Dy" (z)[4))
(()),, ((Dy (,))) (x)),
(G.f(x, yf (z)) Gy ) (400 (x))l
By identification the proof is complete.
C.2 On Lipschitzness of gradients
Proposition 3. If a bifunction (z has L-Lipschitz Riemannian gradient, then it satisfies -
S 2Ld(21,22) 22), where we let 2 (z. If an operator G(x, y) My TrM2 is p-Lipschitz, then it satisfies ([]((1) g D1 [|z I p d( 22) an operator
H(.,y) v) : My + TyM, is p-Lisschitz, then satisfies [[((2) ry H( bn Ily S pd(21,
Proof of Proposition 3
From the definition of Riemannian gradient of product manifold we have
[[]f(z)) 1 ((>))l = [[ G.. + +1,g/(21,,,)) 12 G, (22,
< II M1) Gzf(
[9..f(z1, G. (x2 1/2 Iaa
+ [1,,f(11, - G.f((222))) + 11, v1) 17 f 1/22 1ly1
I Ld(y1, 82) + Ld(x x2) + Ld(x), + Ld(1).22) 2Ld(21, x2)
where we use triangle inequality of Riemannian norm. Similarly, for the other two claims, we verify
[[((z)) (22) = |19(x1, - I ll. + [[[(11 12 G r 1a S pd(yi 322 + pd( x22 22)) The same arguments also hold for H(, v) and hence the proof is omitted.
C.3 Boundedness of ingredients
Lemma 4. Under Assumptions
2 we can show
10 1)l = 10 v) S L holds for any (x, y) E Ur x Ug
any [Dy (x))ll)) S Kι and [Dy - (x1) Dy (T2) (mn) S Lyd((11,2)) 22), any
(x2)
x, 2,,1,, x2 E Us where we let Ly = x₂ p + 2n1tip Kip
f
3 diy (21), (12)) < kfd((1 22), for any I1, x2 Us
For any x, F1, x2 E Eu V. V1, 82 E Uy
Ke =d v2). <
Jh H. g(z, v w - H. g( w2)ll
Kg
| v) - H 1g(gg 1) 1y s
/
17

---

5 Let LF =
+ Then for any 21, Uz GF(x2))l22 S L
x1 (F)))
Proof of Lemma 4
4 First we have for any v My
|| (z. Exp, (tu)) Gug (z, (I, 1) lim 6000
[llellly
Pe Wlll = [Dy9,g(z, [u]l -
lim 1-0
Lllollu.
[4]
|4]
where we use the fact that d(xx,,, (5). y) |k The operator norm is the same between g (x. y) and (z y) is due to the adjointness. This proves the first claim.
We first verify Dy (I) can be bounded as
L $
DDy (x))l> = 1765''(a, 9(2, y
and Dy (x) is also Lipschitz as
[Dyy(x)) ((
( (12) (on) " ( 90 (21))llr(()))
I 117 9(.1)' (z)) 5 (x
(2n
117,)"9(22. |1 (m1) g (
(o2)
(12)
< 4 9(1,,")( (x))) H g(
+ L HI, g(02 (x1)) r H g( y (21)) " (ل (22))
+ - ( (f
(2 ( - (21) + || y (
Cn (x1) - g 002, y (f2))lly (a))
Lp x2)
Lp
d(y" (z1). (x2)) (x2)) P d(l x2) L Ld(y' (z1), (x2))
112
n2
Lip
2pp
E3d 2²
Now suppose we let : [0, Mg defined as c(t) y" (7(1)) where v : [0, + M2 is
geodesic that connects Iı, 22, i.e., (0) 21, Then
((z)) = [.'' (t = / |DD (10)6 s i1 1% (n'(0 = d((2)
P
where we use the fact that the manifold is complete.
For the second claim, we first notice for any (invertible) linear operators A, B, A-1 A B-1
B A)B-1 and thus || - B-1|| < ||l
B|ll B-1 || for some well-defined norm 1 b2 = Hyg(1, y2). we have
|.. Here substituting A = m Hx, y1
H g(x, yn )r - H. g(0, v2) 105,,,99,,,,00 |1Hty g(z, (1) 1l y Hyg(z, 12) va 11, g(2, y2)l
(2)(1)
where we notice H. (I, ) = Hy vn) v1 II and use the isometry property of parallel transport. The same argument applies for ||Z--'((a 9((1,))) - H-1 9(
have
(F(()) - GF(22)l22
[T - -n,((E2, (12))llz
(œ1)) P2yg((z (2z) (e 9(2))) v² (zi)))ly (1) (
(x1)))l
+ 10 (x²)) le, I g( v' (z1) - (02,0(222)) v" (2r) o (m)) (on)) ( [1r (zi))lls
+ 1 12 (12)) l,, ||,,, 9(22, yf ( (m))) (1 (z1)) Gyl(22, 1/¹ (12)))l
18

---

From Assumption 2 and Lemma 4 we can obtain
-
I 19ef(2z1, (a1)) - 9.f(21. (x2) las + -l G. G=f(22, (02)))"(22221111 (22)
Ld(a)) { + L)(o))))))
L²
Ld(y
[nv 9y (x S25(2221) (x2)
(((2))
|Iry Cuf Gzf(21 ( + ( Gy (22, ( (22))la
12 ( +))d(,),,))))
Similarly, we have
(x1)) Png(az v" (22)) (12)
(0)
< 1 (x1)) - Gydlar (12)) ( + Leyt,, - Pzyl(v2, yf (r2)))laz
(x2)
م p)d(01122)
p
and
J s-(2)) g(1 (m1)) - 9 9(22, (x2) r (1z) (
|| ((1,) (z1)) Hn''g(22 g( ('1))l (m)) || ( ( g(02, v' (xi)) (xt) "z'9(x22 v" (22))ll
(
$ { l + plz
112 13 (
Combining all the results together, we can show
[[)gF(>)) - (F(2))l2e $
+ L + PL
+ p) M
+ LM
P 112
+
pL
+
L2
+ 1.) d(2))
11
PM + 4² + dd(11, x2),
PLM
I
+ 1) (L +
which completes the proof.
C.4 On strong convexity of the lower-level problem
Lemma 5 (Convergence under strong convexity). Under Assumptions 7 suppose ny 4
where ( > 1 is curvature constant defined in Lemma 3 then we have d'(uttt .v" (xx)) < (1+ (1+ myU)dd (yx, y (IR)).
Proof of Lemma d d .y" (xx)) < & (vé. y (œx)) + + vi)) 12 + 2ny vi), Exp, v (2x))
We apply the trigonometric distance bound from Lemma to obtain
s (vi. of (xx) + n
+ 2ny (g( y/" (xx)) (( vii
11
d2 (
2
< + CT? n/ dd (vi,, "² (xx)),
where the second inequality is due to geodesic strong convexity and the third inequality is due to ( v1) 119,g( vtt) p (=,, 9 Gog(ak. y² (Ix) I L y (xx)) and the fact that y" (Ix) is optimal. Here, require 7/y thz in order for + n2SL2 Myp < 1.
19

---

C.5 Proof of Lemma 1
Proof of Lemmal
!
Hessian inverse: for the Hessian inverse approximation, we let
= v) Liza(a, v) (1'p"(z"5,, ()),,
It can be seen that (F(xx) Gf(x,, (zx)) and F(an) = Gf(z Then for any r Uz 1/, 4/2 E My. we have
(0)(,,)) - (xx 12)1lx $ [19. v1) (x, (2,, la= + B vn) - 0x) (0, Ml 19 (z, vi
y2) =- |T 9 91 vn)) H g(0, 122 By 32)]
p.M.
$ (L )d(y1: 122) 3 g( v1) Hly (x,, v2 Ia 9 (i)
L LI g(x, v2)
1I, / 31) f v2)ll
$ (L pM L2
LMp
d(y1,
1
n2²
where we use Assumption
2
and Lemma
Conjugate gradient: we let vx Hy" g y² (y v (mx))] E (yr and let 0z
Hy g yk+ ( (x E T My. We first bound
G ((x)) (F(x))l
y (ax)) (lx, + o ue+ ur+ 1llll D " $ G0L (6+1) I y
VEllbwrs
+ 1 1/ (xx)) g yk-
(m() ala 1
I Lddl" (IxR): 11++1) + T1 - T (mx) "llbenn + pdly (Ix), (oR
s (L+ + (..M)d(y (Ix), 1 + 1 (x.)) oflbm
where ||08 < M/n. From standard convergence result eq. 6.19 in
we have
1)
27
116 02
4
+ 10
0g
/Kj
This leads to
oo
B
(m²) 2k lill ve+1
(4)
$ |10
0g 19
+ ( (mx) oz --00 0illami
11
2
151
1102 - Willmass + I1T (ox) 0illmass
/Ka +
V Kl
1
1 ( all + +1++1 Vkl
$ 22
|| (a,
VKl +
Vni
KKl 1 T
Mnp (dd( Met1), (5)
$ 22/
T o" (#x) alll + (1+2/Ni) (N1 + #
VKL
where in the last inequality, Combining the results yield the desired result.
the definition of or and of and the Lipschitzness assumptions.
Neumann series: Hl) = 7 Ei='ii (id "Hygg(x,,,) Then we can bound
((TN)( - (F(w))lm
< Ldy (xx), yx+1)
+ +I ( vi+1))l Hu (V6+1) - (aa)
- Ho 9(22, ((x)))""(=)))) (++) Il 11 [xx, 11 bi-1
2 ( v1+1)l H7, 0(xx,2) v (Ix))llv (ax) (l [|l (mx) ( G, (2x, wh+1) Gy 22 (xx (lv (m))
20

---

+ IZ Uk+1) - Pz (xx)) (xe) (1 11 9( (xx))ll (xx) ff y (xx))lar(s)))
$ (L + k.L + M))(y (k+), + LM I (1 - (mx) Hy g y (xx)) (2x) llw---
Mk41
We now bound
H (v++1) - 18643 (mx) Hyy g( y² (xx)) (xx)))))))) Ilves
ya+1
$ ||Ftx(UK+ (Uk+1) -Hly g( " 8k+ + Hyy g( (k+)) -
(mx) Hog - 1
g(re V² (xxx)) (xx)) (xx)
ua-1 B1ues ae
(id (id- Hyyg Uk + $ 1h Σ
Re "ey (xx)) (xk),
1=
"e (1-2H²"
dl (xx).
P
where we use the lower bound on Hyg(ak. Mk+1). Substituting the results back the bound yields the desired result.
Automatic differentiation: Given ++1 = Expy (n vk)), we can show its differential is
Dzp Us = 1/4 (Dor Vk - my Ny vk) - 7h H H sf Dec + Eg
P y/ (lid my Hyg // Da mo g vi)) + E
where
| y
S Cssl(id N ( af ny Na wie (f yf) llve
$ C3 C1 my (| (ft 1
In addition, we notice LaaF (xx)) Gzf v + (D.. v (9 vi) and we can bound 11 F(x)) (F()))l
| uz) Gos(Er, (xx)) + =1((aav))) (Dzc (1x)) ( (wo) Ioa (n vi Ille
+ ||Dmyy" (xx) llv' (#x) || Vz
(xx) Gys(En, vZ) GyI(ER, v² (xx)) (y
Ln)) + niLLL m & d(y" (xx), v + MI (mx)
(mx)) Daay" (xx) Dan ur 1
(6)
where the second inequality uses Lemma 4
Then we bound
v((z)) Dzay" (xx) Dz "
1/
$ 21 (xx) DzaV" (xx) - -P ((id - n, H,gg vt)) - no To 1lg
C3( 1, + ff ylll
I ((e) (Dag (Ix) + , Hygg '/t ( (xx))) (*x) + Ty v² (xx)))
P a (lid - 11 H F4. 1/ D - n T ( (k vi 1lggtt
+ (0 C + 9 llve
||r r
vi (2x) (Dz.1' (rx) n, H (zx (xx) no No 8 (œx)))
ro
vi 11 p (id n, (
ny Sgn (i(...//)))... v1))
+ Cs - C; +
vi llvg
me x (zx) (id H ( (2
(xa)) my The nut
(mx) (ox)) Sna v/ (xx))
((id - m1 H (x N v - n,, S2 G (o w )1)) (1 C2 + C2 1))
+ +,,)))+++ 1n, v llvl
T (id - Hy 1/" (rx)) (f))) (#x)
r (1x) Dea
21

---

+ Ty (Hyg((a vi) 1ht (=(x)) H y" (xx)))) (1) Do + to vk) 1 S y (=x)))
+ Cn + Ca) (11 Byll + 7,, [) 119 (T l2ll
S - (1-7/(]] 7
4¹ (ax) Dayy' (Ix)) - Dill: + my (M) + pu v² (2xx)
C2 C3) (1 C/ + n,
v
$ (1 I y/' (mt) D (xx) - Dz, uy + (m( + 1)p(1 ( v² (xx))
+ my (C2 + 7/0 Ca) ((1 C1 + 1,, 2f)l
$ (1 - ( (mx) Dayy" (xx) - DD2, 3klle P0
+ 1/4 (Ni + 1)p+ (C2 + my (() - Po + 1n 1.) (1 n/. - Myy14) d(y y (xx)),
where the first equality uses the expression of Da,y" y (Ix) (Proposition D.. The second last in- equality follows from Lemma 5 and the last inequality is due to the smoothness of Riemannian gradient and Lemma 5. i.e., ( (Tx. VEll = |
9- ( gul y (xx)) - guf(ax lbl S
Ld(yie (2x)) S L(I+ + nZCL2 To 2 d(yh: (xxx).
Finally, applying the bound recursively, we obtain
| y (2x)
Dz, (I*) D2, vlll2 S (1 myy1)) || y (2x) Dz, (xx)) Dz,, vblls
8-1
+ C
Σ
(1
S 1-1-1 (1+ m2sL2 - ny d(y). v² (xx))
==0
S k1 nTp1 8 1hy Σ
(1 + LL² - n1p1))811 $ dlyb 32 (xx)
8=0
(1 + SI, nyll)) ليع
S (k moluu 8 + mh Ö
dlyk: (xx))
1 - (1 + n2gLL 2 - nyl4))
S <R((1 - no11) +
11 S (1 + sTa 77114) 8 dlub y/ (rk)).
where we let Ö = (Ni + 1)p + (C2 + 7y Ca) - C1 Ci + myT) and we note that Dar = 0.
Combining the above result with 66 gives
|1 (xxx) - (F
(aa) $ (L + LK)) + (1 m1yH4) d(z" (2x), vh) + M MIr (m Dog V (xx) - D Uh live
2MC
TH < 41 T2 + + + ki)) (1 + n2 nulk) - d(y.. v (xx)) + MMi(1 (1)))))
where we use the fact that + 5L2 1/H4 S
C.6 Proof of Theorem
Proof of Theorem
By smoothness of F(x) (Lemma
45 we have
I. 1GF(#x))13
F( - F(xx) (Zx)
- Fl - LF)
Me + [[)
(7)
Now we consider the different hypergradient estimator separately.
1. Hessian inverse: Let Chiny = L + R + K2L + KI
[xF))))) Grany (xx)) < Chin df (y" (Ix), 11+1) S Chinv (1+ + n ny11) (x)(x))) yk). (8)
22

---

where we notice uk+1 = vi and apply Lemma 5
Furthermore,
d (xx))
2d (vk--111 g 10 (Ix- + 2dl (v" (Iz), y (2x-1))
$ 2((1 + nggLL* - 11144) (dd()(((" (Ik- 1). yk + 2n ( F 1xx
21 + n (1111) s "d(" (1) 11
+ #1 1)1 11ae_.
2(1 + n Tyy S of (y¹
4n P (1 + nZLL² dd (Ix-1),
+ x F 1)11
=2 + 2n xi Chin (1 ngLL² ny ( 9k-1) + ( (9)
where we apply Lemma 5 and 43 in the second inequality.
Construct Lyapunov function Rs F (Ix)= (yk 22 (xx)). Then,
R&+1- R
F(x++1) - Fxx)) + (d (yk+1, 1 (xx+1)) dl (yr, v" (xx)))
s - TE 2 LF G 11:s2 + (/2 + n Lr ( F(xx)12.
+ ((2 + 4n 2 blnv (1 + +ng5 - nu1us 1) 1)d (v/ (k) 4n | (GF(()))22))
s (h Lr 4n x2) || (5x))
The
+ (2 + Chimel 1 + 1/m Lr) + 4n C Binv) (1 + r2, noHis - 1) dP (xx),
s The 5/ |19 [[
2
+ (2 + Chimv
Th _h2 + 5n Lr)) + n22 Tyu1) -1 dss" (xx), yk)
where we combine
7 and (9) in the first inequality and use xi S LF in the third inequality. Now setting The 20LF we can simplify the inequality as
RE+1 Rs
1
3C hins (9F(E)) + ((2
)(1 (1+11) + - - (xx()))) mn)
80LF
807.
1
S
[
80LF
where we choose 2 log( 1601 1 )/ log(1 + -
80Lr
-1y14) = =(x7() for the last inequality.
Summing over k = 0, O...K - 1 yields
K 1 K
80 Roo RK)
< 80L.p
[GF((x))12
K
K
k=0
which suggests mink min.= K- 1 || < $0
Conjugate gradient: Let Ceg L + Ri.M + L(1+
+ 12 2/mm) (m +
Mre). x
Then we can show
[[9F(I))= GegF((x))
1 2T
< 20 (1 + 1 - mop) ( (xx) +x)+ 882
1100
1 11 o" (ax)
12 B (10)
where it follows from Lemma Tand Lemma 5 case
Then following similar analysis as in Hessian inverse
d ( (Fx+1)) s 2(1 + n s (y (Ix), Mw)
(F(x))
23

---

+ +4
(2 + B
( +
( (x)(()))) yk)
22 1
321 m/ L2
l60 (Ta) Villert +422n111gP(C
v ke
+
(11)
Further, noticing 0% as
1 67 k--1 we bound ||0]
1 (22, -11% - 1
y" (mw) lfllve
(
1102
p (xx))
0x 12llya11
(mx))) Vlll + MCOD d(yx+1, (xx))
< 1102-1 - -
ull () + MCOD
S 110T_1 1
I 4" (1=-1) ) // llve + -1u2 - ry (xt) (o
y (0x-1) 2
d(ye+1.""
T 4¹ (1x-1) 1611 + +1u (((x ( 1i- MCOD (d(ws) (xx)) + d(ykt1,y'
$ 11
/
1 17 "1008
1llln 16) + (1 (Ma)(Ia Mip ydly (sx-1). yk)
S 22/KKl Ai
T (o
2MC.D d(yk, y (c=))
+ ||0z
(are)
+
1a
25/
VKi +
"1108
Tye
(ax-1)
( + 22/m((m /
Mr (1 + n//TT? -my1) d(y" 1). yk-1)
+ 110g
pg (Ira)
08-1
+ 2MCoD d(yr, (xx))
(12)
where we use Lemma] in the first and third inequalities. The second last inequality follows from (5)
and d(uk+1, (zx)) S d(yi, y (Ix)). The last inequality follows from Lemma and nq 2 1. Now
we bound
lluğ - (m))
1) "
g( (Itk (9 "|717'0(0,,0
(xx))) pV" (#x)
yn' (ox 1) Hyy g( 1, yf (a (Ik- ( (fx- (xx- (mw)
I MI 9 V (mx)) - (ar) ((x Hly g(Ib-1, 8² (xx 1) 1)) (x())) (1
1 (ox-1) G, "² (Ix)) Gy 1 1/ (I (1)x)
$ MI ( V/' (ax)) - H. g yf (Ix Myy (xx)
+
x1x 9( (xx)) r (o (mx)
(ou H g v (x (c)-1))) (-()) "lly-(
1
pry (an-1) 10 Gyf( (xx)) gy (Fx. y² (Ix (Ix-1)))
P
1 (xx-))))(( Gyf (x (1)
I
+
J
L Mn kid( In 1) + Sxd + - d( #x-1) MKg d(Ix, xx-1)
P
1
/
Co ( F(x - (aan))-F(aaa (F((a=-))l + n2. C, FF
(13)
where we let C. = Mra + MaaRi + K7 + Combining 13) and 12 we obtain
1100
2x
1 2T
(a 1 112. + 20kz (ma + Mrp (1 RTT. " (IR-1); 11-1)
20kl
1000
D
/kl
+
24

---

+ 5n (( - (I + ( + 5M
dl a'n yf (rx)).
(14)
Now we define a Lyapunov function Rs == F((x)) + &(vx. (xx)) + 110g - T (mx) U/ l 2841" Then
RE+1- - -R
= -F((Dx1)) - F((x)) + (d (16+1, (xx+1)) - d (uks 3 (xx)))
+ (116++1 - -o (((++))) / 1la - 11 - 8 (zx) u / ()
S 2 ((LF)) l1: (/200 + t [[GF()x)) - FEx dd y (ox+1)) d &yvr 3" (xxx)
2
/ki- 11-1 1, 2T -1) b/ - potss (xx) ( llll.
+ 20kl(mi + M + n - mn11) &f (x)) (1)
+ 20kz
/ki +
14
+ 5m C. Gog F (xx)) - (F(IN)) + 55 C. + 5M2CRD
&(yN+1,8' (xx+1))
122
Th2 11e Lr - (GFF (xx)) + Th2 + bb + (GF(())))))))))))))))))))))))))) F1x
2
SMCRD D2
+ dd 11)c((1)++1)) ( y" (xx++))) + (20ma(m + Mgg
+ ng - -nnu0s -1) dl (Jx), 9k) (1 2
+
u²
1
Vki
1 22' -/-00- -1)1 - - ( oiller
+ (20ml(
/Kis +
W LF -5n - 2n2 5MCDD² + +1)) (
Ur
$
2
112
+
+ T
1.
+ C) (9F(#*)) - -eg
2
55M
+
(1) + 20kz (Ni + MRg )² )²)
(1 + nisLL - -7g1)5 - -1) (ex. v (xx))
den
112
M
+
(322 31 L
5M2CRD2
+ 1) + 20ki)
VKi -1
/Ki + -1,07 2T -1)1 -
-r (((()) 1 "
112
A
[1
(1 (22 + An ) - Geg FogFF(x F(xx)
+
D2
+ +1) / + 20<z(11 + MKoy ) ))(()) + nisLL² -nau)s 1) - (yn.
of (yn. v" (Ixx)
12
1
(322)
+ +1) + 20ki)
VRai - 1 2T
1/00 -1)1188 - -Mit
(xx) ( 1E
122
/ki /Ki+1
- [GF + (3327) [2 L2
5 C3 D2
+ L2n2) )
/ -1 22 -116 (m2)) o ofllen.
Vki + ,
+ 1) + 20nz
96Λ
112
4A
+1]
8n + 20m
+ M 2 + 16A
cg
(1 + n/LL2 - -nghts - 1) (l)de(ee (2)))
112
P
-
119F(xx))
- 96Λ -996
where we 77.
in the first inequality and (1] in the second inequality. In the third inequality, + SM203 D²
we let A = C2
+ 1) and because I.p O(x?) and A (x+)), ) we can without loss
of generality have LP S A. We also choose 7 = 24x
- and use 10 for the fourth inequality. The last
inequality follows by choosing
5 2 log (
(5MCgDD
+1) + 20 +
Minpg )) /log(1 + +iggI= - 70/1) = -1/h)= -
122
P
3 L2 (32
5MCRD2
+ 1) + 20mi / log
1 =Ō(vm))
T>
log
2
42
Rq +1
25

---

Finally, telescoping the inequality, we obtain
K Σ [9F(()))l
1
96A Ro
96Λ
(F(xz) + (yo.v" v (xo)) + || ( K
K
where we use the fact that og = 0 and the isometry property of parallel transport.
3. Truncated Neumann series: Let Cus == I + kiL. + Mp + R. Here we notice that
Cus Chiuv- Then by Lemma we see
F(rx)) (F(x))) 220 - 5T, - db (1 2n/ M - -11)22 (15)
Similar in the previous analysis,
of (y (IK+1)) S 2(1 + n 11 (Ix)) + 4 (F(x)))
+ 11 112x
(16)
Let the Lyapunov function be Rg = F(xx) + da v² (xx)). Then
R&+1 Rs
< -"1 - - - (1E + +( Uh + n + (9F(())) -
2
(2(1 nnsI. -nop) - 1)d y" (zxx)
1
30 us
$
||9 (x=) + ((2 40L
(1 + -2 nyH13 (yk.
80LF
12
2 - 7 2F
40LF
where we set 1/2 DULF and apply 155 in the second inequality.
Now setting S 2 log( +302 )/ /lgg(1 nz (Tpp)) =((x(t) and telescoping the results yields
K 1CF((x))1 I| - K
K
801 Ro
+ GnG M2 - - S
80 Ro
+
6
K
2
where we set T 2 log( 12mr log(11 711) (olg)))))) 2))
4. Automatic disferentiation: Let Cad == "-- + T.(1 + kt). Then
2M0
(yx. [(g)((.)) Gan Fo S 20 + +12 ny11)s --1 (d y (xx)) + +242) - 12,11)25
and similarly
d (16+11, (#x+1)) s 2(1 + fi - nHH1 (2w)) +
(F
+ |
Let the Lyapunov function be R2 = F(*x) + d² (y 1² (xx)). Then
R&+1 - R& s - Te - n LF 4n m2)) ( (1x) + C + nDLF + 4n - F(x))
2
+ ((2(1+ + 11 - n1144 S-1 1)d (yR, v² (sw))
302 ad 40LF
1CF(( + (112+
)(1 + n - - i)d2(uee 1² (xx))
80LF
3
M2x (1 25 25
40LF
3
s 80LF
(9F((x)))
M2 x² 11y1)225
40LF
26

---

where we set 7 = 2 and choose S 2 log Telescoping the result gives
80LF +30 llgg(1 [ + n2sL2 -myH) + = 6(n75).
40
[GF(x)))
80Lp Ro
+ +6M2xx (1 - K
80Lp
+ 2
E
K
K
by choosing S 1 log( 12
log(1 my/1))
Hence we set (11 (og
for both conditions to hold.
C.7 Proof of Corollary
The computational cost of gradient and Hessian for each method for approximating the hypergradient are as follows.
Corollary 1. The complexities of reaching an c-stationary solution are
Hessian inverse: Gy=O(k}e="), Gg = =(KKGe= JV, = 0 HV NA.
Conjugate gradient: Gj Ol1 G3 Ö(nf(e-1), JV, O(ntc-1), HV,
Truncated Neumann series: G₁ = Oxxe G, = "( JV, = O ) HV, = ÖKK log(=)))
Automatic differentiation: Gf 1kk)t-1), Ö(HicE -oo(( 1)),
O(x((c==log(("))) HV, = (og(
Proof of Corollary [
From the convergence established in Theorem we see the iterations in order to reach e-stationary solution given by
(Hessian inverse) = OLLpE 1)
(Conjugate gradient) K = O(Ac = O(ffe S - Ö(KG)), Ö(
(Truncated Neumann series) K = = Ō(kz log(e-1)
(Automatic differentiation) K
Then based on Algorithm the gradient complexities are Gy = 2K and G, KS and cross- derivative and Hessian product complexities are JV, K,HV, JV, = KS, HV, = KS for AD (which we approximate based on the analysis in Lemma TT. We notice here for the Hessian inverse, because do not compute Hessian vector product, we write NA for Hessian vector product based on the Neumann series. This completes the proof.
KT for CG and NS and
D Proofs for Section 3.3
We first show Lemma 4 holds for each f( y), gi (æ. y). Further, the variance of the estimate be bounded as follows. We here use [-] to denote all possible derivatives, including x, y, zy, ya.
Lemma 6. Under Assumption 4 we have for any x, y E U. (1) [ v) 9 (x, w)lli < M
(2) (1,, 9) Gz v) L (3) LE H (; v) H., g(z, 1) <
μ -
For notation, denote the filtration Fz {vo, va, II,I,,I,...... II, Fk, y and here we let E& == EE With slight abuse notation, we further consider F {yo, To, U1, Fi, ..., M yl,,.., y and correspondingly let Ef- = [[]]
27

---

Lemma 7 (Convergence under strong convexity and stochastic setting). Under stochastic setting and
under the Assumption that is geodesic strongly convex, we can show R&d (vk fu6+, y (rx)) S
(1 + I (1 1" (*+)) +
BB.TT and 1g-1d (UK+1, (mx)) S (1 + ST
7114) (yk, y (xxx) + 1 |811"
MoLM2 1.
Lemma Under Assumption we can bound Es (()))((*))) - (F(xx) S 16M
4M +16M2
18g
8M [Ba1
+ 2C himy de (v y
Bil
D.1 Proofs for the lemmas
Proof of Lemma Here we only prove one and the rest follows exactly. Due the unbiasedness of the stochastic estimate, we have
El ( v El|G.S((x,)))
/ 2111 - (1 = M
where we Assumption 4
Proof of Lemmal
7
Similarly from the proof of Lemma we take expectation over F
( s-1 (zx))
LER[d [f v/ (Ir)) + n Gu9e, (oxs ui)l 2n (9y98, (1x vk)) Exp, 2/ y (ru))
s df u (xx)) + EER 99 (Ies, wf) Gyg(Ees u 11 + 1 wF)
+ 2n (9 ke vk), Exp "/ (xx))igz
< (1 + ngrL² - u u )d (vi. 3² (xx)) + EE
Σ
E 1/6) - Gyg( whe
13, 2
iEB,
CM2
$ (1 + nz - )y (vi. 3² (xx))
|B1
where we use the strong convexity and the fact that Blggyo, (z, - gug(a, vll 1 Bl| Eies, (I, Gyg(æ, 1) 2 Eies, Bllguar(x, y) Gug(2, 1) in the third inequality and Lemma [ the last inequality. Further, we telescope the inequality and taking the expectation E gets
M2 S-1
ER 1e_ (d (xx)) < (1 + n22 n1p1)² &f & (xx)) +
Σ
(1
+
++ 11,1))
B
==0
$ (1 + CLL nup s d( (xx)) + C2M
1
# S2552
|B11'
where we use the fact that E
S-1 0= S 1-8 for < 0 1.
Proof of Lemma 8 Recall
that GF((x)
G=f(z,v(x))
Gzyg(r, 1² (x:) H g(x, u² (a:)) (a. (x))]]] Then
Ex (((z) (F
< 2E 2 Griny F(x)) 2(1 -
< 2F (F Gnime F (xx) 122 + 20 & (Ix))
(17)
where the second inequality uses Lemma l12, follows.
Now we bound the first term E ExI(@F(xH)
-
Gnime ((xx) F
First we bound
Ex H|, 98, (Zx, 1k+1)l9y (#x, Mk+1)l H. g(( Me 19 (Ix. we+ 1ll
28

---

- 2 [| ga vk+1) ( J9 (Zk,, yk+1) of Mo++ I11 10.2
+ 2 g8, (ot, yk+1) H g M++1))l9 fb, (2 yk+ 11l
- 2 H,,
3 "((* yk+ I1 +e+ E& 1 (I yk+1) J (2 yx+ 1 pat
+ 2B 9B. (T 1 Hy ( Mc+1 Ee (I 1
2.
12
Ba
where we notice that [| (0&s yk+11)l $
1 [8a1 Eies, Bh+1))l S M.
Hence, we can bound
Es (FF))) - Gniny F F(xx)
8M
1
1
< 2F 11..119...0 (Tk. 1k++)) Gr (Irk Uk+ 11l +
+
R. Paygs (# Uk+
112
[22
[Bal
+ 41 11 (T Mk+1)
Pey (Ii 3/+11)
1 117,, g ( yk 19 (IT yk+ (ortlee+)
2M2
1
xM
<
M
(18)
[82
B2
14
B3
where Lemma in the last inequality. Combining
18 with 17 yields the desired result.
D.2 Proof of Theorem 2
Proof of Theorem?
From the smoothness of F(x) (i.e., 7 and taking full expectation we obtain,
E[F(TL+1) - F(xx)] s -(( ("T Lr (ELF))((G((WWIIE (T + ((G FF(()))))))))))))))))() (2
Further, we can bound
( (*x+1))
$ 2Ed (ok+(1u11,)) (xx)) 1 +
22
1
< 2(1 my n 2 nuH " EEd (yk: y" (xx)) 1 565L
[1]
+ 4m ((x)) (F
where we use Lemma and
g in the second inequality.
Next, we construct a Lyapunov function as R& = F(*x) + d2 (yk. y² (*x)). Then
E[ - R&]
s F((((++11 - F F(xx)] Eld (16+11 1/ (1++1) - & v (xx))
$ - -thh 2
The - LF
g (xx)) Ie,, + ( 1e + L +
(F(xx))
-
(2(1 -iz5L² -1y14s - 1) Edd (u 1² (xx)) + 21,
11 11, 1. JB,
E|(F((x*)11 + 80LF
E E (F (E
SOLF
+ ((2(1 + TT² - 1//14s - Ed (m). (xx))
no5 |B11 (11)
E|GF([=))1 + ( 4M 16M
&M2R
[B&l + 16M²
3C Edl (y6+1, (Ix))
s SOLF
80LF
[22
|Bal
40LF
+ ((2(1 T2 - 11111s - 1) Jed (yn. yf" (xx)) + M2²
11 T.2 [B,
3C2 hiuv (1 + 5I2
B|GF(<=)11 + ((
Thy11) -1) Ed² (yk, 10 (xx)))
-
80LF
40LF
29

---

3
AM2 + 16
8M
16
3C
7g5M²
+
+
2)
S0LF
[B2]
Ba
|B4a|
40Lr
S/y |B31
4M2
16M k
8M
16M222
3Chin
11yMM2
S
+
+
+ 2)
80LF
80LF
[22
[B3]
|Bal
40LF
[31
where we choose The
40L = 201F the first equality and 2 log( 500.p+330h
hin )/
log( 11 -12111111))
(K}{) for the last inequality. Telescoping the result gives
1 K-
80L, Ro K
12M22 48M2??
24M2
48M
Σ
v
+
K
1B2l
[Bal
B1
0=0
Sn5M²
+ (6C
+ 160LF)
11 17 L
|B1
80LF Ro
+ 2
K
where the last inequality follows from the choice that |51] 2 (24hhnn + -640Lr) SnySM²
p- CT,2 //////E
O /c), [B2] 2 144M/+576MM²2 (nq /c), [Bsl
288M² mi (kk |Bal
576M3
((x7 /e) in the last inequality.
In order to reach e-stationary solution, we require == complexity for f is Gy 22KBB ( and for g is Gg KSB Onp The complexity for cross-derivative is K K|Ba|= O
= O 1) and thus the (stochastic) gradient
E Proofs for Section 3.4
We first give complete proof for the Hessian inverse estimator follows. For the other estimators, we only provide a proof sketch.
Proof of Theorem 3
Proof for HINV. (1) First, we derive the convergence under strong convexity using retraction. By the trigonometric distance bound
d ( y" (xx))
I (vi v² (xx)) + (vé 11 2(Exp, wht Exp, "² (mx))
S df of (zx)) Th Gug (x v 11 2(Exp vi
Retr, bt vit Exp on (#x))
2n (0 v Exp y (Ix))
I dvt. v² (zx)) + < g xx v 11 + 2n (Gvg(ak, vk), Exp vo (xx))
+ 22 U/ +1 Retr
< d (vi. 3/ (xx)) // (Ce 2 2Dc 11 Ville 2n, (9 vi), Exp (au)i
2 (1 + ni (Cc cc + 2Den mn1) d v
where use Assumption 5 in the second inequality and fourth inequality. We require f/y (a+zDeR))= in order achieve linear convergence. For simplicity, we let 7 jmy n (CT+ 2DcR). This leads to d (2y+) yy" (zx)) < (1 r)d y² (xx)).
(2) Next, we notice the bound on hypergradient approximation error still holds as Gnin F(Ix) (F(xx)lla> dyy (x*x), Vo+1), where Chtav L + KpM KiL + M.p.. Further, by L-smoothness,
F(I++1) - -((EK)
I (GF((x), Exprn) retida + LF d (Ix, Fb+1)
30

---

zt 2x+1) - nn (GF(())) (F(x))) ELFIE 1| 1
(FF(xx), Expe 2641 -Rer
(2KMMRR chr ) ( (
I M4R +zL +zIr)) (FF - -FF(x)) + t (
(gF((x)) G The 19
=- 1/e (4n:Men eLr) (I9F (22 (4M)Mn + c ( (F())))
where in the third inequality, we bound [[]F(z)))m,, S MM + LM S 2LM
x
(3) Then we can bound
dyn,, o² (xx)) &
(xx-1) + 2d (v (xx), (#k-1)) 2 --1--7 7) ( (zk-1), 9k 2 F
2 + C (1 - 7) ( (xx 1), yk 11-1)) + (0
where the last inequality follows similarly as
Let Lyapunov function be Rs == F(xx)= + d (vv. y (xx)). Then
R
R
T/e
(4KrMeR + 2LF) ( [[G + (7 (4R + ZLF) [[FF(*)) ( -Fp((=t))(((x)))))) (
2
+ ((2 + An Crinv (1 - 7) - 1)d (1x), mk) +m2
(AmMMCR + 2LF) - - | [9F
2
(2
1 The + (AR:MRR + + + Piny (1 - 7) -1 -1)d(ry (y)())))))) Chinv
The
Lrv²)
|[GF(T)))E ((2 C The L
-1)d² cf(2 (Ix).
2
< |[FF(*)))
16LF
where we use k7o S LFd and let LP = 4k 50LF in the second last inequality, and we choose T/e प,.r S log (.2T. ) log(1 r) @(x?(). in the last inequality. Then telescoping the results yields Finally, we sum over k = 0,... K - 1, which leads to
1
K-1
ICF(ON).. < 16L K
Σ K
Thus in order to achieve e-stationary solution, we require K O(LFe-t) O(k}e-1) and hence the order of gradient and second-order complexities remain unchanged.
Extensions to other estimators. To extend the proof to other hypergradient estimators, first notice that the convergence of inner iterations for solving lower-level problems is agnostic to the choice of hypergradient estimators, i.e.,
d (vit1 y (zx)) (1- (d- ( ( (xx)), T= -m n3(CE 2Den)
(1) Hypergradient approximation error. For hypergradient estimator based on Hessian inverse, conjugate gradient, truncated Neumann series, the hypegradient approximation error remains the same as Lemma] given no retraction involved the computation. That is,
CG: CG:11gePFaaz - (Fxx)) S Cogddyr (Ix), 11+1)++ 2L/ where Cog M + L(1 + (/) + MRe).
/Kc VM++1
1106
(xx) o I
31

---

NS: F(xx) - (F(()))()))
S Cnedly" (xx), Mk+1) + kM 7p)", where Cns = Chiuv-
For hypergradient based automatic differentiation (AD), first show that there exists con- stant C2 (that depends on C3, C, CR) such that (2Retr,(u) = (Retrs(w) (id + Dg + &. with Rottr I [<|/D=u|[=///]]ll Such a result can be derived by bounding the difference between retraction and exponential map. Then we follow the analysis for Lemma 1 as follows. Given u = Retry: ( (-1y Myg9y(Tr. vf)), we have
Dz, ++1 = P w ((id - Bly Hy Hog Dz, Vi nh Do v()) + &
where || 11 S (4 + my my L) ff lloa The rest of the proof follows exactly from Lemma 1, where replace the convergence of inner iteration with the updated rate. This gives
(5 (x+) (F(x))
2MC'
$ 11 Bhy (C2 2DCR)
+ LL1 + +x)))(1 (1))11=)) 8 d(yk v2 (xx)) + Mki(1 nnH))
where C²:== (mi 1)p + (C2 + Ca) -ngl1)11 1)
In summary
AD: -GF S Caa(1-T) d(y y (xx)) (Mi 1111--yph1 1y14), where Cad 2MC
A mylce+2Den) L(1 +K).
(2) Objective decrement. This part also the same across all hypergradient estimators, i.c..
F F(IR+1) - -((x))
1/5 (4m:Mn eLr) (IFF))))))) + (4mMMn + aL (GF(())) FF(m)))
2
(3) Lyapunov function decrement. The definition of Lyapunov function depends on the choice of hypergradient estimator.
For CG, we define R2 == F(xx) analysis, we first bound
(Fx)) + 100 D (=x) Vfll Then following similar
d (y y² (2x+1)) < (2 8n (1 7) dx (xx), yk)
+ +22 1
/kj + 2 100
(19)
and
1106
1 ex okllya
2 1100
(=x-1)) oh-1l + +2/Ti(1i + Mlre) (1 -) d(x""(xxxx1))) (Ix-1),
2 K7
Th
P
+ +1102
r (Ira y (zx-1) 2 1 +
2MCcD
d(yk.y" (xx))
P
Then similarly.
|lug r (mx)
(mx)
(o 1) "h (11 S C,d(ax,
(n -F((-))))ll +
where in the second inequality, we use the bound between retraction and exponential map as well as triangle inequality. Then combining the above two results, gives
||06
D y² (mx) 1x 1 1
32

---

20kr
/K|
1
2T
1
y (0x-1) uk- 1 20ki(ki
(1- d(y (0 1),-1)1)11)
+
1
CM
+ s
(
+ 5
d d (2
4²
Then we can show
R&+1
Rx
1 (GF - + LELF 50 ((F((())))) Goz. (swa)l
02222222
+ +1) (+))((1++11122 (xx+1))) + (20mi(m) + M - 7) - -1) (1ddd&// (2" (xx))
(11++19
1²
1
20mm
VKj - 1 2T
- -1)
11 -1)1 D -
VKj
+
The
G
[[gF((>)) + +(- T2 + (IGF)) - Gog ( (xxx (1xx)
2
+1) B1 Cag) + 20xi(<:) + Mig ) (1 + n - 11111 - -1) (sn: v/ (xx))
+
112
1
52
Vki - )
(32
TT
+ 1) + 20ki)
- -1)1183 110% M ( (o2) r
42
+
|IF
96A
D2 # and without loss of generality LF S A. The second inequality by (9). The last inequality is by appropriately choosing T, which on the same order as the exponential map case. Then telescoping the result yields
where we let A Cha+ kł{
K-
(GF(T)) < 96ARo
K
K
K=
For NS, we define Rs = F(xx)) + d (yk, y² (xx)) and derive d (M1++19 yf (x++1)) 2( -) 7 dd (y" (xx) ) vk) + 4n el Gne F(xx) GF(xw)l.. + 4n
Then we can follow exactly the same proof as HINV that
K
16 Ro (GF( I
+ 2
K
K=
by appropriately choosing S, T and 7/-
For AD, we define R₁. = Fxx)) + d (yx, yn (xx)). Then following the same analysis except for the choice of S, show,
1 K
16L Ro
(gF(x)))
<
+
K
K
2
K=
Thus the proof is now complete.
F Tangent space conjugate gradient
show the tangent space conjugate gradient algorithm for solving the linear system H[r] = G. Similar to [0], we set the initialization to be the transported output of oi"_, from last iteration, where o = 0, which is beneficial for convergence analysis. For practical purposes, we notice setting Up = 0 provides sufficient accurate solution without the expensive parallel transport operation.
In Algorithm 3
33

---

Algorithm 3 Tangent space conjugate gradient T, vn 7)
1: Set To = M, pe ro-
2: for 1= 0, ....
do
3: Compute Fi+1 H(ve]-
||r:11
4: a1+1 Prilo
5:
U+41 = Ut + OctiPt.
6:
Tt+1 ||r,1
a++1 (pe))
7: B[t1
1ir:11
8: PL-1 r141 BL+1Pt-
9: end for
Output: VT
Algorithm 4 Riemannian bilevel solver for min max optimization
1: Initialize 20 M2z yo My.
2: for = 0, ...,
do
3: vl 1/s.
4: for s = 0, ..., S 1 do
5:
Vk
u-+1 ExPyt (- Gug(Ex. v()).
6: end for
7: Update Tk+1 ExPoa ( No where yk-1 = 1/
8: end for
G Extensions: Riemannian mix-max and compositional optimization
The bilevel optimization considered in the paper I generalizes the two other widely studied problems. namely the max optimization and compositional optimization.
G.1
Min-max optimization on Riemannian manifolds
Riemannian min-max problems have gained increasing interest over the recent years 37 which takes the form of
2
67 25 56 35
min max fxx M yEM,
and can be seen as a special case of bilevel optimization problem where g(z, y) = ffx f(..)) v)- Because the problem is nonconvex in the order of minimization and maximization matters [73] 25 Nevertheless, under the assumption where f is geodesic strongly convex in y. the optimal solution x² satisfies G F(*")) 0, where F(x) = G.. (x. y (x)) due Gy (I. y (x:) Gyg(x, y²x)) 0.
Thus Algorithm reduces to alternating gradient descent ascent over Riemannian manifolds, as outlined in Algorithm
Here we adapt the convergence analysis to the min-max optimization setting. Given we no longer require second-order derivatives, we restate assumptions for functions f. below.
Assumption (1) Assumption 1 holds. (2) Function s(x, u), g(x, y) have L-Lipschitz Riemannian gradients. (3) Further, g(x,)) y) is j-geodesic strongly convex in y.
Under the min min-max setup and Assumption 6 we see GF(x) fzf(z, (x)) and thus the Lipschitz constant can be derived as Lr (k: + =) ((r)). ). Further we can directly apply Theorem for the Hessian inverse with Chiuv = L, which leads to the following convergence result.
Theorem 4. Under Assumption
choosing S 2 Ō(*9C).
=
20L, we have
K.... [[F < 50(x(+1)LL) and to reach e-stationary solution, we require gra- dient complexities as Gy O(1cc1) and G, Ö(
K
34

---

G.2 Compositional optimization on Riemannian manifolds
Compositional problems on Riemannian manifolds have been considered in 36 which requires to solve
min (((()), #
(20)
where 3 My - R and Φ .M2 My. It is worth noting that in both works 36 70), the inner function ó MM + He is vector-valued. In contrast, we consider general manifold-valuee function 6. Because the function ø can be potentially complex and may be stochastic, we follow [1 to reformulate 20) into bilevel optimization problem by letting
[(x,0) = (x)), s.t. (x) = arg min(g(r, v) = (6(2),)))) yeM,
As long as the squared Riemannian distance is geodesic strongly convex, the reformulation is equivalent to the original problem (20). As formally stated in Lemma | this is satisfied for non- positively curved space, like Euclidean space, hyperbolic manifold. SPD manifold with affine invariant metric. For positively curved space, the strong convexity guaranteed when restricting the domain relative to the curvature.
Lemma 9. Let U S .M has sectional curvature lower and upper bounded by n and n² respectively. Further U has diameter upper bounded by D, which satisfies D < Vx= if > 0. Then ler 8
VKF D taml taul Vx+D)) when k > O and consider Ç be the same curvalure conslant as in Lemma 3 Then function Hag(I, y) has Riemannian Hessian bounded within [8, () in spectrum.
when TT < 50ana and 8
Proof of Lemma 9
The proof follows from Lemma in 3 Consider an arbitrary curve y [0. +
M, and let fz) dd (Ir, p). for some p .M. From []. we know that H.f (t)
(p) and under the conditions, (1))(()) S (o (p), - (2)) S
(lb'(1)12 where we denote V as the covariant derivative. This immediately leads to
which completes the proof.
Thus, for positively curved manifold, if D < 2V++ we have > 0, which ensures geodesic strong convexity of the inner problem. As shown in Lemma 12 in B] G y) = 22xx,, )(x) and the Riemannian gradient descent on y lead
Expy; 17 Exp uf o(xx)),
which suggests v lies on a geodesic that connects of and d(xx). When S and when the lower-level function is vector-valued, the algorithm recovers the deterministic version of SCGD [66].
However, unlike in the Euclidean space, the Riemannian Hessian does not simplify to the identity operator, but rather the covariant derivative of inverse exponential map and the cross derivatives Glyg(z, v) # (Dó(z))!.
if + 0. (2) Function f(z, y) has 2 Riemannian gradients that are bounded by M and are L-Lipschitz. (3) Function has p-Lipschitz Riemannian Hessian and cross derivatives.
Assumption (1) Assumption holds and further D
We notice that for function y we only require second-order derivatives to be Lipschitz because the first-order Lipschitzness can be inferred from Lemma 9
Theorem 5. Under Assumption Theorem holds with C. 8.
To prove the convergence, we only need to show Lemma 4 holds. It can be readily proved from
Lemma 9 and Assumption that Lemma holds with I, S. Hence the convergence follows directly.
35

---

H Experimental details
H.1 Synthetic problem
We first verify the lower-level problem is geodesic strongly convex.
Proposition 4. For any A,B > O, function f(M) = (M, + (M-1, B) is u-geodesic strongly convex in U C $4 with >a,- + Ab. >- X=, where Au,-, are the minimum eigenvalue of
A, B and X are the bounds for maximum and minimum eigenvalues for M E U.
The inverse of Riemannian Hessian of function I [M) A) + (M-), B) is derived as, for any
symmetric U, H [4-1(M] M GM122 where G is the solution to the Lyapunov equation 1/3) + (M A/ + BMMM//BM--- 1/2) = M UM 1/2
G/ A/
Proof of Proposition We first derive the Euclidean gradient and Hessian
f(() = A - BM V = 'MM M BM
for any U UT. The Riemannian gradient and Hessian derived
G/(M) MAM - B
H/(M UM 'B UU (MV)
UM B + BM BU--UU + 2
(UAM UM BM + MAU BM-
1 (UAM+ MAU + UM- BM-
2
where we let {A}s (A + A 1/2. To show the function is geodesic strongly convex, it suffices
to show H. is positive definite, which is to show (H.(() U)M 2 U > for any U UTUT. To this end, we vectorize the Riemannian Hessian terms of as (M (M] =
(MA I+I& MA + BM- & &I+IS & BM ))vee(()) U), where @ denotes the Kronecker product.
Then, we have
(H](()[U). U)M
= =tr(M UUM''MM (MMMMI)(
2 V (M O M M (MA I I-18 MA BM @ I+ 18 BM (vec((U)
I
= 2 (A& + + & A @ M + M 6 M 'M'B'T veevve (U)
2 (Aa, - Ao Am Vve(UU) (M & Mvec)ve U) = U
where we let dn be the maximum/minimum eigenvalues of A and similarly for An Am.t.
The Hessian inverse can be derived subsequently. This completes the proof.
H.2 Computational time for each hypergradient estimator
This section report the average runtime in seconds (over 10 runs) for single evaluation of hypergradient using four different strategies (for Hypergradient estimation) for the synthetic problem (Section 1.1 Fi. The hyper-parameten are set be the same the main experiment. We in general automatic differentiation (AD) is the most efficient strategy. Nevertheless, according to Figure T(a), less accurate compared other strategies.
Table Comparison runtime for single computation of hypergradient.
HINV
CG
NS
AD
0.0154
0.1037
0.1030
0.0053
36

---

H.3 Hyperparameter selection
The selection of syper-parameters is performed reflect the best performance. The stepsize is selected from the range [le-3, 5e-3, 1e-2, 5e-2, le-1, 5e-1, 1] and for Neumann series is selected from [0.1, 0.5, 1.0, 1.5, 2.0] and number of inner iterations is selected from [S. 10, 30, 50, 100]. Figure d) shows the sensitivity of hypergradient error as we vary and the number of inner iterations.
H.4 Computational Complexity
This section lists the computational complexity for each task considered in the experiment section. In Table 4 we present an estimate of the per-iteration complexity of computing the gradient, Hessian/Jacobian-veerr vector products. We highlight that we only provide estimates of the complexities given that there may not exist closed form expressions for the gradient and second-orde derivatives.
Here, y R4 denote the size of validation set and training set respectively. For meta learning, denotes the number tasks and n denotes the number of samples for each task. For domain adaptation, m,n denote the number of samples for two domains, denotes the number of Sinkhorn iterations.
Table 4: Per-iteration complexity estimate each task
Hyper-rep (shallow)
Meta learning
Domain adaptation
x size
d x
d X r dxr O(mnd'r)
m X
y size
r(r + 1)/2
d
On + (
(smn)
Gg
O(nr*))
O(mnd*r)
Old mat nd")
JV
O(n,(r* + drr))
O(mnd'r)
Old + mdd + nd + smn)
HVg
O
O(mndr²)
Olds + mdd + +dd)
I Experiment Configurations
All the experiments are conducted on a single NVIDIA RTX 4060 GPU. All datasets used in the paper publicly available, which properly cited in the main paper. We include detailed setups for the experiments in the main paper as well as documented in code (provided as supplementary material).
J Broader Impact
This paper proposes new algorithms and are of theoretical in nature. We do not foresee any immediate negative societal impact that we feel obliged report.
37

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: Please refer to Section 3 for theoretical developments and Section for empirical studies.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer:
[Yes]
Justification:
We include discussions of limitations in the conclusion section, i.e., Section 5
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
38

---

We have included all the assumptions and complete proof in the paper. Please refer to Section and Section
Justification:
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: We have provided the implementation details for reproducing the experiments, along with the code in the supplementary material.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that appropriate to the research performed.
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
39

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
We have provided the code and data supplementary material ensure reproducibility.
Justification:
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nnis.cc
aublic/guudees/Codeebbbmsss
www.sssinnooiiillo for more details.
While encourage the release of code and data, understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
/ wwtp:/cccppplgguuuiiaeeeae
Stteee//CooSSSemsssS.Si.a www.ssssiicggyoo for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification:
We have included the details in the experiment section.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Does the report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: We have reported standard deviation in Figure 2(b) and (c). For Figure the results barely change as we vary the random seed. For Riemannian meta learning, it becomes time costly to report standard deviation.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
40

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not show in tables figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: We have included the detailed configuration in Appendix ! We report the time exccution in Appendix H.2
Guidelines:
The answer NA means that the paper does include experiments.
The paper should indicate type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (c.g., preliminary or failed experiments that didn't 't make into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics Ethics/https:
Answer: [Yes]
Justification:
The paper conform, in every respect, with the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (c.g.. if there special consid- eration due to laws or regulations in their jurisdiction..
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification:
We have included broader impact section, i.c., Section
Guidelines:
The answer NA means that there is no societal impact of the work performed.
41

---

If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), faimess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: This paper does not release data or models that have high risks for misuse.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but encourage authors to take this into account and make best. faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: In the experiment section, this paper has cited the papers when using exisiting assets, inluding code and data.
Guidelines:
The answer NA means that the paper does existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset used and, possible, include a URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
42

---

If assets are released, the license, copyright information, and terms of use in the thcode has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
package should be provided. For popular datasets, paperswi
latasets.paperssit
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information not available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: This paper does not introduce new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moode as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country institution) were obtained?
Answer: [NA]
Justification: This paper does not involve human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
43

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
44


Reviews:
[{"review_id": "c5NRujuRBd", "reviewer": "Reviewer_CDbu", "summary": "This paper addresses bilevel optimization problems where variables of both the lower and upper-level problems are constrained on Riemannian manifolds. To solve such problems, the authors propose a hypergradient descent algorithm under the key assumption that the lower function is geodesically strongly convex and establish its iteration complexity. To efficiently estimate the hypergradient, several methods with error bounds are also proposed. Additionally, the paper provides several interesting application examples of bilevel optimization on Riemannian manifolds and presents numerical results for these examples.", "questions": "1. Why are the performances of “hinv-20” and “hin-50” similar in Figures 1(a) and 2(a)?\n\n2. Why is “hinv” more efficient, as shown in Figure 1(b)?", "limitations": "Yes, the limitations are acknowledged.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- This paper generalizes the widely used hypergradient descent algorithm to the Riemannian manifold setting and provides a comprehensive analysis of the proposed algorithms.\n- The paper includes four approaches to estimating the inverse of the Riemannian Hessian and a stochastic version of the proposed algorithm.\n- The theoretical analysis is solid, with estimation error bounds and the iteration complexity of the proposed algorithm presented.\n- The numerical results are extensive, showcasing results on four interesting Riemannian bilevel optimization application examples.", "weaknesses": "- A closely related work is [42], and although the authors claimed key differences in Section 1, some points should be clarified. \n\n     - The various hypergradient estimators beyond conjugate gradient and Neumann might be standard. Do the other estimators show better performance? \n\n     - While the paper presents interesting applications (see instances given in Section 4.2), these problems can also be solved by the methods in [42]. The authors should provide a numerical comparison of the methods. \n\n     - Shall the authors include the complexities of the methods in [42] in Table 1?\n\n- The estimations of the inverse Hessian are also widely used in Euclidean space.\n\n- The assumption of geodesic strong convexity for the lower function is too strong for the Riemannian setting, making the application examples in the paper less typical. Additionally, it should be noted that the lower-level problem in Section 4.2 is constrained in Euclidean space.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for positive comments on recognizing our solid theoretical analysis and extensive numerical results. We also greatly appreciate  the constructive feedback and comments. \n\n\n**1. (W1) Clarification regarding [42]: (1) Do other estimators show better performance? (2) Provide numerical comparisons with [42]. (3) Include complexities of methods in [42].**\n\nThank you for raising these points. Below, we address each one individually.\n\n(1) There exists a accuracy-runtime trade-off for these estimators. The trade-off can be observed from Figure 1 where we compared various estimators in terms of convergence, hypergradient error and runtime. We have also additionally compared per-hypergradient runtime for these estimators in Table 3 (Appendix H.4). From the figures and table, we see hypergradient based on automatic differentiation requires the least runtime while suffers from a poor hypergradient approximation accuracy. Neumann series based estimators balance the trade-off while conjugate gradient requires the most runtime but has the lowest approximation error. Choosing a suitable estimator depends on specific applications.\n\n(2-3) We highlight that the methods presented in [42] are covered by the framework that we propose, i.e., Algorithm 1 in [42] coincides with RHGD-CG and Algorithm 2 in [42] is RSHGD-NS that we propose, except that they use exponential map for the update, while our framework works for general retraction. Thus, the methods in [42] can be considered as special cases in our framework where we have already included the complexities. \n\n*Please find  numerical comparisons with [42] in the one-page supplementary material.* In particular, we compare on the synthetic problem (Figure 1) and shallow and deep hyper-representation (Figure 2). We observe the use of retraction in our case results in more efficient update and thus reduces the runtime for convergence.\n\n\n\n\n\n**2. (W2) The estimations of the inverse Hessian are also widely used in Euclidean space.**\n\nWhile it is true that the estimation strategies for inverse Hessian have been explore in the Euclidean space, in more general manifold setting, certain more care is needed when using second-order derivatives. For example, $\\mathcal{G}^2_{xy}$ (used in Algorithm 1) is a second-order cross derivative relating two different manifolds (one for x and the other for y) and defining this is tricky. Additionally, Hess is now an *operator on the tangent space* instead of a matrix in the Euclidean space. Consequently, its inverse needs to be properly characterized and computed.\n\n\n\n\n\n\n**3. (W3) (1) Assumption of geodesic strong convexity is strong. (2) Lower-level problem in Section 4.2 is in Euclidean space.**\n\n(1) Analyzing bi-level problem with non-(geodesic)-strongly-convex lower-level problem is tricky (even in the Euclidean setting) because of non-unique solution of the lower-level problem. \nNon-unique solution would give rise to several challenging cases such hypergradient becoming ill-defined, invertibility of the Hessian, etc. However, it is possible to relax this assumption to geodesic convex lower objective, by adding a strongly convex regularizer term to the objective function.  \n\n(2) Since Euclidean space is a special case of Riemannian manifold, it is not unreasonable to consider problems where lower-level problem is in the Euclidean space. However, we would like to emphasize that Sections 4.1 and 4.4 explore lower-level problems constrained to the SPD manifold. \n\n\n**4. (Q1) Why performances of hinv-20 and hinv-50 similar?**\n\nIt should be noted that the numbers 20 and 50 refers to the number of inner iterations for solving the lower-level problem. The similar performance of hinv-20 and hinv-50 is because the lower-level problem for such synthetic problem is easy to solve and 20 iterations already provides a good close-to-optimal solutions. \n\n\n**5. (Q2) Why hinv more efficient as in Figure 1(b)?**\n\nIn this synthetic example, we can derive an analytic form of the Hessian inverse. In addition, the problem dimension is not large (d = 50) and thus computing the Hessian inverse directly is more efficient in this setup.", "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": "DmRJHHPqlM", "reviewer": "Reviewer_kqCt", "summary": "This paper introduces a novel approach for solving bilevel optimization problems where both upper and lower-level variables are constrained on Riemannian manifolds. The authors propose four hypergradient estimation strategies (HINV, CG, NC, AD), analyze their estimation errors with convergence and complexity analysis, and generalize to stochastic bilevel optimization and general retraction. The paper also shows many applications of the framework including hyper-representation over SPD matrices, Riemannian meta-learning, and unsupervised domain adaptation.", "questions": "1. What are the computational trade-offs when choosing between the different hypergradient estimation strategies?", "limitations": "The authors acknowledge that the framework relies on some strong mathematical properties that limit its broader applicability and may be improved in the future.", "rating": 6, "confidence": 3, "soundness": 4, "presentation": 2, "contribution": 3, "strengths": "1. The proposed framework extends bilevel optimization to Riemannian manifolds, which fills a gap in this field.\n2. The paper introduces multiple hypergradient estimation strategies and provides a thorough theoretical analysis.\n3. Practical Relevance: Demonstrates the applicability of the framework through several machine learning applications, including hyper-representation and meta-learning.", "weaknesses": "1. Assumptions: The framework relies on assumptions such as geodesic strong convexity and Lipschitz continuity, which may limit its applicability in some practical scenarios.\n2. The practical performance evaluation is limited to synthetic and small-scale problems, which may not fully represent its utility in large-scale real-world applications.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for acknowledging that our theoretical analysis is thorough and our work has practical relevance. We also appreciate the constructive feedback. \n\n**1. (W1) Assumption such as geodesic strong convexity and Lipschitz continuity.**\n\nThe assumption of (geodesic) strong convexity is common for analyzing bilevel optimization, which is also the case even in the Euclidean space [35,40,53]. This is because without strong convexity, the lower-level problem can be ill-defined, where there may exist many local solutions. Studying bilevel problem with non-(geodesic)-strongly-convex lower-level problem is even challenging in the Euclidean space. Nevertheless, it is possible to relax such assumption to (geodesic) convexity with an extra strongly convex regularizer, or to (Riemannian) PL condition where a global minimizer exists. \nOn the other hand, the assumption on Lipschitzness is also standard in the bilevel literature such as [35,40,46]. \n\n\n**2. (W2) Empirical evaluation is limited to synthetic and small-scale problems.**\n\nWe would like to emphasize that in addition to the synthetic and small-scale setting, we  evaluate our approach on real-world scenarios such as unsupervised domain adaptation on Office-Caltech dataset (Section 4.4) and meta-learning where we optimize a convolutional neural network (CNN) with orthogonal constraints on mini-Imagenet dataset (Section 4.3). In particular, the meta-learning example leads to a large-scale problem that involves training a 3-layer CNN,  with each layer kernel parameters constrained to Stiefel manifold. This results in a parameter set of $3 \\times 144 \\times 16$ parameters. In addition, we train over 20,000 tasks where each task contains 50 images of size $84 \\times 84$ (5-ways, 5-shots). We will make sure to include the details on experiment setup in the revised version. \n\n\n**3. (Q1) Computational trade-offs when choosing between the different hypergradient estimators.**\n\nWe have discussed computational trade-offs in Table 1 and Appendix H.4 (together with convergence plots in Figure 1), where we compare the runtime for evaluating hypergradient for different estimation strategies. In summary, CG is the most computational demanding but yields the best approximation while AD requires the least computational effort but gives a poor hypergradient estimation. NS balances the computation and estimation error.", "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": "IIQTb00Zyy", "reviewer": "Reviewer_Bnza", "summary": "This paper presents a framework for addressing bilevel optimization problems in which variables of both lower- and upper-level problems are constrained on Riemannian manifolds. It introduces multiple hypergradient estimation strategies on manifolds and investigates their estimation error. The paper includes convergence and complexity analyses for the proposed hypergradient descent algorithm on manifolds. Furthermore, it extends these advancements to stochastic bilevel optimization and the utilization of general retraction on the manifolds. The paper also demonstrates the practicality of the proposed framework through several applications.", "questions": "1. The notations $C_{hinv}$, $C_{cg}$, $c_{ns}$, and $C_{ad}$ in Theorem 1 are not utilized in the main body of the paper.\n\n2. Can the stochastic bilevel problem studied in this paper be extended to the general stochastic setting, i.e.,  $\\underset{{x\\in \\mathcal{M}_x}}{\\min} F(x) = \\mathbb{E}[f(x,y)]$ and $y^*(x) = \\underset{y\\in\\mathcal{M}_y}{\\text{argmin}}~ \\mathbb{E}[g(x,y)]$?\n\n3. For other questions, please see the weaknesses above.", "limitations": "NA", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper is well-written and is easy to read.\n\n2. This paper derives the intrinsic Riemannian hypergradient using the implicit function theorem, proposes four strategies for estimating the hypergradient, and provides estimation error bounds for all the proposed strategies.\n\n3. This paper presents the Riemannian hypergradient descent algorithm for tackling bilevel optimization problems on manifolds, providing convergence guarantees. Furthermore, it extends the framework to the stochastic setting and incorporates retraction mapping on the manifolds.", "weaknesses": "1. Firstly, my primary concern pertains to the novelty of this paper. The algorithm proposed in this paper appears to mirror the framework introduced in a prior work [42]. Therefore, it is crucial to delineate the similarities and disparities between this paper and the referenced study [42].\n\n2. The experimental results should include a comparison with the methods presented in [42].\n\n3. The analysis of retraction on manifolds lacks depth and should include the convergence results of three additional hypergradient estimators.\n\n4. Assumption 1 seems to be excessively strict, which is unusual in common bilevel optimization. Moreover, the notation $D$ in $d(x_1,x_2)\\le D$ may be confusing. Furthermore, Assumptions 2 and 3 could be inferred from Assumption 1, indicating that this paper may not necessitate Assumptions 2 and 3.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the general appreciation of our work as well as constructive comments that we address below.\n\n\n**1. (W1)  Similarities and disparities compared to [42].**\n\nWe would like to emphasize that [42] is a *concurrent* work and not a prior work. It became available (publicly) after we had completed the draft. Moreover, we would like to highlight several notable differences compared to [42], some of which have already been discussed in lines 48-61. \n\n(1) Our proposed framework is more general and subsumes the developments of [42]. In particular, Algorithms 1 and 2 in [42] are special cases of our proposed RHGD-CG and RSHGD-NS algorithms, respectively. While [42] use exponential map in their Algorithms 1 and 2, we study more general retraction, of which exponential map is one specific instance. It should be noted that our contributions is not limited to CG variants of proposed RHGD algorithm as we study other hypergradient estimation strategies as well. Overall, we provide hypergradient estimation error bound for a variety of estimation strategies and numerically compare them in terms of convergence, estimation error (Figure 1) and runtime (Table 3). \n\n(2) We present many applications that require to solve a Riemannian bilevel optimization problem. Many of them are of independent interest. \n\n(3) In Appendix G, we have shown how the proposed framework also leads to a convergence analysis for Riemannian min-max optimization and Riemannian compositional optimization problems, where they are special instances. \n\nTo summarize, our work provides a more general and practical framework than [42].\n\n\n\n**2. (W2) Experiment comparison with the methods presented in [42].**\n\nWe have now included numerical comparisons in the one-page supplementary material. Specifically, we compare the performance on the synthetic problem (Figure 1) and hyper-representation task (Figure 2). We see the use of retraction in our case improves the efficiency of the solvers, which is reflected in the reduced runtime to achieve convergence. \n\n\n\n**3. (W3) The analysis of retraction should include the results of three additional hypergradient estimators.**\n\nThank you for the suggestion. Below we include the convergence results for the three additional hypergradient estimators and will make sure to include them in the final revision.\n\n**Theorem.** (Extension of Theorem 3) Under Assumptions 1,2,3,5 and let $\\tilde L_F = 4 \\kappa_l c_R M + 5 \\bar c L_F$. Consider Algorithm 1 with the following settings. \n\n- HINV: Suppose $\\eta_x = \\Theta(1/\\tilde L_F), S \\geq \\tilde \\theta(\\kappa_l^2 \\zeta)$. Then $\\min_{k=0,..., K-1} \\| \\mathcal{G} F(x_k) \\|^2_{x_k} \\leq 16 \\tilde L_F \\Delta_0/K$. \n\n- CG: Suppose $\\eta_x = \\Theta(1/\\tilde \\Lambda), S \\geq \\tilde \\theta(\\kappa_l^2 \\zeta), T \\geq \\tilde \\Theta(\\sqrt{\\kappa_l})$, where $\\tilde \\Lambda = C_v^2 \\bar c + \\kappa_l^2 ( \\frac{5M^2 C_0^2 D^2}{\\mu} + \\bar c)$, then $\\min_{k=0,..., K-1} \\| \\mathcal{G} F(x_k) \\|^2_{x_k} \\leq 96 \\tilde \\Lambda (\\Delta_0 + \\| v_0^* \\|^2_{y^*(x_0)})/K$. \n\n- NS: Suppose $\\eta_x = \\Theta(1/\\tilde L_F), S \\geq \\tilde \\Theta(\\kappa_l^2 \\zeta)$. Then for any $\\epsilon>0$, $T \\geq\\tilde \\Theta(\\kappa_l \\log(1/\\epsilon))$, we have $\\min_{k=0,..., K-1} \\| \\mathcal{G} F(x_k) \\|^2_{x_k} \\leq 16 \\tilde L_F \\Delta_0/k + \\epsilon/2$. \n\n- AD: Suppose $\\eta_x = \\Theta(1/\\tilde L_F)$, $S \\geq \\tilde \\Theta(\\kappa_l^2 \\zeta \\log(1/\\epsilon))$ for any $\\epsilon > 0$, we have $\\min_{k=0,..., K-1} \\| \\mathcal{G} F(x_k) \\|^2_{x_k} \\leq 16 \\tilde L_F \\Delta_0/k + \\epsilon/2$. \n\n\nThe proof basically contains three parts, where we first show the convergence of inner iterations in terms of retraction, which is agnostic to the choice of hypergradient estimator. Then we bound the hypergradient estimation error where only AD based hypergradient needs different treatment when using retraction because it requires to differentiate through the retraction. For other hypergradient estimators, the bounds in Lemma 1 still hold given no retraction is used in their computation. Lastly we show the defined Lyapunov function is monotonically decreasing. \n\n\n\n\n\n**4.  (1) Assumption 1 is strict and unusual in bilevel optimization. (2) The notation $D$ in $d(x_1, x_2) \\leq D$ is confusing. (3) Assumptions 2 and 3 could be inferred from Assumption 1 and the necessity to include Assumption 2, 3.**\n\n(1) Assumption 1 is often unavoidable in the literature of Riemannian optimization, especially when dealing with geodesic (strongly) convex objectives [21, 37, 62, 63, 64]. This is required to ensure bounded curvature in order to apply the trigonometry distance bound (Lemma 3), which is crucial for achieving linear convergence of the inner iterations. (2) For notation clarity on $d(x_1, x_2) \\leq D$, we will change the notation in the revised manuscript. (3) It is correct that Assumption 2 and 3 can be inferred from Assumption 1. Nevertheless, we explicitly state such assumptions for clarity and for properly defining the Lipschitz constants. \n\n\n\n**5. (Q1) The notations $C_{hinv}$, etc in Theorem 1 are not used in the main text.**\n\nWe will move such notations to appendix.\n\n\n**6. (Q2) Extension to the general stochastic setting?**\n\nYes, the proposed framework can accommodate the general stochastic setting. To this end, we require to construct unbiased gradient estimate and properly control its variance. A more formal treatment of such general case is left for future exploration.", "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": "orOs6u4eRp", "reviewer": "Reviewer_4rfq", "summary": "The paper proposes an RHGD algorithm (Algo 1) to solve the bilevel optimization problems (line 11). \n- Thm 1 proves the convergence to a stationary point of $F$ when using different approximations of the hypergradient.\n- Thm 2 shows the convergence under stochastic setting.\n- Thm 3 proves the convergence with cheaper retraction.\nBoth synthetic numerical tests and applications to real problem are provided.", "questions": "1. Could you please clarify the synthetic problem in Sec 4.1 satisfies all the assumptions? I only find that Assumption 3 is verified. Assumption 2 is automatically satisfied due to the compactness of the Stiefel manifold. However, could you clarify how to verify Assumption 1?\n1. Maybe it is a dumb question, but in Sec. 3.1, why HINV and AD are different?", "limitations": "Small typo: In line 112 there are 2 negative signs.", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "The paper is well-written. The theoretical results are solid and the contributions are new to me. The algorithm is fully implementable. What I like is the convergence is not only established for the Hessian inverse but also its cheaper approximation. The synthetic experiment is carefully designed.", "weaknesses": "1. Assumption 1 is too strong: it requires the trajectory of $z$ to be in the unique geodesic neighborhood of $z^*$.\n1. Assumption 2 is slightly too strong, does this mean many simple functions like quadratic functions: $f(x, y)=|x|^2+|y|^2$ do not satisfy this assumption?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the positive feedback on our work, particularly acknowledging that our paper is well-written with solid theoretical results. We also greatly appreciate the constructive comments. \n\n\n**1. (W1) Assumption 1 is strong.**\n\nWe would like to highlight that Assumption 1 is often unavoidable for analyzing Riemannian optimization algorithms, especially when dealing with geodesic strongly convex functions, see for example [21, 37, 62, 63, 64] that also makes use of such an assumption. This assumption provides bound for the curvature, which is crucial for deploying a trigonometry distance bound (Lemma 3) and achieving linear convergence. This assumption can be satisfied for compact manifolds and by restricting the domain of interest. \n\n\n\n**2. (W2) Assumption 2 is strong.**\n\nAssumption 2 is not strong and can be inferred from Assumption 1 where we consider a compact domain. In this regard, the quadratic function you provided *satisfies* such an assumption when the domain is bounded. \n\n\n\n\n**3. (Q1) Verify the synthetic problem in Sec 4.1 satisfies Assumption 1.**\n\nThank you for the question. Assumption 1 on the bounded domain can be satisfied directly for the Stiefel manifold.  For the SPD manifold, we can take the maximal domain encompassing all the iterates considered. \n\n\n\n\n**4. (Q2) Why HINV and AD are different.**\n\nHINV computes the analytic expression of hypergradient using the Hessian inverse, which does not depend on the inner iterations. However, AD computes an estimate of hypergradient by differentiating through the inner iterations.", "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": "HIbYLzwmyA", "reviewer": "Reviewer_JBnb", "summary": "This paper studies Riemannian bilevel optimization, where variables of both lower and upper level problems are constrained on Riemannian manifolds. The authors propose several hypergradient based algorithms via Neumann series and automatic differentiation. \nConvergence analysis is provided for the proposed approaches. Experiments on synthetic problems, hyper-representation learning, meta-learning, unsupervised domain adaptation are provided to demonstrate the effectiveness of the proposed methods.", "questions": "See weakness for information.", "limitations": "See weakness for information.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1.\tThe studied topic is of interest to bilevel optimization research. As far as I know, this should be the first work to study this type of research.\n2.\tThe analysis is comprehensive. The algorithms cover AID and ITD based hypergradient approximations. \n3.\tQuite a few experiments are provided to support the theory.", "weaknesses": "1.\tMore motivation examples should be provided to validate the importance of bilevel optimization on manifolds. In the experiments, for example, in manifold meta-learning, it would be good to compare the method with meta-learning methods without the orthogonality constraint. \n2.\tAssumption 1 may be quite strong. It requires all iterates z_1,…, are bounded in a compact space. The analysis in previous standard bilevel optimization does not require this assumption. \n3.\tMain results contain complex constants and are a little bit hard to parse. Some simplifications may be helpful here.\n4.\tThe extension in 3.3 will be more interesting if Hessian-vector rather than Hessian-inverse is considered.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for appreciating the contributions of our work in terms of comprehensive analysis and supportive experiments, as well as providing the constructive comments. \n\n\n\n**1. (W1) More examples to validate the importance of bilevel optimization on manifolds.**\n\nThank you for the suggestion. We believe  we have provided sufficient motivating applications, showcasing the importance of bilevel optimization on manifolds. One such example is in Section 4.4 that motivates the use of bilevel optimization formulation for unsupervised domain adaptation (under the optimal transport framework). As far as we know, this is a novel application of bilevel optimization to domain adaptation, where we learn a whitening matrix $\\textbf{M}$ in the lower-level problem. Our numerical experiments show that employing Riemannian bilevel optimization yields more suitable adaptation by accounting for the metric structure (through whitening) in unsupervised domain adaptation.\n\n\n\n\n\n**2. (W2) Assumption 1 is strong and is non-standard in previous bilevel optimization.**\n\nStandard bilevel optimization in the Euclidean space does not require this assumption. However, for optimization on manifolds, especially when dealing with (geodesic) strongly convex objectives, such an assumption is standard and often unavoidable. See for example [21, 37, 62, 63, 64]. This assumption is to ensure bounded curvature, thus allowing to show fast convergence under strong convexity. \n\n\n\n**3. (W3) Simplification regarding complex constants.**\n\nThank you for the suggestion. We will simplify the notations in our revision.\n\n\n\n**4. (W4) Extend Hessian-inverse to Hessian-vector in Section 3.3.**\n\nWe agree that a full-spectrum analysis covering other hypergradient estimators (in particular the ones using Hessian-vector products) would be interesting. The analysis techniques would be similar to what we already did for Hessian inverse. In order to make the contents concise, we had decided to omit the analysis for other hypergradient estimation strategies. In the revised version, however, we would include a discussion in this regard.", "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:
We sincerely thank the reviewers for the reviews. We are particularly encouraged by the numerous positive comments on our work, including, **well-written** (Reviewer Bnza, Reviewer 4rfq),  **comprehensive analysis** (Reviewer CDbu, Reviewer kqCt, Reviewer JBnb),  **solid theoretical results** (Reviewer CDbu, Reviewer 4rfq), **interesting applications** (Reviewer CDbu), **extensive numerical results** (Reviewer CDbu, Reviewer kqCt, Reviewer JBnb). 

Additionally, we greatly appreciate your constructive feedback, which has significantly helped us improve the paper. In response to your comments, we have

- addressed each point individually, and 

- have attached a supplementary one-page PDF containing numerical comparisons with [42]. 


We hope we have addressed all your concerns and questions through our responses as well as the additional experimental results. We look forward to more interactive discussions during the discussion phase.

Meta-review: This paper proposes a framework for solving bilevel optimization problems with both lower and upper level variables constrained on Riemannian manifolds. It introduces hypergradient estimation strategies on manifolds, analyzes their estimation error, and provides convergence and complexity analysis for the proposed hypergradient descent algorithm. The framework is extended to stochastic bilevel optimization and general retraction. Additionally, the paper demonstrates the utility of the framework in several applications.

The paper makes notable contributions by developing techniques for hypergradient estimation on manifolds, providing theoretical analysis for the hypergradient descent algorithm and its stochastic variant, and demonstrating their efficiency through concrete applications.  The reviewers agreed that the paper was strong enough to merit publication, especially after considering the rebuttal. I recommend acceptance.
Decision: Accept (poster)