summaries: They authors show that for linear optimization on strongly convex domains the lazy subgradient method achieves a best-of-both-worlds behavior of a expected regret bound of O(\sqrt{N}) and O(\log(N)) for adversarial and stochastic adversaries, respectively. The best known prior bounds were O(\sqrt{log(N)}) adversarial regret with a more complicated algorithm (based on A,B-prod). On a technical side, the analysis for lazy subgradient is much more intricate and considers loss contributions in the direction of the comparator separately from other directions, and may be of interest in other online linear optimization settings.  Simulations are provided that convincingly argue for the proposed method in the proposed setting.

summaries: The authors study online linear optimization in the case of known strongly convex constraints sets. They show that the simple method of online subgradient descent, which is known to achieve $O(\sqrt{N})$ regret in the adversarial case, also achieves $O(\log(N))$ regret in the i.i.d. case. This result is more elegant than and slightly improves upon existing results in the strongly convex case.

summaries: The paper considers the problem of online linear optimization under a full information feedback setting on a strongly convex domain. In this setting, the authors show that the classical lazy subgradient method is optimal both in the adversarial and stochastic settings.  In particular, the authors show that for the same choice of hyper-parameters, the lazy sub-gradient method achieves O(\sqrt{N}) regret in the adversarial setting and O(log{N}) regret in the stochastic setting. To show this result, the authors have to do a more nuanced analysis of the algorithm (that relies on differential geometry methods) than existing works.

summaries: The authors study online linear optimization and show that the Lazy Online Subgradient Algorithm is *universal* on strongly convex domains: it simultaneously achieves the optimal $O(\sqrt T)$ regret when the costs are adversarial and the optimal $O(\log T)$ regret when the costs are stochastic. This result is a clear improvement over previous work of Huang et al, who described a more complex algorithm and didn't achieve the optimal rates. The proof is based on an intricate analysis, drawing from differential geometry, of the points selected by the algorithm, and in particular, how much of each step affects regret against the minimizer $x^*$ and how much is orthogonal to the direction $x_i - x^*$ and doesn't affect regret. This technique may be of interest to the broader ML community.

abstractText: We study Online Lazy Gradient Descent for optimisation on a strongly convex domain. The algorithm is known to achieve O( √ N) regret against adversarial opponents; here we show it is universal in the sense that it also achieves O(logN) expected regret against i.i.d opponents. This improves upon the more complex metaalgorithm of Huang et al [20] that only gets O( √ N logN) and O(logN) bounds. In addition we show that, unlike for the simplex, order bounds for pseudo-regret and expected regret are equivalent for strongly convex domains.

authors:
Daron Anderson
Douglas Leith

title: Optimality of the subgradient algorithm in the stochastic setting

venue: Arxiv E-prints,

title: Universal algorithms: Beyond the simplex

authors:
Andreas Argyriou
Charles A Micchelli
Massimiliano Pontil

title: On spectral learning

venue: Journal of Machine Learning Research,

authors:
Peter Auer
Chao-Kai Chiang

title: An algorithm with nearly optimal pseudo-regret for both stochastic and adversarial bandits

authors:
Aditya Bhaskara
Ashok Cutkosky
Ravi Kumar
Manish Purohit

title: Online Learning with Imperfect Hints

venue: Proceedings of the 37th International Conference on Machine Learning, Vienna,

authors:
Sébastien Bubeck
Nicolo Cesa-Bianchi
Sham M Kakade

title: Towards minimax policies for online linear optimization with bandit feedback

venue: In Conference on Learning Theory,

authors:
Sébastien Bubeck
Michael Cohen
Yuanzhi Li

title: Sparsity, variance and curvature in multiarmed bandits

venue: In Algorithmic Learning Theory,

authors:
Sébastien Bubeck
Aleksandrs Slivkins

title: The best of both worlds: stochastic and adversarial bandits

authors:
Nicolo Cesa-Bianchi
Yishay Mansour
Gilles Stoltz

title: Improved second-order bounds for prediction with expert advice

venue: Machine Learning,

authors:
Mohammed Dahleh
Munther A Dahleh
George Verghese

title: Lectures on dynamic systems and control

authors:
Steven De Rooij
Tim Van Erven
Peter D Grünwald
Wouter M Koolen

title: Follow the leader if you can, hedge if you must

venue: The Journal of Machine Learning Research,

authors:
Ofer Dekel
Nika Haghtalab
Patrick Jaillet

title: Online learning with a hint

venue: In Advances in Neural Information Processing Systems,

authors:
V.F. Demyanov
A.M. Rubinov

title: Approximate Methods in Optimization Problems. (Modern Analytic and Computational Methods in Science and Mathematics)

venue: IX + 256 S. New York

title: Rates of convergence for conditional gradient algorithms near singular and nonsingular extremals

venue: SIAM Journal on Control and Optimization,

title: General topology. Sigma Series in Pure and Applied Mathematics, Volume 6

venue: Heldermann Verlag Berlin,

authors:
Pierre Gaillard
Gilles Stoltz
Tim Van Erven

title: A second-order bound with excess losses

authors:
Dan Garber
Elad Hazan

title: Faster rates for the frank-wolfe method over strongly-convex sets

venue: In 32nd International Conference on Machine Learning,

title: Regret bounds for prediction problems

venue: In Proceedings of the twelfth annual conference on Computational learning theory,

authors:
Elad Hazan
Nimrod Megiddo

title: Online learning with prior knowledge

venue: In International Conference on Computational Learning Theory,

authors:
Ruitong Huang
Tor Lattimore
András György
Csaba Szepesvári

title: Following the leader and fast rates in online linear prediction: Curved constraint sets and other regularities

venue: J. Mach. Learn. Res.,

authors:
Sham M Kakade
Shai Shalev-Shwartz
Ambuj Tewari

title: Regularization techniques for learning with matrices

authors:
Ehsan Kazemi
Thomas Kerdreux
Liqiang Wang

title: Trace-norm adversarial examples

venue: arXiv preprint arXiv:2007.01855,

authors:
Thomas Kerdreux
Alexandre d’Aspremont
Sebastian Pokutta

title: Projection-Free Optimization on Uniformly Convex Sets

venue: arXiv e-prints, page arXiv:2004.11053,

authors:
Evgeny S Levitin
Boris T Polyak

title: Constrained minimization methods

venue: USSR Computational mathematics and mathematical physics,

title: Curvature of feasible sets in offline and online optimization

venue: arXiv preprint arXiv:2002.03213,

authors:
Jaouad Mourtada
Stéphane Gaïffas

title: On the optimality of the Hedge algorithm in the stochastic regime

title: Optimum bounds for the distributions of martingales in Banach spaces

venue: The Annals of Probability, 22(4):1679–1706,

title: Geodesics revisited. Chaotic Modeling and Simulation (CMSIM

authors:
Paat Rusmevichientong
John N Tsitsiklis

title: Linearly parameterized bandits

venue: Mathematics of Operations Research,

authors:
Amir Sani
Gergely Neu
Alessandro Lazaric

title: Exploiting easy data in online optimization

venue: Proceedings of the 27th International Conference on Neural Information Processing Systems - Volume

authors:
Yevgeny Seldin
Gábor Lugosi

title: An improved parametrization and analysis of the EXP3++ algorithm for stochastic and adversarial bandits

authors:
Yevgeny Seldin
Aleksandrs Slivkins

title: One practical algorithm for both stochastic and adversarial bandits

venue: In Proceedings of the 31st International Conference on Machine Learning,

title: Online learning and online convex optimization

venue: Found. Trends Mach. Learn.,

title: Strong convexity of sets and functions

venue: Journal of Mathematical Economics,

authors:
Chen-Yu Wei
Haipeng Luo

title: More adaptive algorithms for adversarial bandits

title: Differential geometry: Lecture notes

venue: Trinity College, Dublin,

authors:
Julian Zimmert
Tor Lattimore

title: Connections between mirror descent, thompson sampling and the information ratio

authors:
Julian Zimmert
Yevgeny Seldin

title: An optimal algorithm for stochastic and adversarial bandits

title: Online Convex Programming and Generalized Infinitesimal Gradient Ascent

heading: 1 Introduction

text:
Online linear optimisation is a repeated game with one opponent. On turn n we know the cost vectors a1, a2, . . . , an−1 ∈ Rd and select an action xn from the domain X ⊂ Rd. The opponent observes xn and chooses the next cost vector an ∈ Rd and we pay cost an · xn. Our goal is to select actions to make the regret ∑N i=1 ai · (xi − x∗) small for x∗ ∈ argmin{ ∑N i=1 ai · x : x ∈ X}.
If the cost vectors are selected by an adversarial opponent then a1, a2, . . . , an−1 give no information about an and selecting a good action seems like an impossible task. Despite this, there are algorithms that get O( √ N) regret against such opponents. The loss in performance against x∗ is O(1/ √ N) per turn which vanishes as N → ∞. The central algorithms in the field are Hedge [18, 34]; Prod [9]; and (sub)Gradient Descent, which has Greedy [41] and Lazy [34] variants. All three algorithms get O( √ N) regret. The first two are specialised to X the simplex. The latter works for any compact convex domain. The bulk of online optimisation literature is refinements of the central algorithms.
In many cases however the cost vectors are non adversarial. For example they might be determined by fluctuations in traffic, the weather or the stock market over a short period of time. The standard model of easy opponents is cost vectors drawn independently from some fixed distribution. In the easy setting the history a1, . . . , an provides information about the next cost vector. Indeed ãn+1 = a1+...+ann is an increasingly reliable sequence of estimates for a = E[an+1] and it seems reasonable to play xn+1 to minimise ãn+1 · x. This is called Follow-the-Leader (FTL). For example if X is the simplex FTL gives pseudo-regret E [ ∑N i=1 a·(xi−y∗)] ≤ O(1) for y∗∈argmin{a·x : x∈X}. Note this is
weaker than bounding the expected regret E [∑N i=1ai · (xi − x∗) ] . On the other hand FTL gives no guarantee against adversarial opponents, and we might not know in advance whether the opponent is i.i.d or adversarial. Hence we would like an algorithm that is universal in that it always gets O( √ N) regret but specialises to a much stronger bound against i.i.d opponents.
Recently Huang et al [20] solved this problem for smooth, strongly convex domains. For these domains they showed FTL gives O(logN) regret if ‖ ∑n i=1 ai‖ grows linearly.1 The result was generalised to nonsmooth domains by [22, 26]. Huang et al then combine FTL with the (A,B)-Prod meta-algorithm of [31] to ensure O( √ N logN) regret against adversarial opponents. In this work
1 Without smoothness they show the growth condition occurs with high probability against i.i.d opponents.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Algorithm 1: Online Lazy Gradient Descent Data: Cost vectors a1, a2, . . . ∈ Rd. Parameter η > 0. Domain X ⊂ Rd.
1 for n = 0, 1, . . . do 2 xn+1 =ΠX ( − η√
n
∑n i=1 ai ) 3 Receive an+1 and pay cost an+1 · xn+1
we improve and extend their results. In particular we show the meta-algorithm is unnecessary as the much simpler Online Lazy Gradient Descent already achieves O( √ N) and O(logN) bounds.
Terminology and Notation
Throughout ‖ · ‖ is the Euclidean norm. Given a vector subspace V ⊂ Rd and linear operator A : V → V write the operator norm as ‖A‖ = max{‖Av‖ : v ∈ V, ‖v‖ = 1}. It is known (see [10] Theorems 4.1 and 4.2) that ‖A‖ = maxj √ |λj | for λj the eigenvalues of ATA. Moreover if A is symmetric we have ‖A‖ = maxj |µj | for µj the eigenvalues of A. The set Z ⊂ Rd is called λ-strongly convex to mean for each x, y ∈ Z and α ∈ [0, 1] the ballB(αx+(1−α)y, r) is contained in Z for r = λ2α(1− α)‖x− y‖ 2. The pseudo-regret P[RN ] and expected regret E[RN ] are defined:
P[RN ] = E [ N∑ i=1 a·xi −min x∈X N∑ i=1 a·x ] E[RN ] = E [ N∑ i=1 ai ·xi −min x∈X N∑ i=1 ai ·x ] (1)
The Jensen Inequality implies P[RN ] ≤ E[RN ]. We write ΠX for the Euclidean projection onto the domain X . The Online Lazy Subgradient algorithm is given as Algorithm 1. Note this is different from the so-called Greedy Subgradient algorithm with starting point x1 = 0 and recursive update xn+1 =ΠX(xn−ηan/ √ n). Section 5 Figure 2(b) suggests Greedy Subgradient is not universal.
Summary of Contributions
Throughout we make the following assumptions on the domain and cost vectors. Both assumptions are nontrivial. For example Assumption 1 holds on the Euclidean and `p balls for p∈(1, 2] but not for the simplex or polytopes. Assumption 2 on the cost vectors is a typical model of an easy opponent and fails if the opponent is intelligent and adapts to our past moves. Assumption 1. The domain X⊂Rd is compact, m-strongly convex and contains the origin. Write D = max{‖x‖ :x, y∈X}. The boundaryM = ∂X is a (d−1)-dimensional C2 manifold. Namely each z∈M has a neighborhood U in Rd and C2 function F :U→R with nonzero gradient such that M∩U={x ∈ U : F (x)=0}. Such a function is called a coordinate patch at z.
Henceforth write N(x) for the outwards unit normal at x∈M. To represent N :M→Rd locally we can choose a coordinate patch F :U→R at x∗ and write∇N= ∇F‖∇F‖ . Since F is C
2 the matrix ∇N of partial derivatives exists. Assumption 2. The cost vectors a1, a2, . . . are i.i.d with expectation E[an] = a. For all n we have ‖an‖ ≤ L and ‖an−a‖ ≤ R and E‖an−a‖ ≤ δ. Let x∗ ∈ argmin{a ·x : x ∈ X} be the expected minimiser. The domain X being strongly convex implies x∗ is unique.
Our main result is Theorem 3 which says Online Lazy Gradient Descent achieves O ( L2 m‖a‖ logN ) expected-regret under Assumptions 1 and 2. This bound is essentially tight by [20] Theorem 9. The improved behaviour against i.i.d opponents is remarkable since the algorithm was designed with the adversarial case in mind, and predates the recent interest in universal algorithms. Our second contribution is Theorem 2 which states the gap between expected regret and pseudo-regret is O ( L2 m‖a‖ ) . Theorem 2 is of independent interest since it fails if the domain is the simplex.
From a technical standpoint, our Gradient Descent analysis goes beyond that of [20, 22, 26] for FTL. Their central idea is that the normal vector toM at xn ∈ argmin{(a1+ . . .+an−1) · x : x ∈X} points along An−1 =−a1+...+an−1n . Using ‖An−An−1‖=O(1/n) and how strong convexity says the unit normal changes quickly as we vary the basepoint, they prove ‖Θn−Θn+1‖=O(1/n) for the normal vectors Θn,Θn+1 at xn, xn+1. Summing the resulting series gives a logarithmic regret
bound. Unfortunately this method is insufficient to analyse Gradient Descent instead of FTL as it only gives ‖Θn−Θn+1‖=O(1/ √ n). and the resulting series give only an O( √ N) bound. Thus we
follow a new approach where we prove a high probability bound ‖N(xn)−N(x∗)‖ = O(1/ √ n) for the unit normals. For comparison existing works do not mention the expected minimiser x∗ at all. Then we use strong convexity to see the boundary surface locally looks like a quadratic with axis pointing along −a. Hence ‖xn−x∗‖ = O(1/ √ n) and only O(1/n) of the displacement is in the a-direction; and the larger orthogonal component does not contribute to pseudo-regret. Summing the series gives an O(logN) bound for pseudo-regret. To convert this into a bound for expected regret we use our Theorem 2 which says the gap between expected regret and pseudo-regret is finite.
Our differential geometry methods are more nuanced than existing works. Existing works consider only the magnitude ‖xn − xn+1‖ but we consider both the magnitude and direction of xn − x∗. The relevant object is∇N(z)v in Lemmas 1-5 and Proposition 1 which captures the rate of change of the unit normal as we perturb the basepoint in the v-direction.

heading: 2 Related Work

text:
The first strand of related work is universal algorithms. That is to say online algorithms with Õ( √ N) regret against adversarial opponents that specialise to get stronger bounds against i.i.d opponents. The original research focus was on the simplex and bandit feedback [4, 8, 32, 33, 37, 40]. In the full information setting the FlipFlop algorithm [11] interpolates between FTL and Hedge on the simplex to get O(1) and O( √ N logN) regret against i.i.d and adversarial opponents respectively. The (A,B)-Prod meta-algorithm [31] does something similar in much greater generality. Another result comes from [16]. They design a variant of the Prod algorithm [9] with different learning rates for each arm. As a corollary they get O(1) expected regret against i.i.d opponents with a unique best arm. Apart from [31] all the above algorithms apply only to the simplex.
However recent results suggest such complicated algorithms are unnecessary. In the full information setting it was recently proved that two of the central algorithms do much better against i.i.d opponents. Hedge gets O(log(d)L2∞/∆) pseudo-regret on the simplex [27]; and Online Lazy Gradient Descent gets O(L22/∆) pseudo-regret on polytopes
2 [1, 2].The latter is notable since it works for domains other than the simplex. Here d is the dimension; ∆ the suboptimality gap; and we assume ‖an‖ ≤ L2 and ‖an‖∞ ≤ L∞ where appropriate. The above bounds are in addition to the O( √ N) regret bounds of Gradient Descent and Hedge (see [34] Chapter 2) against adversarial opponents.
The second strand of related work is the study of online algorithms on strongly convex domains. Typically these are unit balls with respect to some norm. For example [39] consider the `p unit ball with p ∈ (1, 2] in the full-information setting. Another setting of interest [5, 12, 19] is when predictions are available for the direction of the next cost vector. In the bandit setting [7] consider `p balls; [6] consider the Euclidean ball; and [30] consider general strongly convex sets. They allow infinitely many arms but require the cost vectors be stochastically generated. In the offline setting [17] show the Frank-Wolfe method has O(1/n2) convergence rates − and hence finite regret − when the domain is strongly convex. [23] prove similar results for more general so-called uniformly convex sets. Under additional conditions the older works [13, 14, 25] show an exponential convergence rate.

heading: 3 Bound for Pseudo-Regret

text:
Here we show Online Lazy Gradient Descent gets O(logN) pseudo-regret against i.i.d cost vectors on a C2 strongly convex domain. This is in addition the O( √ N) regret bound for adversarial costs.
Theorem 1. Under Assumptions 1 and 2, Algorithm 1 with parameter η > 0 achieves pseudo-regret
P[RN ] ≤ 144
m‖a‖
(√ 8πDR
η + 4R2
) logN +O(1).
In particular for hyperparameter η = Θ(D/L) we have P[RN ] ≤ O ( L2 m‖a‖ logN ) .
By [34] Theorem 2.4 hyperparameter η=Θ(D/L) is optimal for both adversarial and i.i.d bounds. 2Huang et al [20] already proved an O(1) bound for polytopes with unique optimal vertex. However dependence on dimension and suboptimality gap is not obvious.
We begin with some geometric intuition for our proof strategy. To interpret Lemma 1 it is convenient to imagine the domain is translated to put z at the origin and the outward unit normal N(z) points along the d-th coordinate direction (henceforth called the height). Since N(z) has unit length, the left-hand-side of Lemma 1 is the height of the boundary at x. The lemma says the height changes quadratically with coefficient depending on the direction v as defined below.3
Definition 1. Given points x, y ∈ M a geodesic from x to y is a C2 path inM with endpoints x and y with ‖γ′(t)‖ = 1 for each t, and such that γ has minimal length among such paths. That length is called the geodesic distance d(x, y). The unit vector v ∈ Rd is called a direction from x to y to mean there is a geodesic γ : [0, d(x, y)]→M from x to y with γ′(0) = v.
In Lemma 6 we plug the expected minimiser z = x∗ into Lemma 1. The height then points along a and the left-hand-side of the bound becomes the pseudo-regret of x. We then (a) remove the dependence on direction and bound the coefficient by 1/m and (b) use strong convexity to convert the right-hand-side into the squared difference ‖N(x)−N(x∗)‖2. Finally we plug x = xn into Lemma 6 and use Lemma 7 to see ‖N(xn) − N(x∗)‖2 = O(1/n) with high probability. The idea is to write the unprojected action yn = − η√n ∑n i=1 ai as yn = −η √ na + η√
n
∑n i=1(a − ai). The first term points along − a ‖a‖ = N(x ∗) and has length Ω( √ n).
The second term is O(1) with high probability due to the i.i.d sum. Hence after normalisation the component of yn orthogonal to − a‖a‖ has size O(1/ √ n). Since the domain is bounded the
same holds if we instead normalise yn − w for any w ∈ X . In particular taking w = xn we get∥∥ yn−xn ‖yn−xn‖ − ( − a‖a‖ )∥∥ = O(1/√n). Since xn = ΠX(yn) we have yn−xn‖yn−xn‖ = N(xn). Since − a‖a‖ = N(x
∗) we get ‖N(xn)−N(x∗)‖2 = O(1/n) with high probability. Sum the series to get PN = O (∑N n=1 1 n ) = O(logN). Lemma 1. For each z∈M there exists a coordinate patch F : U → R such that for each x ∈ U ∩M and direction v from z to x we have
|N(z)(x− z)| ≤ |v T∇2F (z)v| ‖∇F (z)‖ d(x, z)2.
Lemma 1 is proved using Lemmas 2-5 and Proposition 1 which are in turn proved in Appendix II. The arguments there use only routine differential geometry, but we provide self-contained proofs, since the literature focuses on intrinsically defined surfaces, rather than surfaces embedded in Rd. Lemma 2. For each z ∈M and unit vector v tangent toM at z we have ‖∇N(z)v‖ ≥ m.
Lemma 3. For each coordinate patch F at z and unit tangent v toM at z we have v T∇2F (z)v ‖∇F (z)‖ ≥m.
Lemma 4 says the unit normal function has bounded second derivative, despite any complications caused by the function being only locally defined. Lemma 4. There exist constants M, r > 0 such that each z ∈ M has a coordinate patch F : B(z, r)→ R with ‖∇N(y)‖, ‖∇
2F (y)‖ ‖∇F (y)‖ ≤M for all y ∈ B(z, r).
Lemma 5. Let M be the constant from Lemma 4. For each geodesic γ : [0, d(x, y)]→M from x to y we have for d = d(x, y) and all t ≤ d the inequalities
(1) ‖γ′′(t)‖ ≤M (2) ‖γ(t)− γ(0)− tγ′(0)‖ ≤ M2 t 2 (3) ∣∣‖y − x‖ − d∣∣ ≤ M2 d2 Proposition 1. Each z ∈M has a neighborhood U in Rd such that for all x ∈ U ∩M and direction v from z to x we have
‖N(x)−N(z)‖ ≥ ‖∇N(z)v‖ 2 d(x, z)
Proof of Lemma 1. Let F,M, r be from Lemma 4. Since F is C2 Taylor’s theorem says there is a function ψ(t) defined on some neighborhood of the origin with limt→0 ψ(t) = 0 and
F (x) = F (z) +∇F (z)(x− z) + (x− z) T∇2F (z)(x− z)
2 + ψ(z − x)‖z − x‖2
3It is of course possible to define geodesics and directions using a non-Euclidean notion of distance. However in this paper we only consider Euclidean geodesics.
Hence there is r > 0 such that for all x ∈ B(z, r) we have |ψ(z − x)| ≤ m4 ‖∇F (z)‖ and so∥∥∥F (x)− F (z)−∇F (z)(x− z)− (x−z)T∇2F (z)(x−z)2 ∥∥∥ ≤ m4 ‖∇F (z)‖‖x− z‖2. Shrink r if necessary to ensure B(z, r) has geodesic diameter at most min{m/2M2, √ 2m/M3/2}. By assumption F (x) = F (z) = 0. Hence for d = d(x, z) the above gives∥∥∥∇F (z)(x− z)− (x−z)T∇2F (z)(x−z)2 ∥∥∥ ≤ m4 ‖∇F (z)‖‖x− z‖2 ≤ m4 ‖∇F (z)‖d2. Let v be a direction from z to x. For w = (x− z)− dv. The reverse triangle inequality gives∥∥∥ |∇F (z)(x− z)|− |(dv)T∇2F (z)(dv)|2 ∥∥∥≤ m4 ‖∇F (z)‖d2+|wT∇2F (z)(dv)|+ |wT∇2F (z)w|2 . (2) Lemma 5 says ‖w‖ ≤Md2/2 and Lemma 4 says ‖∇
2F (z)‖ ‖∇F (z)‖ ≤M . Hence we get
|wT∇2F (z)(dv)| ≤ ‖w‖‖∇2F (z)‖‖dv‖ ≤ M2 ‖∇2F (z)‖ ‖∇F (z)‖ ‖∇F (z)‖d 3
≤ M 2 2 ‖∇F (z)‖d 3 ≤ ‖∇F (z)‖4 md 2.
where the last inequality uses the assumption d ≤ m/2M2. Similarly d ≤ √
2m/M3/2 gives |wT∇2F (z)w| ≤ ‖∇F (z)‖2 md
2. Hence (2) becomes∥∥∥ |∇F (z)(x− z)| − |vT∇2F (z)v|2 d2∥∥∥ ≤ m4 ‖∇F (z)‖d2 + m4 ‖∇F (z)‖d2 = m2 ‖∇F (z)‖d2. and so |∇F (z)(x − z)| ≤ ( |vT∇2F (z)v|
2 + m 2 ‖∇F (z)‖
) d2. Since the unit normal is
N(z) = ∇F (z)‖∇F (z)‖ we can divide both sides to get N(z)(x− z)| ≤ 1 2 ( |vT∇2F (z)v| ‖∇F (z)‖ +m ) d2. Lemma 3 says m ≤ |v T∇2F (z)v| ‖∇F (z)‖ . This completes the proof.
Lemma 6. The minimiser x∗ has a neighborhood U in Rd such that each x ∈ U ∩M has
a · (x− x∗) ≤ 4‖a‖ m ‖N(x)−N(x∗)‖2.
Proof. Let F : U → R be a coordinate patch at z = x∗ that meets the conditions in Lemma 1 and Proposition 1. Since x∗ = argmin{a · x : x ∈ X} the outwards unit normal is N(x∗)=− a‖a‖ and the lemma gives
|a · (x− x∗)| ≤ ‖a‖ |v T∇2F (x∗)v| ‖∇F (x∗)‖ d(x, x∗)2.
For all x ∈ U and direction v from x∗ to x the proposition gives d(x, x∗)≤2‖N(x)−N(x ∗)‖
‖∇N(x∗)v‖ . Combine the two bounds to get
|a · (x− x∗)| ≤ 4‖a‖ ( |vT∇2F (x∗)v| ‖∇F (x∗)‖
1
‖∇N(x∗)v‖2
) ‖N(x)−N(x∗)‖2.
To complete the proof Lemma 3 says |v T∇2F (x∗)v| ‖∇F (x∗)‖ = |v T∇N(x∗)v|≤‖v‖‖∇N(x∗)v‖≤‖N(x∗)v‖. Hence the middle factor above is at most 1‖N(x∗)v‖ which is at most 1/m by Lemma 2.
The proofs of Lemmas 7 and 8 in Appendix III follow easily from the vector martingale theorems of Pinelis [28]. The use of these theorems − stated only for bounded martingales − is the major obstruction to generalising our main results to unbounded but subgaussian cost vectors. Lemma 7. For each n≥16D2/η2‖a‖2. The Online Lazy Gradient Descent actions x1, x2... give
E‖N(x∗)−N(xn+1)‖2 ≤ 1
n
36
‖a‖2
(√ 8πDR
η + 4R2
) + 8 exp ( − ‖a‖ 2
32R2 n
) .
Lemma 8. For each neighborhoodU of x∗ in Rd the Online Lazy Gradient Descent actions x1, x2, . . . give ∑∞ i=1 P (xi /∈U) <∞.
Proof of Theorem 1. It is enough to prove a tail bound E [∑N
i>n0 a·(xi − x∗)
] ≤ O(logN) for some
fixed n0. To use Lemma 7 let n0 =b16D2/η2‖a‖2c.4 Let U be a neighborhood of x∗ that meets the conditions of Lemma 6. Since each a · (xn+1 − x∗) ≤ 2LD we see
E|a · (xn+1 − x∗)| ≤ 2LD · P (xn+1 /∈ U) + ∫ xn+1∈U |a · (xn+1 − x∗)|dP.
Lemma 6 says the integral is at most∫ xn+1∈U 4‖a‖ m ‖N(xn+1)−N(x∗)‖2dP ≤ 4‖a‖ m E|N(xn+1)−N(x∗)‖2.
and we see E|a · (xn+1 − x∗)| is at most
2LD · P (xn+1 /∈ U) + 4‖a‖ m E‖N(xn+1)−N(x∗)‖2.
Lemma 8 says the series ∑∞ i=1P (xi /∈ U) converges. Lemma 7 gives constants k ≥ 0 and K =
36 ‖a‖2 (√ 8πDR η + 4R 2 ) such that E‖N(xn+1)−N(x∗)‖2 ≤ K/n+8e−kn. Hence we get
E
[ N∑
i>n0
a·(xi − x∗) ] = O(1)+
4‖a‖ m N∑ i=n0 ( K/i+ 8e−ki ) = O(1)+ 4‖a‖ m ∫ N n0 ( K/x+ 8e−kx ) dx
= O(1) + 4‖a‖ m K logN = O(1) + 4‖a‖ m K logN = 144 m‖a‖
(√ 8πDR
η + 4R2
) +O(1).
4 Equivalence of E[RN ] and P[RN ]
Our Theorem 1 gives an O(logN) bound for pseudo-regret. This seems weaker than the O(logN) expected regret bound of [20]. Here we show the two types of order bounds are in fact equivalent for strongly convex domains. Theorem 2. Under Assumptions 1 and 2, Algorithm 1 with parameter η > 0 gives
E[RN ]− P[RN ] ≤ 31R2
m‖a‖ = O
( L2
m‖a‖
) .
Theorem 2 holds even for non smooth domains. We need only the following result of Vial [35]. Lemma 9. (Vial Theorem 1(v)) Suppose Z is m-strongly convex. For all x, y ∈ ∂Z and outward unit normals p, q at x, y respectively we have ‖x− y‖ ≤ 1m‖p− q‖.
To prove Theorem 2 we first put the left-hand-side in a convenient form. Lemma 10. Let y∗ ∈ argmin {∑N i=1 ai · x : x ∈ X } . Under Assumption 2 we have
E[RN ]− P[RN ] ≤ √ E‖x∗ − y∗‖2 √ E ∥∥∑N i=1(ai − a) ∥∥2.
Proof. By definition E[RN ]− P[RN ] equals E [∑N i=1 ai · (xi − y∗) ] − E [∑N i=1 a · (xi − x∗) ] =
E [∑N i=1(ai−a)·xi ] +E [∑N i=1 ai ·(x∗−y∗) ] . To see the first term is zero recall each xi is a function of a1, . . . , ai−1. Since a1, a2, . . . are independent we see each xi is independent of ai and independent of ai − a. Thus the expectation distributes over the product. Since E[ai − a] = 0 the first term vanishes. Write the second term as E [∑N i=1(ai−a) · (x∗−y∗) ] −E [∑N i=1 a · (y∗−x∗) ] . Since x∗
minimises a·x the second term is positive and we get E[RN ]−P[RN ] ≤ E [∑N i=1(ai−a)·(x∗−y∗) ] . To complete the proof apply Cauchy-Schwarz to the L2 norm. 4According to the adversarial bound the regret before turn n0 is O( √ n0)≤O(D/η‖a‖).
To prove Theorem 2 we need the following three results from Appendix III. Like in Section 3 the first two follow from the vector martingale theorems of Pinelis [28]. The third result is well known but we were unable to find a sufficiently general proof in the literature. Proposition 2. Suppose X1, X2, . . . ∈ Rd are independent random variables with each E[Xn] = 0 and ‖Xn‖ ≤ R. For each t > 0 we have
P (∥∥∥ 1 n n∑ i=1 Xi ∥∥∥ > t) ≤ 2 exp(− t2 2R2 n ) .
Lemma 11. For y∗=argmin {∑N i=1ai ·x : x∈X } the Online Lazy Gradient Descent actions give
E‖N(x∗)−N(y∗)‖2 ≤ 144R 2
‖a‖2N + 4 exp
( −‖a‖ 2
8R2 N
) .
Lemma 12. For X a nonnegative random variable we have E[X] = ∫∞ 0 P (X > x)dx.
Proof of Theorem 2. Proposition 2 gives P ( ‖ ∑N i=1(ai − a)‖ > t ) ≤ 2 exp ( − t 2 2NR2 ) for all
t ≥ 0. Replace t with t2 to see P (∥∥∑N i=1(ai−a) ∥∥2 > t) ≤ 2 exp (− t2NR2 ) . Lemma 12 says
E‖ ∑N i=1(ai − a)‖2 ≤ ∫∞ 0 P ( ‖ ∑N i=1(ai − a)‖2 > t ) dt ≤ 4NR2. Hence Lemma 10 gives
E[RN ]− P[RN ] ≤ 2R √ N √ E‖x∗ − y∗‖2.
For the second factor Lemma 9 gives √ E‖x∗ − y∗‖2 ≤ 1m √ E‖N(x∗)−N(y∗)‖2. Lemma 11 gives E‖N(x∗) − N(y∗)‖2 ≤ K/N + 4e−kN for the constants K = 144R2/‖a‖2 and k = ‖a‖2/8R2. Hence √ E‖x∗ − y∗‖2 ≤ 1m √ K/N + 4e−kN ≤ 1m √ K N + 2 me −(k/2)N and so
E[RN ]− P[RN ] ≤ 2R √ K
m +
4R
m
√ Ne−(k/2)N = 24R2
m‖a‖ +
4R
m
√ Ne−(k/2)N . (3)
To bound the last term consider the function √ xe−βx for β ≥ 0. It is straightforward to see the function is increasing and then decreasing. By differentiating we see the maximiser is x = 1/2β and the function is at most 1/ √ 2eβ. For β = k/2 we see √ Ne−(k/2)N ≤ 1/ √ ke = 2 √ 2R/ √ e‖a‖.
Thus the right-hand-side of (3) is at most (
24 + 8 √ 2√ e
) R2
m‖a‖ ≤ 31R2 m‖a‖ .
Finally we can combine Theorems 1 and 2: Theorem 3. Under Assumptions 1 and 2, Algorithm 1 with parameter η > 0 gives
E[RN ]≤ 144
m‖a‖
(√ 8πDR
η + 4R2
) logN +O(1).

heading: 5 Numerical Simulations

text:
Here we present numerical simulations for Online Lazy Gradient Descent, using a range of strongly convex domains and i.i.d opponents. The higher dimensional simulations ran on the order of minutes, due to use of an all-purpose Python package to compute minimisers. Performance time can be improved by using the Franke Wolfe method of [17] to simplify the minimisation problems.
`p-Norm Balls. Here we examine performance for domains X= { x∈Rd : (∑d i=1 |x(i)|p )1/p≤1}. For p ∈ (1, 2 ] Corollary 1 of [17] says X is (p− 1)d 1 2− 1 p strongly convex.
For simplicity we use stepsize η = 1. When searching for worst-case performance it is enough, by symmetry of the domain, to only consider E[an] = a with nonnegative nondecreasing components. We consider a = b/‖b‖ for b(i) = (i/d)r and a range of parameters r ≥ 0. For example three degenerate cases are a = 1√
d (1, . . . , 1) for r = 0; a =
√ 2
d(d+1) (1, 2, . . . , d) for r = 1 and
a = (0, . . . , 0, 1) for r →∞. We use cost vectors an = a+ µn for two types of noise: (1) µn(j) =
1√ d (Bn1 , . . . , B n d ) (2) µn(j) = (0, . . . , 0, B n d ).
Here Bnj are independent and take values ±1 with equal probability. In both cases ‖µn‖ = 1. Noise type (1) perturbs all components equally and individual perturbations are small for large d. Type (2) perturbs only the largest with size 1. We run several instances (samples) of each online problem and plot the average pseudo-regret over instances.
Figure 1(a) shows that for noise type (1) performance degrades as r shrinks. This is consistent with the picture for d = 2 and r = 0 where the minimiser (2−1/p, 2−1/p) also minimises the curvature, and how Theorem 1 gives a worse bound for small curvature. For noise type (2) Figure 1(b) shows larger r gives worse performance. This can be attributed to how r is approximately (0, . . . , 0, 1) and the noise is focused on the final component rather than spreading over all components evenly.
For smaller of d, p the pattern is the same but takes longer to emerge. In Figure 3, Appendix I, we consider the worst cases in Figure 1 and vary the dimension. For r = 0 and noise type (1) the performance improves slightly with large d. This is perhaps due to the higher dimension having a regularising effect on the noise. For r = 8 and noise type (2) the performance deteriorates with large d. This is expected since the curvature (p− 1)d 1 2− 1 p shrinks in high dimensions. In Figure 4, Appendix I, we vary the parameter p. In both worst cases the performance deteriorates as p → 1 though it also takes longer for the worse performance to emerge.
For all cases the growth rates are orders of magnitude below those suggested by Theorem 1. For example in the setup of Figure 1 we have p = 1.5, d = 32 and the convexity parameter is m = (p− 1)d 1 2− 1 p = 0.5 · 32− 14 ' 0.21. This makes the coefficient in Theorem 1 approximately 6180.
Comparison with Greedy Gradient Descent. Theorem 3 says Online Lazy Gradient Descent gives logarithmic expected regret on a strongly convex domain against i.i.d cost vectors. One can ask does Online Greedy Gradient Descent with update xn+1 =ΠX ( xn−ηan/ √ n )
give the same bound.
Figure 2(b) suggests otherwise. The worst-case regret of Greedy Gradient Descent seems to be
Ω(
√
n). The coefficients D1, C1 in Figure 2(b) are chosen to match the regret on turns 1000 and 500,
and the Θ( √ n) fit is very close over the entire range. The coefficient C2 is chosen to match on turns 0 and 1000 and is a worse fit.
Schatten Matrix Norms. Schatten matrix norms are well-studied as regularisers in the context of machine learning. For example see [3, 21] and the citing works.
Definition 2. Given a real n×m matrix A there is a unique decomposition A= UΣV such that V,U are real orthogonal, U is n×n, V is m×m, and Σ is n×m with real diagonal entries σ1≥σ2≥ . . .≥σmin(n,m) and all other entries zero. The σi are called the singular values of A. For any p≥0 define the Schatten p-norm ‖A‖S(p) =‖σ‖p as the p-norm of the vector of singular values.
For simplicity we only consider square matrices. Corollary 2 of [17] says the Schatten unit balls Bp(1) = {X ∈ Rd×d : ‖X‖S(p) ≤ 1} are (p− 1)d 1 2− 1 p strongly convex. When searching for worst case performance the following lemma says it is enough to only consider when E[an] is a diagonal matrix. See Appendix I for proof.
(a) Lazy Gradient Descent on the Schatten unit ball (b) Greedy Gradient Descent on the Euclidean unit ball
Figure 2
Lemma 13. Suppose a1, a2, . . . ∈ Rd×d are i.i.d cost vectors with E[an] = a. There is another set of i.i.d cost vectors b1, b2, . . . such that E[bn] = b is diagonal and ‖an − a‖ = ‖bn − b‖; ‖an‖ = ‖bn‖; ‖a‖ = ‖b‖ and the expected regret of playing Online Lazy Gradient Descent (Algorithm 1) on Bp(1) against a1, a2, . . . with x0 = 0 is the same as playing against b1, b2, . . . with x0 = 0.
Likewise by symmetry we can assume the diagonal entries are nonnegative and nondecreasing. For simplicity we use stepsize η = 1. We consider a = b/‖b‖ for b = diag ( (1/d)r, (2/d)r, . . . , 1 ) and a range of parameters r ≥ 0. We consider cost vectors an = a+ µn and noise µn = νn/‖νn‖ for three types of νn. Here Bnij are independent and take value ±1 with equal probability.
(1) νnij = B n ij for all i, j (2) νnij = B n ij for all i and some uniformly chosen j; and ν n ik = 0 for all k 6= j. (3) νnij = B n ij for some uniformly chosen i, j; and ν n lk = 0 for all (l, k) 6= (i, j).
Compared to the examples for `p balls, fewer samples are needed to obtain smooth curves. For noise type (1) Figure 2(a) shows the worst performance occurs for r = 0. The different types of noise make minimal difference to performance. See Figure 5, Appendix I. Thus when varying p and d we only consider noise type (1). In Figure 6, Appendix I, we see large d leads to slightly better performance, and small p leads to worse performance. Though the infuence of p and d is minimal compared to that in Figures 5 and 2 for `p-balls.
For Schatten balls diagonal noise matrices seem to be the most difficult to deal with. In that case all an are also diagonal and the Schatten norm reduces to the p-norm of the diagonal entries.

heading: 6 Conclusion

text:
Online Lazy Gradient descent on a strongly convex domain is universal in the sense that it achieves O(logN) regret against i.i.d opponents. This is in addition to the well-known O( √ N) regret bound against adversarial opponents. Moreover for strongly convex domains, order bounds for pseudo-regret and expected regret are equivalent. For various types of noise on `p and Schatten unit balls, the growth coefficients are orders of magnitude smaller than those suggested by our analysis.
Disclosure of Funding. This work was supported by Science Foundation Ireland grant 16/IA/4610.
Broader Impact. This is theoretical work on special cases of online optimisation, with no particular application in mind. We do not forsee the work having any direct impact on society. The only ethical concerns are those related to the use of online optimisation in general.

title: Online Lazy Gradient Descent is Universal on Strongly Convex Domains


summaries: The authors consider the Riemannian optimization with Bures-Wasserstein (BW) geometry for symmetric positive definite (SPD) matrix manifold. The authors compare the proposed approach with its counterpart with the popular Affine-Invariant (AI) geometry.  The authors illustrate that the proposed approach (with BW geometry) is more suitable than its counterpart with AI geometry when SPD matrices are ill-conditioned.

summaries: This paper studies the Bures-Wasserstein (BW) geometry of symmetric positive-definite (SPD) matrices and compares it to the commonly-used affine-invariant (AI) geometry in the context of Riemannian optimization. It is shown theoretically that the BW geometry has several advantages. Particularly, it is shown that (i) the BW metric is linearly dependent on the SPD matrix, (ii) the condition number of the Riemannian Hessian of the BW metric is smaller than that of the AI metric, and (iii) the curvature constant of the BW metric is smaller than that of the AI metric. Based on these properties, it is shown that the convergence of Riemannian steepest descent and Riemannian trust-region based on the BW metric is faster than based on the AI metric. The theoretical results are demonstrated on simulations of six problems: weighted least squares, Lyapunov equations, trace regression, metric learning, log-det maximization, and Gaussian mixture model.

summaries: This paper provides several results about the Bures-Wasserstein  (BW) geometry on the positive definite matrices. A thorough comparison with the Affine Invariant (AI) metric is made. The main focus of the paper is on optimization: is it faster to optimize a function under the AI or under the BW metric? The authors argue that when the solution of the problem $X^*$ is ill-conditioned, then algorithms derived under BW metric have faster convergence theoretical convergence. This is due to the conditioning of the Riemannian Hessian at the optimum, which is better for the BW metric when $X^*$ has poor conditioning. The authors also demonstrate that a variety of functions that are convex for the AI geometry are also convex for the BW geometry. Experiments on an array of synthetic problems validate that both Riemannian gradient descent and Riemannian trust region methods are faster in the BW geometry.

summaries: This paper investigates the advantages of the Bures-Wasserstein (BW) metric over the widely used Affine Invariant (AI) one, for the manifold of symmetric positive definite (SPD) matrices. Rewriting the two metrics in a compatible manner reveals that for a wide variety of objectives, the BW metric is a linear function of the input matrix whereas the AI has a quadratic dependence. One exception is the log-det function which is included in the analysis. The paper experimentally validates that for ill-conditioned problems BW is a better choice compared to AI.

abstractText: In this paper, we comparatively analyze the Bures-Wasserstein (BW) geometry with the popular Affine-Invariant (AI) geometry for Riemannian optimization on the symmetric positive definite (SPD) matrix manifold. Our study begins with an observation that the BW metric has a linear dependence on SPD matrices in contrast to the quadratic dependence of the AI metric. We build on this to show that the BW metric is a more suitable and robust choice for several Riemannian optimization problems over ill-conditioned SPD matrices. We show that the BW geometry has a non-negative curvature, which further improves convergence rates of algorithms over the non-positively curved AI geometry. Finally, we verify that several popular cost functions, which are known to be geodesic convex under the AI geometry, are also geodesic convex under the BW geometry. Extensive experiments on various applications support our findings.

authors:
P-A Absil
Christopher G Baker
Kyle A Gallivan

title: Trust-region methods on Riemannian manifolds

venue: Foundations of Computational Mathematics,

authors:
P-A Absil
Robert Mahony
Rodolphe Sepulchre

title: Optimization algorithms on matrix manifolds

authors:
Naman Agarwal
Nicolas Boumal
Brian Bullins
Coralia Cartis

title: Adaptive regularization with cubics on manifolds

venue: Mathematical Programming,

authors:
Awad H Al-Mohy
Nicholas J Higham

title: Computing the Fréchet derivative of the matrix exponential, with an application to condition number estimation

venue: SIAM Journal on Matrix Analysis and Applications,

authors:
Jesús Alcalá-Fdez
Luciano Sánchez
Salvador Garcia
Maria Jose del Jesus
Sebastian Ventura
Josep Maria Garrell
José Otero
Cristóbal Romero
Jaume Bacardit
Victor M Rivas

title: Keel: a software tool to assess evolutionary algorithms for data mining problems

venue: Soft Computing,

authors:
Luigi Ambrosio
Nicola Gigli
Giuseppe Savaré

title: Gradient flows: in metric spaces and in the space of probability measures

venue: Springer Science & Business Media,

authors:
Vincent Arsigny
Pierre Fillard
Xavier Pennec
Nicholas Ayache

title: Geometric means in a novel vector space structure on symmetric positive-definite matrices

venue: SIAM journal on matrix analysis and applications,

authors:
David Arthur
Sergei Vassilvitskii

title: k-means++: The advantages of careful seeding

venue: Technical report,

authors:
Christian Berg
Jens Peter Reus Christensen
Paul Ressel

title: Harmonic analysis on semigroups: theory of positive definite and related functions, volume

title: Optimisation on Manifolds in Julia, 2019. https://github.com/ kellertuer/Manopt.jl

title: Positive definite matrices

venue: Princeton university press,

authors:
Rajendra Bhatia
Tanvi Jain
Yongdo Lim

title: On the Bures-Wasserstein distance between positive definite matrices

venue: Expositiones Mathematicae,

authors:
Mukul Bhutani
Pratik Jawanpuria
Hiroyuki Kasai
Bamdev Mishra

title: Low-rank geometric mean metric learning

venue: ICML workshop on Geometry in Machine Learning (GiMLi),

title: Stochastic gradient descent on Riemannian manifolds

venue: IEEE Transactions on Automatic Control,

title: Riemannian trust regions with finite-difference Hessian approximations are globally convergent

venue: In International Conference on Geometric Science of Information,

title: An introduction to optimization on smooth manifolds

venue: Available online,

authors:
Nicolas Boumal
Bamdev Mishra
P-A Absil
Rodolphe Sepulchre

title: Manopt, a Matlab toolbox for optimization on manifolds

venue: The Journal of Machine Learning Research,

authors:
Daniel A Brooks
Olivier Schwander
Frédéric Barbaresco
Jean-Yves Schneider
Matthieu Cord

title: Exploring complex time-series representations for Riemannian machine learning of radar data

venue: In IEEE International Conference on Acoustics, Speech and Signal Processing,

authors:
T Tony Cai
Anru Zhang

title: ROP: Matrix recovery via rank-one projections

venue: Annals of Statistics,

authors:
Zeineb Chebbi
Maher Moakher

title: Means of Hermitian positive-definite matrices based on the log-determinant α-divergence function

venue: Linear Algebra and its Applications,

authors:
Calin Cruceru
Gary Becigneul
Octavian-Eugen Ganea

title: Computationally tractable Riemannian manifolds for graph embeddings

venue: In AAAI Conference on Artificial Intelligence,

authors:
Henri De Plaen
Michaël Fanuel
Johan AK Suykens

title: Wasserstein exponential kernels

venue: In International Joint Conference on Neural Networks,

authors:
Ian L. Dryden
Alexey Koloydenko
Diwei Zhou

title: Non-Euclidean statistics for covariance matrices, with applications to diffusion tensor imaging

venue: The Annals of Applied Statistics,

authors:
Matthieu Guillaumin
Jakob Verbeek
Cordelia Schmid

title: Is that you? Metric learning approaches for face identification

venue: In International Conference on Computer Vision,

authors:
Andi Han
Junbin Gao

title: Improved variance reduction methods for Riemannian non-convex optimization

venue: IEEE Transactions on Pattern Analysis and Machine Intelligence,

authors:
Andi Han
Bamdev Mishra
Pratik Jawanpuria
Junbin Gao

title: Generalized Bures-Wasserstein geometry for positive definite matrices

title: On Riemannian optimization over positive definite matrices with the Bures-Wasserstein geometry

authors:
Mehrtash T Harandi
Mathieu Salzmann
Richard Hartley

title: From manifold to manifold: Geometry-aware dimensionality reduction for SPD matrices

venue: In European conference on computer vision,

authors:
Mehrtash T Harandi
Conrad Sanderson
Richard Hartley
Brian C Lovell

title: Sparse coding and dictionary learning for symmetric positive definite matrices: A kernel approach

venue: In European Conference on Computer Vision,

authors:
Reshad Hosseini
Mohamadreza

title: Mash’al. MixEst: An estimation toolbox for mixture models

authors:
Reshad Hosseini
Suvrit Sra

title: An alternative to EM for Gaussian mixture models: batch and stochastic Riemannian optimization

authors:
Wen Huang
Pierre-Antoine Absil
Kyle A. Gallivan
Paul Hand

title: Roptlib: an object-oriented C++ library for optimization on Riemannian manifolds

venue: Technical Report FSU16-14.v2,

authors:
Wen Huang
Kyle A Gallivan
P-A Absil

title: A Broyden class of quasi-Newton methods for Riemannian optimization

venue: SIAM Journal on Optimization,

authors:
Zhiwu Huang
Luc Van Gool

title: A Riemannian network for SPD matrix learning

authors:
Zhiwu Huang
Ruiping Wang
Xianqiu Li
Wenxian Liu
Shiguang Shan
Luc Van Gool
Xilin Chen

title: Geometry-aware similarity learning on SPD manifolds for visual recognition

venue: IEEE Transactions on Circuits and Systems for Video Technology,

authors:
Zhiwu Huang
Ruiping Wang
Shiguang Shan
Xianqiu Li
Xilin Chen

title: Log-Euclidean metric learning on symmetric positive definite manifold with application to image set classification

venue: In International Conference on Machine Learning,

authors:
Pratik Jawanpuria
Arjun Balgovind
Anoop Kunchukuttan
Bamdev Mishra

title: Learning multilingual word embeddings in latent metric space: A geometric approach

venue: Transactions of the Association for Computational Linguistics,

authors:
Pratik Jawanpuria
Maksim Lapin
Matthias Hein
Bernt Schiele

title: Efficient output kernel learning for multiple tasks

venue: In Conference on Neural Information Processing Systems,

authors:
Pratik Jawanpuria
Bamdev Mishra

title: A unified framework for structured low-rank matrix learning

venue: In ICML,

authors:
Pratik Jawanpuria
Satya Dev N T V
Anoop Kunchukuttan
Bamdev Mishra

title: Learning geometric word meta-embeddings

venue: In ACL Workshop on Representation Learning for NLP,

authors:
Pratik Jawanpuria
Satya Dev N T V
Bamdev Mishra

title: Efficient robust optimal transport with application to multi-label classification

venue: In IEEE CDC,

authors:
Sadeep Jayasumana
Richard Hartley
Mathieu Salzmann
Hongdong Li
Mehrtash Harandi

title: Kernel methods on the Riemannian manifold of symmetric positive definite matrices

venue: In Conference on Computer Vision and Pattern Recognition,

authors:
Michel Journée
Francis Bach
P-A Absil
Rodolphe Sepulchre

title: Low-rank optimization on the cone of positive semidefinite matrices

title: Explicit solutions of linear matrix equations

venue: SIAM Review,

authors:
Yiding Lin
Valeria Simoncini

title: A new subspace iteration method for the algebraic riccati equation

venue: Numerical Linear Algebra with Applications,

title: Riemannian geometry of symmetric positive definite matrices via cholesky decomposition

authors:
Ilia A. Luchnikov
Alexander Ryzhov
Sergey N. Filippov
Henni Ouerdane

title: QGOpt: Riemannian optimization for quantum technologies

venue: SciPost Phys.,

authors:
Sridhar Mahadevan
Bamdev Mishra
Shalini Ghosh

title: A unified framework for domain adaptation using metric learning on manifolds

venue: In ECML-PKDD,

authors:
Luigi Malagò
Luigi Montrucchio
Giovanni Pistone

title: Wasserstein Riemannian geometry of Gaussian densities

venue: Information Geometry,

authors:
Mayank Meghwanshi
Pratik Jawanpuria
Anoop Kunchukuttan
Hiroyuki Kasai
Bamdev Mishra

title: Mctorch, a manifold optimization library for deep learning

authors:
Gilles Meyer
Silvère Bonnabel
Rodolphe Sepulchre

title: Regression on fixed-rank positive semidefinite matrices: a Riemannian approach

title: Toponogov’s theorem and applications

venue: Lecture Notes,

authors:
Bamdev Mishra
Hiroyuki Kasai
Pratik Jawanpuria

title: Riemannian optimization on the simplex of positive definite matrices

venue: NeurIPS workshop on Optimization for Machine Learning,

authors:
Madhav Nimishakavi
Pratik Jawanpuria
Bamdev Mishra

title: A dual framework for low-rank tensor completion

venue: In NeurIPS,

title: Manifold-valued image processing with SPD matrices. In Riemannian Geometric Statistics in Medical Image Analysis, pages 75–134

authors:
Xavier Pennec
Pierre Fillard
Nicholas Ayache

title: A Riemannian framework for tensor computing

venue: International Journal of computer vision,

authors:
Minh Ha Quang
Marco San Biagio
Vittorio Murino

title: Log-hilbert-schmidt metric between positive definite operators on hilbert spaces

authors:
D Rothschild
A Jameson

title: Comparison of four numerical algorithms for solving the Liapunov matrix equation

venue: International Journal of Control,

authors:
Hiroyuki Sato
Hiroyuki Kasai
Bamdev Mishra

title: Riemannian stochastic variance reduced gradient algorithm with retraction and vector transport

authors:
Bernhard Schölkopf
Alexander J Smola

title: Learning with kernels: support vector machines, regularization, optimization, and beyond

venue: MIT press,

authors:
Martin Slawski
Ping Li
Matthias Hein

title: Regularization-free estimation in trace regression with symmetric positive semidefinite matrices

title: The positive definite completion problem revisited

title: A new metric on the manifold of kernel matrices with application to matrix geometric means

authors:
Suvrit Sra
Reshad Hosseini

title: Conic geometric optimization on the manifold of positive definite matrices

authors:
Juan Luis Suárez
Salvador García
Francisco Herrera

title: A tutorial on distance metric learning: Mathematical foundations, algorithms, experimental analysis, prospects and challenges

authors:
Yann Thanwerdas
Xavier Pennec

title: Is affine-invariance well defined on SPD matrices? a principled continuum of metrics

authors:
James Townsend
Niklas Koep
Sebastian Weichwald

title: Pymanopt: A Python toolbox for optimization on manifolds using automatic differentiation

authors:
Koji Tsuda
Gunnar Rätsch
Manfred K. Warmuth

title: Matrix exponentiated gradient updates for on-line learning and Bregman projection

venue: Journal of Machine Learning Research,

authors:
Oncel Tuzel
Fatih Porikli
Peter Meer

title: Pedestrian detection via classification on Riemannian manifolds

title: Convex functions and optimization methods on Riemannian manifolds, volume 297

authors:
Bart Vandereycken
Stefan Vandewalle

title: A Riemannian optimization approach for computing low-rank solutions of Lyapunov equations

title: Geodesic convex optimization: Differentiation on manifolds, geodesics, and convexity

authors:
Chengjing Wang
Defeng Sun
Kim-Chuan Toh

title: Solving log-determinant optimization problems by a Newton-CG primal proximal point algorithm

authors:
Hongyi Zhang
Sashank J Reddi
Suvrit Sra

title: Riemannian SVRG: Fast stochastic optimization on Riemannian manifolds

authors:
Hongyi Zhang
Suvrit Sra

title: First-order methods for geodesically convex optimization

venue: In Conference on Learning Theory,

authors:
Jianjia Zhang
Lei Wang
Luping Zhou
Wanqing Li

title: Learning discriminative Stein kernel for SPD matrices and its applications

venue: IEEE Transactions on Neural Networks and Learning Systems,

heading: 1 Introduction

text:
Learning on symmetric positive definite (SPD) matrices is a fundamental problem in various machine learning applications, including metric and kernel learning [69, 24, 38, 13, 66], medical imaging [57, 56], natural language processing [37, 40], computer vision [28, 35, 34], multi-task learning [39, 54], domain adaptation [48, 41], modeling time-varying data [18], object detection [70], and quantum mechanics [53, 47], among others.
The set of SPD matrices of size n× n, defined as Sn++ := {X : X ∈ Rn×n,X> = X, and X 0}, has a smooth manifold structure with a richer geometry than the Euclidean space. When endowed with a metric (inner product structure), the set of SPD matrices becomes a Riemannian manifold [11]. Hence, numerous existing works [57, 7, 42, 36, 34, 46, 56] have studied and employed the Riemannian optimization framework for learning over the space of SPD matrices [2, 16].
Several Riemannian metrics on Sn++ have been proposed such as the Affine-Invariant [57, 11], the Log-Euclidean [7, 58], the Log-Det [64, 20], and the Log-Cholesky [46], to name a few. One can additionally obtain different families of Riemannian metrics on Sn++ by appropriate parameterizations based on the principles of invariance and symmetry [23, 20, 67, 56]. However, to the best of our knowledge, a systematic study comparing the different metrics for optimizing generic cost functions defined on Sn++ is missing. Practically, the Affine-Invariant (AI) metric seems to be the most widely used metric in Riemannian first- and second-order algorithms (e.g., steepest descent, conjugate gradients, and trust regions) as it is the only Riemannian SPD metric available in manifold optimization toolboxes, such as Manopt [17], Manopt.jl [10], Pymanopt [68], ROPTLIB [32], and McTorch [50]. Moreover, many interesting problems in machine learning are found to be geodesic convex (generalization of Euclidean convexity) under the AI metric, which allows fast convergence of optimization algorithms [76, 31].
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Recent works have studied the Bures-Wasserstein (BW) distance on SPD matrices [49, 12, 72]. It is a well-known result that the Wasserstein distance between two multivariate Gaussian densities is a function of the BW distance between their covariance matrices. Indeed, the BW metric is a Riemannian metric. Under this metric, the necessary tools for Riemannian optimization, including the Riemannian gradient and Hessian expressions, can be efficiently computed [49]. Hence, it is a promising candidate for Riemannian optimization on Sn++. In this work, we theoretically and empirically analyze the quality of optimization with the BW geometry and show that it is a viable alternative to the default choice of AI geometry. Our analysis discusses the classes of cost functions (e.g., polynomial) for which the BW metric has better convergence rates than the AI metric. We also discuss cases (e.g., log-det) where the reverse is true. In particular, our contributions are as follows.
• We observe that the BW metric has a linear dependence on SPD matrices while the AI metric has a quadratic dependence. We show this impacts the condition number of the Riemannian Hessian and makes the BW metric more suited to learning ill-conditioned SPD matrices than the AI metric.
• In contrast to the non-positively curved AI geometry, the BW geometry is shown to be non-negatively curved, which leads to a tighter trigonometry distance bound and faster convergence rates for optimization algorithms.
• For both metrics, we analyze the convergence rates of Riemannian steepest descent and trust region methods and highlight the issues arising from the differences in the curvature and condition number of the Riemannian Hessian.
• We verify that common optimization problems that are geodesic convex under the AI metric are also geodesic convex under the BW metric.
• We support our analysis with extensive experiments on applications such as weighted least squares, trace regression, metric learning, and Gaussian mixture model.

heading: 2 Preliminaries

text:
The fundamental ingredients for Riemannian optimization are Riemannian metric, exponential map, Riemannian gradient, and Riemannian Hessian. We refer readers to [2, 16] for a general treatment on Riemannian optimization.
A Riemannian metric is a smooth, bilinear, and symmetric positive definite function on the tangent space TxM for any x ∈ M. That is, g : TxM× TxM −→ R, which is often written as an inner product 〈·, ·〉x. The induced norm of a tangent vector u ∈ TxM is given by ‖u‖x = √ 〈u, u〉
x . A
geodesic on a manifold γ : [0, 1] −→M is defined as a locally shortest curve with zero acceleration. For any x ∈ M, u ∈ TxM, the exponential map, Expx : TxM −→ M is defined such that there exists a geodesic curve γ with γ(0) = x, γ(1) = Expx(u) and γ ′(0) = u.
First-order geometry and Riemannian steepest descent. Riemannian gradient of a differentiable function f :M−→ R at x, denoted as gradf(x), is a tangent vector that satisfies for any u ∈ TxM, 〈gradf(x), u〉x = Duf(x), where Duf(x) is the directional derivative of f(x) along u. The Riemannian steepest descent method [71] generalizes the standard gradient descent in the Euclidean space to Riemannian manifolds by ensuring that the updates are along the geodesic and stay on the manifolds. That is, xt+1 = Expxt(−ηt gradf(xt)) for some step size ηt.
Second-order geometry and Riemannian trust region. Second-order methods such as trust region and cubic regularized Newton methods are generalized to Riemannian manifolds [1, 3]. They make use of the Riemannian Hessian, Hessf(x) : TxM−→ TxM, which is a linear operator that is defined as the covariant derivative of the Riemannian gradient. Both the trust region and cubic regularized Newton methods are Hessian-free in the sense that only evaluation of the Hessian acting on a tangent vector, i.e., Hessf(x)[u] is required. Similar to the Euclidean counterpart, the Riemannian trust region method approximates the Newton step by solving a subproblem, i.e.,
min u∈TxtM:‖u‖xt≤∆
mxt(u) = f(xt) + 〈gradf(xt), u〉xt + 1
2 〈Hxt [u], u〉xt ,
whereHxt : TxtM−→ TxtM is a symmetric and linear operator that approximates the Riemannian Hessian. ∆ is the radius of trust region, which may be increased or decreased depending on how model value mxt(u) changes. The subproblem is solved iteratively using a truncated conjugate gradient algorithm. The next iterate is given by xt+1 = Expxt(u) with the optimized u.
Next, the eigenvalues and the condition number of the Riemannian Hessian are defined as follows, which we use for analysis in Section 3. Definition 1. The minimum and maximum eigenvalues of Hessf(x) are defined as λmin = min‖u‖2x=1〈Hessf(x)[u], u〉x and λmax = max‖u‖2x=1〈Hessf(x)[u], u〉x. The condition number of Hessf(x) is defined as κ(Hessf(x)) := λmax/λmin.
Function classes on Riemannian manifolds. For analyzing algorithm convergence, we require the definitions for several important function classes on Riemannian manifolds, including geodesic convexity and smoothness. Similarly, we require the definition for geodesic convex sets that generalize (Euclidean) convex sets to manifolds [65, 74]. Definition 2 (Geodesic convex set [65, 74]). A set X ⊆M is geodesic convex if for any x, y ∈ X , the distance minimizing geodesic γ joining the two points lies entirely in X .
Indeed, this notion is well-defined for any manifold because a sufficiently small geodesic ball is always geodesic convex. Definition 3 (Geodesic convexity [65, 74]). Consider a geodesic convex set X ⊆ M. A function f : X −→ R is called geodesic convex if for any x, y ∈ X , the distance minimizing geodesic γ joining x and y satisfies ∀ t ∈ [0, 1], f(γ(t)) ≤ (1− t)f(x) + tf(y). Function f is strictly geodesic convex if the equality holds only when t = 0, 1. Definition 4 (Geodesic strong convexity and smoothness [65, 33]). Under the same settings in Definition 3. A twice-continuously differentiable function f : X −→ R is called geodesic µ-strongly convex if for any distance minimizing geodesic γ in X with ‖γ′(0)‖ = 1, it satisfies d
2f(γ(t)) dt2 ≥ µ,
for some µ > 0. Function f is called geodesic L-smooth if d 2f(γ(t)) dt2 ≤ L, for some L > 0.

heading: 3 Comparing BW with AI for Riemannian optimization

text:
This section starts with an observation of a linear-versus-quadratic dependency between the two metrics. From this observation, we analyze the condition number of the Riemannian Hessian. Then, we further compare the sectional curvature of the two geometries. Together with the differences in the condition number, this allows us to compare the convergence rates of optimization algorithms on the two geometries. We conclude this section by showing geodesic convexity of several generic cost functions under the BW geometry. The proofs for this section are in our extended report [27].
AI and BW geometries on SPD matrices. When endowed with a Riemannian metric g, the set of SPD matrices of size n becomes a Riemannian manifoldM = (Sn++, g). The tangent space at X is TXM := {U : U ∈ Rn×n and U> = U}. Under the AI and BW metrics, the Riemannian exponential map, Riemannian gradient, and Hessian are compared in Table 1, where we denote {A}S := (A + A>)/2 and exp(A) as the matrix exponential of A. LX[U] is the solution to the matrix linear system LX[U]X + XLX[U] = U and is known as the Lyapunov operator. We use
∇f(X) and ∇2f(X) to represent the first-order and second-order derivatives, i.e., the Euclidean gradient and Hessian, respectively. The derivations in Table 1 can be found in [56, 12]. In the rest of the paper, we useMai andMbw to denote the SPD manifolds under the two metrics. From Table 1, the computational costs for evaluating the AI and BW ingredients are dominated by the matrix exponential/inversion operations and the Lyapunov operator L computation, respectively. Both at most cost O(n3), which implies a comparable per-iteration cost of optimization algorithms between the two metric choices. This claim is validated in Section 4.
A key observation. From Table 1, the Affine-Invariant metric on the SPD manifold can be rewritten as for any U,V ∈ TXM,
〈U,V〉ai = tr(X−1UX−1V) = vec(U)>(X⊗X)−1vec(V), (1)
where vec(U) and vec(V) are the vectorizations of U and V, respectively. Note that we omit the subscript X for inner product 〈·, ·〉 to simplify the notation. The specific tangent space where the inner product is computed should be clear from contexts.
The Bures-Wasserstein metric is rewritten as, for any U,V ∈ TXM,
〈U,V〉bw = 1
2 tr(LX[U]V) =
1 2 vec(U)>(X⊕X)−1vec(V), (2)
where X⊕X = X⊗ I + I⊗X is the Kronecker sum. Remark 1. Comparing Eq. (1) and (2) reveals that the BW metric has a linear dependence on X while the AI metric has a quadratic dependence. This suggests that optimization algorithms under the AI metric should be more sensitive to the condition number of X compared to the BW metric.
The above observation serves as a key motivation for further analysis.

heading: 3.1 Condition number of Riemannian Hessian at optimality

text:
Throughout the rest of the paper, we make the following assumptions. Assumption 1. (a). f is at least twice continuously differentiable with a non-degenerate local minimizer X∗. (b). The subset X ⊆M (usually as a neighbourhood of a center point) we consider throughout this paper is totally normal, i.e., the exponential map is a diffeomorphism.
Assumption 1 is easy to satisfy. Particularly, Assumption 1(b) is guaranteed for the SPD manifold under the AI metric because its geodesic is unique. Under the BW metric, for a center point X, we can choose the neighbourhood such that X = {ExpX(U) : I + LX[U] ∈ Sn++} as in [49]. In other words, X is assumed to be unique-geodesic under both the metrics. We now formalize the impact of the linear-versus-quadratic dependency, highlighted in Remark 1. At a local minimizer X∗ where the Riemannian gradient vanishes, we first simplify the expression for the Riemannian Hessian in Table 1.
OnMai, 〈Hessaif(X∗)[U],U〉ai = tr(∇2f(X∗)[U]U) = vec(U)>H(X∗)vec(U), with H(X) ∈ Rn2×n2 is the matrix representation of the Euclidean Hessian ∇2f(X) and U ∈ TX∗Mai. The maximum eigenvalue of Hessaif(X∗) is then given by λ∗max = max‖U‖2ai=1 vec(U)
>H(X∗)vec(U), where ‖U‖2ai = vec(U)>(X∗ ⊗X∗)−1vec(U). This is a generalized eigenvalue problem with the solution to be the maximum eigenvalue of (X∗ ⊗X∗)H(X∗). Similarly, λ∗min corresponds to the minimum eigenvalue of (X∗ ⊗X∗)H(X∗).
On Mbw, 〈Hessbwf(X∗)[U],U〉bw = vec(U)>H(X∗)vec(U) and the norm is ‖U‖2bw = vec(U)>(X∗ ⊕X∗)−1vec(U). Hence, the minimum/maximum eigenvalue of Hessbwf(X∗) equals the minimum/maximum eigenvalue of (X∗ ⊕X∗)H(X∗). Let κ∗ai := κ(Hessaif(X
∗)) = κ((X∗ ⊗X∗)H(X∗)) and κ∗bw := κ(Hessbwf(X∗)) = κ((X∗ ⊕ X∗)H(X∗)). The following lemma bounds these two condition numbers. Lemma 1. For a local minimizer X∗ of f(X), the condition number of Hessf(X∗) satisfies
κ(X∗)2/κ(H(X∗)) ≤κ∗ai ≤ κ(X∗)2κ(H(X∗)) κ(X∗)/κ(H(X∗)) ≤κ∗bw ≤ κ(X∗)κ(H(X∗)).
It is clear that κ∗bw ≤ κ∗ai when κ(H(X∗)) ≤ √ κ(X∗). This is true for linear, quadratic, higher-order polynomial functions and in general holds for several machine learning optimization problems on the SPD matrices (discussed in Section 4). Case 1 (Condition number for linear and quadratic optimization). For a linear function f(X) = tr(XA), its Euclidean Hessian matrix is H(X) = 0. For a quadratic function f(X) = tr(XAXB) with A,B ∈ Sn++, H(X) = A ⊗ B + B ⊗ A. Therefore, κ(H(X∗)) is a constant and for ill-conditioned X∗, we have κ(H(X∗)) ≤ √ κ(X∗), which leads to κ∗ai ≥ κ∗bw. Case 2 (Condition number for higher-order polynomial optimization). For an integer α ≥ 3, consider a function f(X) = tr(Xα) with derived H(X) = α ∑α−2 l=0 (X
l ⊗ Xα−l−2). We get κ∗ai = α ∑α−1 l=1 ((X ∗)l ⊗ (X∗)α−l) and κ∗bw = α(X∗ ⊕X∗)( ∑α−2 l=0 (X
l ⊗Xα−l−2)). It is apparent that κ∗ai = O(κ(X∗)α) while κ∗bw = O(κ(X∗)α−1). Hence, for ill-conditioned X∗, κ∗ai ≥ κ∗bw.
One counter-example where κ∗bw ≥ κ∗ai is the log-det function. Case 3 (Condition number for log-det optimization). For the log-det function f(X) = − log det(X), its Euclidean Hessian is∇2f(X)[U] = X−1UX−1 and H(X) = X−1⊗X−1. At a local minimizer X∗, Hessaif(X∗)[U] = U with κ∗ai = 1. While onMbw, we have κ∗bw = κ((X∗ ⊕X∗)((X∗)−1 ⊗ (X∗)−1)) = κ((X∗)−1 ⊕ (X∗)−1) = κ(X∗). Therefore, κ∗ai ≤ κ∗bw.

heading: 3.2 Sectional curvature and trigonometry distance bound

text:
To study the curvature ofMbw, we first show in Lemma 2, the existence of a matching geodesic between the Wasserstein geometry of zero-centered non-degenerate Gaussian measures and the BW geometry of SPD matrices. Denote the manifold of such Gaussian measures under the L2-Wasserstein distance as (N0(Σ),W2) with Σ ∈ Sn++. Lemma 2. For any X,Y ∈ Sn++, a geodesic betweenN0(X) andN0(Y) on (N0(Σ),W2) is given by N0(γ(t)), where γ(t) is the geodesic between X and Y onMbw.
The following lemma builds on a result from the Wasserstein geometry [6] and uses Lemma 2 to analyze the sectional curvature ofMbw. Lemma 3. Mbw is an Alexandrov space with non-negative sectional curvature.
It is well-known thatMai is a non-positively curved space [21, 56] while, in Lemma 3, we show that Mbw is non-negatively curved. The difference affects the curvature constant in the trigonometry distance bound of Alexandrov space [77]. This bound is crucial in analyzing convergence for optimization algorithms on Riemannian manifolds [77, 76]. In Section 3.3, only local convergence to a minimizer X∗ is analyzed. Therefore, it suffices to consider a neighbourhood Ω around X∗. In such a compact set, the sectional curvature is known to be bounded and we denote the lower bound as K−.
The following lemma compares the trigonometry distance bounds under the AI and BW geometries. This bound was originally introduced for Alexandrov space with lower bounded sectional curvature [77]. The result for non-negatively curved spaces has been applied in many work [76, 60, 25] though without a formal proof. We show the proof in [27], where it follows from the Toponogov comparison theorem [52] on the unit hypersphere and Assumption 1. Lemma 4. Let X,Y,Z ∈ Ω, which forms a geodesic triangle onM. Denote x = d(Y,Z), y = d(X,Z), z = d(X,Y) as the geodesic side lengths and let θ be the angle between sides y and z such that cos(θ) = 〈Exp−1X (Y),Exp −1 X (Z)〉/(yz). Then, we have
x2 ≤ ζy2 + z2 − 2yz cos(θ),
where ζ is a curvature constant. Under the AI metric, ζ = ζai = √ |K−ai |D
tanh( √ |K−ai |D) with D as the
diameter bound of Ω, i.e. maxX1,X2∈Ω d(X,Y) ≤ D. Under the BW metric, ζ = ζbw = 1.
It is clear that ζai > ζbw = 1, which leads to a tighter bound under the BW metric.

heading: 3.3 Convergence analysis

text:
We now analyze the local convergence properties of the Riemannian steepest descent and trust region methods under the two Riemannian geometries. Convergence is established in terms of the
Riemannian distance induced from the geodesics. We begin by presenting a lemma that shows in a neighbourhood of X∗, the second-order derivatives of f ◦ ExpX are both lower and upper bounded. Lemma 5. In a totally normal neighbourhood Ω around a non-degenerate local minimizer X∗, for any X ∈ Ω, it satisfies that λ∗min/α ≤ d 2 dt2 f(ExpX(tU)) ≤ αλ ∗ max, for some α ≥ 1 and ‖U‖ = 1. λ∗max > λ ∗ min > 0 are the largest and smallest eigenvalues of Hessf(X ∗).
For simplicity of the analysis, we assume such an α is universal under both the Riemannian geometries. We, therefore, can work with a neighbourhood Ω with diameter uniformly bounded by D, where we can choose D := min{Dai, Dbw} such that α is universal. One can readily check that under Lemma 5 the function f is both µ-geodesic strongly convex and L-geodesic smooth in Ω where µ = λ∗min/α and L = αλ ∗ max. We now present the local convergence analysis of the two algorithms, which are based on results in [77, 1]. Theorem 1 (Local convergence of Riemannian steepest descent). Under Assumption 1 and consider a non-degenerate local minimizer X∗. For a neighbourhood Ω 3 X∗ with diameter bounded by D on two Riemannian geometriesMai,Mbw, running Riemannian steepest descent from X0 ∈ Ω with a fixed step size η = 1αλ∗max yields for t ≥ 2,
d2(Xt,X ∗) ≤ α2D2κ∗ (1−min{ 1
ζ ,
1
α2κ∗ })t−2 .
Theorem 2 (Local convergence of Riemannian trust region). Under the same settings as in Theorem 1, assume further in Ω, it holds that (1) ‖HXt − Hessf(Xt)‖ ≤ `‖gradf(Xt)‖ and (2) ‖∇2(f ◦ ExpXt)(U)−∇
2(f ◦ ExpXt)(0)‖ ≤ ρ‖U‖ for some `, ρ universal onMai,Mbw. Then running Riemannian trust region from X0 ∈ Ω yields,
d(Xt,X ∗) ≤ (2√ρ+ `)(κ∗)2d2(Xt−1,X∗).
Theorems 1 and 2 show thatMbw has a clear advantage compared toMai for learning ill-conditioned SPD matrices where κ∗bw ≤ κ∗ai. For first-order algorithms,Mbw has an additional benefit due to its non-negative sectional curvature. As ζai > ζbw = 1, the convergence rate degrades on Mai. Although the convergence is presented in Riemannian distance, it can be readily converted to function value gap by noticing µ2 d 2(Xt,X ∗) ≤ f(Xt)− f(X∗) ≤ L2 d 2(Xt,X ∗). Additionally, we note that these local convergence results hold regardless of whether the function is geodesic convex or not, and similar comparisons also exist for other Riemannian optimization methods.

heading: 3.4 Geodesic convexity under BW metric for cost functions of interest

text:
Finally we show geodesic convexity of common optimization problems onMbw. Particularly, we verify that linear, quadratic, log-det optimization, and also certain geometric optimization problems, that are geodesic convex under the AI metric, are also geodesic convex under the BW metric. Proposition 1. For any A ∈ Sn+, where Sn+ := {Z : Z ∈ Rn×n,Z> = Z, and Z 0}, the functions f1(X) = tr(XA), f2(X) = tr(XAX), and f3(X) = − log det(X) are geodesic convex onMbw.
Based on the result in Proposition 1, we also prove geodesic convexity of a reparameterized version of the Gaussian density estimation and mixture model onMbw (discussed in Section 4). Similar claims onMai can be found in [31]. We further show that monotonic functions on sorted eigenvalues are geodesic convex onMbw. This is an analogue of [65, Theorem 2.3] onMai. Proposition 2. Let λ↓ : Sn++ −→ Rn+ be the decreasingly sorted eigenvalue map and h : R+ −→ R be an increasing and convex function. Then f(X) = ∑k j=1 h(λ ↓ j (X)) for 1 ≤ k ≤ n is geodesic convex onMbw. Examples of such functions include f1(X) = tr(exp(X)) and f2(X) = tr(Xα), α ≥ 1.

heading: 4 Experiments

text:
In this section, we compare the empirical performance of optimization algorithms under different Riemannian geometries for various problems. In addition to AI and BW, we also include the Log-Euclidean (LE) geometry [7] in our experiments.
The LE geometry explores the the linear space of symmetric matrices where the matrix exponential acts as a global diffeomorphism from the space to Sn++. The LE metric is defined as
〈U,V〉le = tr(DU log(X)DV log(X)) (3)
for any U,V ∈ TXM, where DU log(X) is the directional derivative of matrix logarithm at X along U. Following [69, 51, 58], for deriving various Riemannian optimization ingredients under the LE metric (3), we consider the parameterization X = exp(S), where S ∈ Sn, i.e., the space of n× n symmetric matrices. Equivalently, optimization on the SPD manifold with the LE metric is identified with optimization on Sn and the function of interest becomes f(exp(S)) for S ∈ Sn. While the Riemannian gradient can be computed efficiently by exploiting the directional derivative of the matrix exponential [4], deriving the Riemannian Hessian is tricky and we rely on finite-difference Hessian approximations [15].
We present convergence mainly in terms of the distance to the solution X∗ whenever applicable. The distance is measured in the Frobenius norm, i.e., ‖Xt−X∗‖F. When X∗ is not known, convergence is shown in the modified Euclidean gradient norm ‖Xt∇f(Xt)‖F. This is comparable across different metrics as the optimality condition X∗∇f(X∗) = 0 arises from problem structure itself [43]. We initialize the algorithms with the identity matrix for the AI and BW metrics and zero matrix for the LE metric (i.e., the matrix logarithm of the identity).
We mainly present the results on the Riemannian trust region (RTR) method, which is the method of choice for Riemannian optimization. Note for RTR, the results are shown against the cumulative sum of inner iterations (which are required to solve the trust region subproblem at every iteration). We also include the Riemannian steepest descent (RSD) and Riemannian stochastic gradient (RSGD) [14] methods for some examples. The experiments are conducted in Matlab using the Manopt toolbox [17] on a i5-10500 3.1GHz CPU processor.
In our extended report [27], we include additional experiments comparing convergence in objective function values for the three geometries. We also present results for the Riemannian conjugate gradient method, and results with different initializations (other than the identity and zero matrices) to further support our claims.
The code can be found at https://github.com/andyjm3/AI-vs-BW.
Weighted least squares. We first consider the weighted least squares problem with the symmetric positive definite constraint. The optimization problem is minX∈Sn++ f(X) = 1 2‖A X−B‖ 2 F, which is encountered in for example, SPD matrix completion [63] where A is a sparse matrix. The Euclidean gradient and Hessian are∇f(X) = (A X−B) A and∇2f(X)[U] = A U A, respectively. Hence, at optimal X∗, the Euclidean Hessian in matrix representation is H(X∗) = diag(vec(A A)). We experiment with two choices of A, i.e. A = 1n1>n (Dense) and A as a random sparse matrix (Sparse). The former choice for A leads to well-conditioned H(X∗) while the latter choice leads to an ill-conditioned H(X∗). Also note that when A = 1n1>n , κ ∗ ai = κ(X ∗)2 and κ∗bw = κ(X ∗). We generate X∗ as a SPD matrix with size n = 50 and exponentially decaying eigenvalues. We consider two cases with condition numbers κ(X∗) = 10 (LowCN) and 103 (HighCN). The matrix B is generated as B = A X∗. Figure 1 compares both RSD and RTR for different metrics. When A is either dense or sparse, convergence is significantly faster onMbw than bothMai andMle. The advantage of usingMbw becomes more prominent in the setting when condition number of X∗ is high. Figure 1(e) shows thatMbw is also superior in terms of runtime.
Lyapunov equations. Continuous Lyapunov matrix equation, AX + XA = C with X ∈ Sn++, is commonly employed in analyzing optimal control systems and differential equations [59, 45]. When A is stable, i.e., λi(A) > 0 and C ∈ Sn++, the solution X∗ 0 and is unique [44]. When C ∈ Sn+ and is low rank, X∗ ∈ Sn+ is also low rank. We optimize the following problem for solving the Lyapunov equation [73], i.e., minX∈Sn++ f(X) = tr(XAX)− tr(XC). The Euclidean gradient and Hessian are respectively ∇f(X) = AX + XA − C and ∇2f(X)[U] = AU + UA with H(X) = A⊕A. At optimal X∗, the condition number κ(H(X∗)) = κ(A). We experiment with two settings for the matrix A, i.e. A as the Laplace operator on the unit square where we generate 7 interior points so that n = 49 (Ex1), and A is a particular Toeplitz matrix with n = 50 (Ex2). The generated A matrices are ill-conditioned. The above settings correspond to Examples 7.1 and 7.3 in [45]. Under each setting, X∗ is set to be either full or low rank. The matrix C is generated as C = AX∗ + X∗A. The full rank X∗ is generated from the full-rank Wishart distribution while the low rank X∗ is a diagonal matrix with r = 10 ones and n− r zeros in the diagonal. We label the four cases as Ex1Full, Ex1Low, Ex2Full, and Ex2Low. The results are shown in Figures 2(a)-(d), where we observe that in all four cases, the BW geometry outperforms both AI and LE geometries.
Trace regression. We consider the regularization-free trace regression model [62] for estimating covariance and kernel matrices [61, 19]. The optimization problem is minX∈Sd++ f(X) =
1 2m ∑m i=1(yi − tr(A>i X))2, where Ai = aia>i , i = 1, ...,m are some rank-one measurement
matrices. Thus, we have∇f(X) = ∑m i=1(a > i Xai − yi)Ai and∇2f(X)[U] = ∑m i=1(a > i Uai)Ai. We create X∗ as a rank-r Wishart matrix and {Ai} as rank-one Wishart matrices and generate yi = tr(AiX
∗) + σ i with i ∼ N (0, 1), σ = 0.1. We consider two choices, (m, d, r) = (1000, 50, 50) and (1000, 50, 10), which are respectively labelled as SynFull and SynLow. From Figures 2(e)&(f), we also observe that convergence to the optimal solution is faster for the BW geometry.
Metric learning. Distance metric learning (DML) aims to learn a distance function from samples and a popular family of such distances is the Mahalanobis distance, i.e. dM(x,y) =√
(x− y)>M(x− y) for any x,y ∈ Rd. The distance is parameterized by a symmetric positive semi-definite matrix M. We refer readers to this survey [66] for more discussions on this topic. We particularly consider a logistic discriminant learning formulation [24]. Given a training sample {xi, yi}Ni=1, denote the link tij = 1 if yi = yj and tij = 0 otherwise. The objective is given by minM∈Sd++ f(M) = − ∑ i,j (tij log pij + (1− tij) log(1− pij)) , with pij = (1 + exp(dM(xi,xj))) −1. We can derive the matrix Hessian as H(M) = ∑ i,j pij(1 − pij)(xi − xj)(xi−xj)>⊗(xi−xj)(xi−xj)>. Notice κ(H(M∗)) depends on M∗ only through the constants pij . Thus, the condition number will not be much affected by κ(M∗).
We consider two real datasets, glass and phoneme, from the Keel database [5]. The number of classes is denoted as c. The statistics of these two datasets are (N, d, c) = (241, 9, 7) for glass
(5404, 5, 2) for phoneme. In Figures 2(g)&(h), we similarly see the advantage of using the BW metric compared to the other two metrics that behave similarly.
Log-det maximization. As discussed in Section 3.1, log-det optimization is one instance where κ∗bw ≥ κ∗ai. We first consider minimizing negative log-determinant along with a linear function as studied in [75]. That is, for some C ∈ Sn++, the objective is minX∈Sn++ f(X) = tr(XC)− log det(X). The Euclidean gradient and Hessian are given by ∇f(X) = C − X−1 and ∇2f(X)[U] = X−1UX−1. This problem is geodesic convex under both AI and BW metrics. We generate X∗ the same way as in the example of weighted least square with n = 50 and set C = (X∗)−1. We consider two cases with condition number cn = 10 (LowCN) and 103 (HighCN). As expeted, we observe faster convergence of AI and LE metrics over the BW metric in Figures 3(a)&(b). This is even more evident when the condition number increases.
Gaussian mixture model. Another notable example of log-det optimization is the Gaussian density estimation and mixture model problem. Following [31], we consider a reformulated problem on augmented samples y>i = [x > i ; 1], i = 1, ..., N where xi ∈ Rd are the original samples. The density is parameterized by the augmented covariance matrix Σ ∈ Rd+1. Notice that the log-likelihood of Gaussian is geodesic convex onMai, but not onMbw. We, therefore, define S = Σ−1 and the reparameterized log-likelihood is pN (Y; S) = ∑N i=1 log ( (2π)1−d/2 exp(1/2) det(S)1/2 exp(− 12y > i Syi) ) , which is now geodesic convex onMbw due to Proposition 1. Hence, we can solve the problem of Gaussian mixture model similar as in [31].
Here, we test on a dataset included in the MixEst package [30]. The dataset has 1580 samples in R2 with 3 Gaussian components. In Figure 3(c), we observe a similar pattern with RTR as in the log-det example. We also include performance of RSGD, which is often preferred for large scale problems. We set the batch size to be 50 and consider a decaying step size, with the best initialized step size shown in Figures 3(d)&(e). Following [8], the algorithms are initialized with kmeans++. We find that the AI geometry still maintains its advantage under the stochastic setting.

heading: 5 Conclusion and discussion

text:
In this paper, we show that the less explored Bures-Wasserstein (BW) geometry for SPD matrices should often be the preferable choice than the Affine-Invariant geometry for optimization, particularly for learning ill-conditioned matrices. A systematic analysis shows that the BW metric preserves geodesic convexity of some popular cost functions and leads to better convergence rates in optimization.
We also theoretically discuss a ‘counter-example’ of log-det optimization where the Affine-Invariant (AI) geometry enjoys a better second-order conditioning and validate our findings empirically. This issue is addressed in our recent work [26], where we propose a generalized Bures-Wasserstein (GBW) geometry that is built on a generalization of the Lyapunov operator in the metric:
〈U,V〉gbw = 1
2 vec(U)>(X⊗M + M⊗X)−1vec(V), (4)
where U and V are symmetric matrices and M is a given SPD matrix. When M = I, the metric (4) reduces to the special BW metric (2). The use of the parameter M in (4) allows great flexibility for optimization algorithms. For one, choosing a particular M allows preconditioning the Hessian by locally approximating the AI geometry. This leads to improved convergence for log-det optimization.
Our comparisons between AI and BW geometries are based on optimization over generic cost functions. For specific problems, however, there may exist other alternative metrics that potentially work better. This is an interesting research direction to pursue. We also remark that optimization is not the only area where the AI and BW geometries can be compared. It would be useful to qualitatively compare the two metrics for other learning problems on SPD matrices, such as barycenter learning. In addition, kernel methods have been studied on the SPD manifold [42, 29, 78] that embed SPD matrices to a high dimensional feature space, known as the Reproducing Kernel Hilbert Space (RKHS). Such representations are used for subsequent learning tasks such as clustering or classification. But only a positive definite kernel provides a valid RKHS. We show in appendix Section A that the induced Gaussian kernel based on the BW distance is a positive definite kernel unlike the case for the AI metric. This difference highlights a potential advantage of the BW metric for representation learning on SPD matrices.

heading: A Positive definite kernel from BW geometry

text:
In this section, we show the existence of a positive definite Gaussian kernel on Mbw. Existing works [55, 22] have studied the Wasserstein distance kernel. First, we present the definition of a positive (resp. negative) definite function as in [9]. Definition 5. Consider X to be a nonempty set. A function f : X ×X −→ R is called positive definite if and only if k is symmetric and for all integers m ≥ 2, {x1, ..., xm} ⊆ X and {c1, ..., cm} ⊆ R, it satisfies ∑m i,j=1 cicjf(xi, xj) ≥ 0. A function f is called negative definite if and only if under the
same conditions, it satisfies ∑m i,j=1 cicjf(xi, xj) ≤ 0 with ∑m i=1 ci = 0.
The following theorem shows the Gaussian kernel induced from the BW distance is positive definite on the SPD manifold. Theorem 3. The induced Gaussian kernel k(·, ·) := exp(−d2bw(·, ·)/(2σ2)) is positive definite.
Proof of Theorem 3. From Theorem 4.3 in [42], it suffices to prove the BW distance d2bw is negative definite. Indeed for any m ≥ 2, {X1, ...,Xm} ⊆ Mbw, {c1, ..., cm} ⊆ R with ∑m i=1 ci = 0, we have m∑
i,j=1
cicjd 2 bw(Xi,Xj)
= m∑ j=1 cj m∑ i=1 citr(Xi) + m∑ i=1 ci m∑ j=1 cjtr(Xj)− 2 m∑ i,j=1 cicjtr(X 1/2 i XjX 1/2 i ) 1/2
=− 2 m∑
i,j=1
cicjtr(X 1/2 i XjX 1/2 i ) 1/2 ≤ 0.
This shows d2bw is negative definite, and therefore, the exponentiated Gaussian kernel is positive definite.
Broader Impact
This paper is of a theoretical nature and does not foresee any immediate societal impacts.

title: On Riemannian Optimization over Positive Definite Matrices with the Bures-Wasserstein Geometry