Projected Riemannian Gradient Descent for the Bures-Wasserstein Barycenter: Dimension-Independent Linear Convergence at Unit Step Size

Projected Riemannian Gradient Descent for the Bures–Wasserstein Barycenter: Dimension-Independent Linear Convergence at Unit Step Size A. Afham∗ Centre for Quantum Technologies, National University of Singapore, Singapore

arXiv:2609.03762v1 [cs.LG] 3 Sep 2026

Abstract The computation of the Bures–Wasserstein (BW) barycenter of an ensemble of positive definite matrices arises throughout machine learning, optimal transport, and quantum information. Riemannian gradient descent (RGD) at unit step size—the fixed-point iteration used in practice—converges rapidly, yet existing analyses present a dichotomy: unit-step guarantees carry worst-case exponential dependence on the dimension, while dimension-independent guarantees require small step sizes that forfeit the empirical speed. We resolve this dichotomy, not by improving the guarantees for unit-step RGD, but by proposing a Projected RGD algorithm that achieves dimension-independent linear convergence at unit step size. The achieved rate, (1 − κ−3/2 ), where κ is the condition number of the ensemble, also polynomially improves on the best small-step guarantee (κ3/2 versus κ5/2 iteration complexity). The crux is a novel Projection Lemma: clipping the eigenvalues of a positive matrix to an interval [α, β] is the closed-form, non-expansive (1-Lipschitz) BW-metric projection onto the set {S ∶ αI ≤ S ≤ βI}—a statement which, unlike its known one-sided counterpart, does not follow from convexity. The projection is moreover free: it reuses an eigendecomposition the next iteration must perform in any case, so the projected and unprojected iterations cost the same per step. The same analysis covers the invariant matrix projection problem of Brahmachari et al. (2025), whose fixed-point algorithm we identify as unit-step RGD on a totally geodesic submanifold, thereby extending the dimension-independent guarantee to that setting verbatim.

1

Introduction

We consider two optimization problems on the Bures–Wasserstein (BW) manifold of positive definite matrices. The first is the computation of the BW barycenter of a probability distribution supported on the positive definite cone [AC11]. The problem is prevalent in machine learning and optimal transport—averaging covariance matrices, interpolating Gaussian distributions, and processing positive-definite features in computer vision [AC11; CNR25; CMRS20; BJL19]—and features in quantum information through Bayesian quantum tomography [AKF22], quantum cryptography & Shannon theory [KRS09; Tom15; Bei13; GW15], and quantum resource theories [WY16; Liu+17; WG03]. We first define the squared BW distance between positive semidefinite matrices P, Q ∈ PSD(d): B(P, Q) ≡ d2BW (P, Q) ∶= Tr[P + Q]−2F(P, Q),

√ 1 1 1 1 F(P, Q) ∶= Tr[ P 2 QP 2 ] = ∥P 2 Q 2 ∥ , (1) 1

where F(P, Q) is the quantum fidelity [Uhl76; Joz94] between P and Q, and ∥A∥1 ∶= √ Tr[ A† A] is the trace norm; see Section 2 for notation and preliminaries. ∗

Email: afham@nus.edu.sg

1

(a) unit against small step (d = 2, n = 5, κ = 102)

(b) transient exit from [αI, βI] (κ = 300) 10−1 10−2

10−2

λmin (S ⋆ ) − α

10−3 10−4 unit step RGD (η = 1)

10−6 10

λmin (St) − α

[f(St) − f(S ⋆ )] / [f(S0) − f(S ⋆ )]

10

0

small step RGD (η = α 0/2β 0)

−8

10−4

λmin (St) = α

0

β

−10−4

102

λmax (St)

−10−3 -4.43e-03

−10−2

10−10

−10−1 10

−12

100

101

102

−10

103

101

λmin (St)

100

-2.98e-01

0

0

2

0

4

6

2

4

8

6

8

10

α 10

12

12

iteration t

iteration t

Figure 1: (a) Unit-step against small-step RGD on 30 instances of the pinned family at d = 2, n = 5, κ = 100 (every member with spectrum {α, β} in a Haar-random eigenbasis), from the initialization S0 = α+β I: the 2 relative optimality gap per iteration, translucent per instance with the medians opaque. At the median, the unit step converges to machine precision in 16 iterations, the small step η = α′ /(2β ′ ) of [ACGS21] (marked every 200th iterate) in about 3300. (b) The transient exit of Example 13 (κ = 300): the signed deviation λmin (St ) − α of the unprojected unit-step iterates, on a symmetric log scale. The iterates at t = 1, 2 (circled) fall below the floor α (dotted), by 2.98 × 10−1 and 4.43 × 10−3 , before re-entering at t = 3 and converging to the strictly interior limit λmin (S⋆ ) − α ≈ 2.52 × 10−3 (dash-dotted). The inset shows the full trajectory of both eigenvalues against the well-conditioned band [αI, βI] ∶= {S ∶ αI ≤ S ≤ βI}.

Problem 1. Let (Ri )i∈[n] ⊆ PD(d) be a collection of d × d positive definite matrices and w ∶= (w1 , . . . , wn ) an associated probability vector. Compute the BW barycenter of the distribution (Ri , wi )i∈[n] : 1 n ∑ wi B(Ri , S) ≡ argmin f1 (S). S∈PD(d) 2 i=1 S∈PD(d)

S⋆ ∶= argmin

(2)

Throughout, the ensemble is contained in a spectral interval: αI ≤ Ri ≤ βI for all i ∈ [n], for some 0 < α ≤ β < ∞, and the ratio κ ∶= β/α is the ensemble’s condition number. For a finite ensemble one may take α = mini∈[n] λmin (Ri ) and β = maxi∈[n] λmax (Ri ), and the same containment covers distributions of infinite support. No closed form for the solution is known outside special cases (e.g., √ pairwise commuting Ri ), but the optimum satisfies the fixed-point equation S⋆ = ∑ni=1 wi BJL19]. An algorithm and a dichotomy. (FP) iteration of [ÁDCM16], ⎡

1

1

S⋆2 Ri S⋆2 [AC11;

The standard method for Problem 1 is the fixed-point

−1 ⎢ St+1 ∶= St 2 ⎢⎢∑ wi ⎢ i ⎣

√

⎤2 1 ⎥ − St Ri St ⎥⎥ St 2 , ⎥ ⎦ 1 2

1 2

S0 ∶= I,

(3)

which coincides with Riemannian gradient descent (RGD) on the BW manifold at unit step BW size [ZP19]: St+1 = ExpBW is St [−∇f (St )], where ∇f is the Riemannian gradient and Exp the Riemannian exponential map on the BW manifold. Empirically, the unit-step algorithm converges rapidly, requires no tuning, and outperforms Euclidean gradient methods, SDP solvers, and RGD with smaller step sizes [ACGS21; BRT25]. Current theory does not match this behavior; instead, it offers a dichotomy. Analyses of the unit-step algorithm carry an iteration complexity that grows as κd in the worst case [CMRS20]. Dimension-free guarantees 2

exist, but only at the small step sizes η = O(1/κ) [ACGS21],1 which on typical instances cost two to three orders of magnitude more iterations (Fig. 1(a)). One may thus have the practical step size or the dimension-free guarantee, but not both; Section 1.1 details both lines. This article side-steps this dichotomy by presenting a Projected BW Gradient Descent algorithm which achieves a dimension-independent linear convergence rate at unit step size: each step of the practical algorithm is composed with an exact, non-expansive projection, and the composition retains the empirical speed of the unit step while admitting the dimension-free guarantee. The dimension-free analysis of [ACGS21] requires the iterates to remain in the wellconditioned set [αI, βI] ∶= {S ∶ αI ≤ S ≤ βI} determined by the ensemble. The unit-step iterates cannot escape [αI, βI] through the ceiling:2 S ≤√βI ⇒ ExpBW S [−η∇f (S)] ≤ βI for any η ∈ [0, 1], a consequence of the convexity of S ↦ λmax (S) along the generalized geodesics traced by the gradient √ step [CMRS20; ACGS21; BRT25]. However, the floor enjoys no such protection: S ↦ λmin (S) fails to be concave along generalized geodesics— in contrast to its concavity along barycenters [ACGS21, Thm. 6, Rem. 3]—so no step-sizeindependent lower eigenvalue bound is available. Example 13 (Appendix A) gives an instance (d = 2, n = 3, κ = 300) in which the ensemble, the current iterate, and the barycenter all lie in the interior of [αI, βI], yet the next two unit-step iterates have minimal eigenvalue strictly below α. Thus at unit step size, the arguments underpinning existing dimensionfree rates fail. The exit is transient as the ultimate convergence of unit-step RGD to the barycenter is known. In the aforementioned example, the iterate permanently returned to the well-conditioned set after two steps, converging to a limit that lies strictly above the floor α (see Fig. 1(b)). Moreover, such instances had to be adversarially constructed; on random instances the exit never occurred from the canonical initialization.3 Such transient exits are avoided in the analysis of [ACGS21] by choosing sufficiently small step sizes (η = O(κ−1 )). The cost of the small step sizes is demonstrated in Fig. 1(a) on the smallest instructive case, d = 2: unit-step RGD gains nearly a decimal digit of accuracy per iteration and reaches the double-precision floor within 16 iterations at the median, whereas the small step η = α′ /(2β ′ ) associated with the dimension-free analysis of [ACGS21]4 needs about 3300 iterations for the same accuracy. The factor between the two is essentially independent of the target accuracy and across 30 random instances of this family its median is 202 (range 141 to 399). A step of RGD at any step size evaluates the same map—n transport maps St−1 #Ri , each an eigendecomposition—and the step size only scales the resulting displacement, so a unit step and a small step cost the same to compute. An iteration count is therefore a wall-clock count, and the practitioner who chooses η = 1 saves two to three orders of magnitude of computation, on every instance. This is why the unit step is the algorithm used in practice, and why a guarantee for it is highly desired. As elaborated upon later, we supply this guarantee through a remedy that is standard in constrained optimization: if the iterates do not natively stay in the region of interest, enforce membership by projection—provided the projection is cheap and does not undo the progress of the gradient step. We show that both requirements are met on the BW manifold, at essentially zero additional cost. The projection is constructed in Section 3.2, the full algorithm is stated as Algorithm 1 in Section 3.3, and its convergence is proved in Section 4. Their results allow for a refined condition number κ′ ∶= β ′ /α′ ≤ κ; see Remark 11. Provided the ensemble satisfies Ri ≤ βI, as it does by the choice of β. 3 In randomized trials (Section 5) no exit was ever observed from the prescribed initialization of Section 3.3. From randomly rotated feasible initializations, exits occurred on a fraction 10−5 to 2 × 10−4 of the instances at n = 3, the fraction growing as κ ranges over 102 to 103 . At n = 5 no exit was ever observed. 4 ′ ′ Strictly, their analysis prescribes the smaller step α′ /(2 ∑n i=1 wi λmax (Ri )) ≤ α /(2β ); every small-step comparison in this article runs at the larger α′ /(2β ′ ), which is conservative in their favor. 1

2

3

The second problem we consider is the invariant matrix projection problem of [BRT25], which generalizes Problem 1 and unifies certain problems of interest in quantum information [BRT25, Sec. 7]. Problem 2. Given a projective unitary representation U ∶= {Ui }i∈[n] ⊆ U(d) and R ∈ PD(d), 1 min B(R, S) ≡ min f2 (S), S∈P 2 S∈P

(4)

where P ∶= {P ∈ PD(d) ∶ [P, Ui ] = 0 for all Ui ∈ U} is the set of invariant positive matrices. Brahmachari, Rubboli, and Tomamichel [BRT25] introduced a modified fixed-point iteration to solve the above problem: ⎡

√

−1 ⎢ ⎛ St+1 ∶= St 2 ⎢⎢Φ

⎢ ⎝ ⎣

⎤2 ⎞⎥⎥ − 12 St RSt ⎥ St , ⎠⎥ ⎦ 1 2

1 2

S0 = I,

(5)

where Φ(X) ∶= n1 ∑i Ui XUi† is the twirling map associated with U . In [Afh25], the author showed that this update is unit-step BW gradient descent for the modified objective 1 2 B(R, Φ(S)) over the full positive definite manifold; Appendix E gives an arguably more elegant interpretation as unit-step BW gradient descent on the submanifold of invariant matrices, which we prove is totally geodesic. Brahmachari, Rubboli, and Tomamichel [BRT25] proved linear convergence of their algorithm, but—echoing the barycenter situation—with a worst-case dependence on the dimension. Remark 3. The invariant matrix projection problem generalizes the BW barycenter problem [BRT25]. We nonetheless analyze the two separately: the barycenter is the more widely used of the two problems, and treating it on its own keeps the constants and the notation simpler. Projection is all you need. In this article we present a Projected Bures–Wasserstein Gradient Descent (Projected BW-GD) algorithm: unit-step BW gradient descent composed, at each iteration, with an exact BW-metric projection onto the well-conditioned set [αI, βI]. We prove that this algorithm converges linearly at rate (1 − κ−3/2 ) for both Problems 1 and 2, with no dependence on the dimension. To our knowledge this is the first dimensionfree guarantee for any unit-step method on these problems; moreover, the iteration complexity κ3/2 log(1/ϵ) improves polynomially on the best known dimension-free guarantee, O(κ5/2 log(1/ϵ)) at step size 1/(2κ) [ACGS21]. Note that we do not claim to show the dimension-independent convergence of the FP algorithm. What we guarantee is that by merely incorporating an eigenvalue-clipping operation at every iteration, we can construct a Projected RGD algorithm, which reduces to the standard FP algorithm whenever the iterates land in [αI, βI], and converges at a dimension-independent rate to the barycenter S⋆ (see Theorem 10 for a precise statement). Note that the guarantee comes at no additional cost: clipping reuses an eigendecomposition the following iteration must perform in any case, so the projected and unprojected iterations cost essentially the same per step, except for an O(d) difference (see Section 3.3) which is dominated by the O(d3 ) cost of the eigendecompositions.

1.1

Related Works

The Bures–Wasserstein barycenter and Riemannian gradient descent. The BW barycenter has been extensively studied for its role in optimal transport, statistics, and 4

machine learning [AC11; CNR25; KSS21]; see also the references therein. Existence and consistency go back to [AC11; LL17]; statistical properties of empirical barycenters are studied in [PZ19; KSS21; LPRS22]. The FP algorithm of Álvarez-Esteban et al. [ÁDCM16], who also established its asymptotic convergence, was identified by Zemel and Panaretos [ZP19] as unit-step RGD on the BW manifold, with an independent convergence proof. Its non-asymptotic analysis has developed along two lines. Chewi et al. [CMRS20] gave the first linear rate at unit step size via a Polyak–Łojasiewicz (PL) inequality, with constants governed by the conditioning of the iterates; absent control of the smallest eigenvalue, these degrade as κd in the worst case. Altschuler et al. [ACGS21] then proved a dimension-free rate by establishing eigenvalue control along the trajectory, at the cost of the step-size restriction η ≤ 1/(2κ) and iteration complexity O(κ5/2 log(1/ϵ)). No step-size-independent eigenvalue control is possible as {S > 0 ∶ λmin (S) ≥ α} is not closed under generalized geodesics [ACGS21, Rem. 3], and Example 13 realizes the failure inside the feasible set, for two consecutive iterations. Recent work accelerates the same fixed-point iteration by Riemannian Anderson mixing [AEO26], with convergence established locally—on suitable BW balls around the solution—which leaves the global, unit-step regime of the practical algorithm without a guarantee. Our work is best read against this backdrop: rather than seeking finer control of the unprojected dynamics, we restore the invariance of the well-conditioned set by an exact, non-expansive projection. Invariant Matrix Projection. Brahmachari, Rubboli, and Tomamichel [BRT25] introduced the invariant matrix projection problem and a generalized fixed-point algorithm for it, which corresponds to unprojected unit-step RGD on the invariant submanifold (Appendix E). The framework encompasses, among others, the fidelity of asymmetry and of coherence in quantum resource theories, the max-conditional entropy and the order-1/2 sandwiched Rényi mutual information in quantum Shannon theory, and the geometric measure of entanglement of maximally correlated states [BRT25, Sec. 7]. Their linear convergence rate likewise carries a worst-case dependence on the dimension. Projected and safeguarded methods on the BW manifold. Projection steps with respect to the BW distance have appeared in two roles. As a constraint mechanism, Fan et al. [Fan+24] analyze projected BW-GD over a BW ball B(Q; r) ∶= {S ∶ B(S, Q) ≤ r2 }. Their setting differs from ours in both the feasible set and the role it plays: their BW ball encodes an external constraint imposed on the problem. Our feasible set is instead a spectral interval chosen for the analysis, and it enjoys three properties their ball does not: it provably contains the unconstrained optimum, projection onto it is closed-form eigenvalue clipping, and the projection is non-expansive. It is this combination that yields unit-step, dimension-free rates, which are unavailable in their setting. As a safeguard, one-sided eigenvalue clipping has been used to control iterates in BW-space analyses: Altschuler et al. [ACGS21, Prop. 3] prove that clipping from above is a BW contraction, and Lambert et al. [Lam+22] employ the same device in variational inference. Our Projection Lemma extends this in two ways: two-sided clipping to [α, β] is non-expansive, and it is the exact BW-metric projection onto [αI, βI]. As explained later, the extension requires a new argument.5 On the structural side, closedform BW projections onto sets of interest in quantum information—such as sets of bipartite matrices with a given marginal—were recently derived in [AT26]. That the one-sided clip is moreover the exact BW projection onto {S ≥ 0 ∶ S ≤ βI} is not stated in [ACGS21; Lam+22], whose analyses require only the contraction property; it follows readily from their argument in the total-manifold picture, since the lifted constraint set is a unitarily invariant convex ball and Euclidean projections onto such sets descend through the quotient. 5

5

Step size Chewi et al. [CMRS20] Altschuler et al. [ACGS21] Brahmachari et al. [BRT25]† This work†

η=1 η = 2κ1 ′ η=1 η=1

Iteration complexity

Dimension-free

κ -dependent (worst case) 5/2 O(κ′ log 1ϵ ) d-dependent (worst case) 3/2 κ′ log 1ϵ

✗ ✓ ✗ ✓

d

Table 1: Linear convergence guarantees for first-order methods for the BW barycenter (and, where indicated, the invariant matrix projection problem). Iteration complexity is the number of iterations to reach f (St ) − f (S⋆ ) ≤ ϵ(f (S0 ) − f (S⋆ )), up to absolute constants; d is the matrix dimension. Two condition numbers appear: κ = β/α from the spectral bounds of the individual inputs, and the refined κ′ = β ′ /α′ ≤ κ from the ′ ′ ensemble-average quantities of Eq. (7). The refined constant of [ACGS21] is ∑n i=1 wi λmax (Ri )/α ≥ κ , so writing their row in κ′ is conservative in their favor. κ and κ′ coincide in the worst case over ensembles. †

1.2

Also covers the invariant matrix projection problem.

Our contributions

The central object of our study is Projected BW-GD: unit-step RGD composed, at every iteration, with the BW projection onto [αI, βI], where α, β are read off the inputs (α = mini λmin (Ri ), β = maxi λmax (Ri ) for the barycenter; α = λmin (R), β = λmax (R) for the invariant projection) and are known to contain the optimum, S⋆ ∈ [αI, βI] [ACGS21; BRT25]. Our contributions are as follows. 1. A negative result: the well-conditioned set is not invariant at unit step size. Example 13 certifies that the well-conditioned set is not invariant under the unit step, and the construction extends to every n ≥ 3 and d ≥ 2. This sharpens the geodesic-convexity failure of [ACGS21, Rem. 3], whose configuration starts outside the well-conditioned set, into a statement about the algorithm itself. It pinpoints the quantity a fix must control— the spectral floor. 2. The Projection Lemma. We prove (Lemma 6) that for any positive definite P , clipping the eigenvalues of P to [α, β] is the unique minimizer of B(P, ⋅) over [αI, βI], and that this projection is non-expansive in the BW distance. The known one-sided result [ACGS21, Prop. 3] appeals to convexity of the spectral-norm ball and the 1-Lipschitz property of Euclidean projection; the two-sided result (specifically, clipping from below) does not follow from such an argument, and is supplied instead by the Alberti–Uhlmann variational formula [Alb83; Wat18], operator monotonicity of the matrix geometric mean [Bha09], and a pinching argument [Tom15]. We are not aware of another closed-form BW projection identified as such, and we expect it to serve as a constrained-optimization primitive on this manifold. 3. Dimension-free linear convergence at unit step size, with improved κdependence. Because the projection is exact and non-expansive, composing it with the unit gradient step (a) confines all iterates to [αI, βI], on which strong convexity and a PL inequality hold with dimension-free constants, and (b) preserves monotone descent. The result, stated as the unified theorem below and proved as Theorem 10, is the rate (1 − κ−3/2 ) for both problems at unit step size; Table 1 situates it among existing guarantees. An a posteriori certificate (Corollary 12) bounds the optimality gap of any iterate by κ3/2 2 B(St , K(St )), a quantity the iteration has already computed, allowing termination on the fly.

6

4. A unified square-root-manifold analysis, with structural by-products. We identify BW-GD on PD(d) with Euclidean GD on the square-root manifold GL(d) (Burer– Monteiro factored GD [BM03]). Euclidean strong convexity of f on [αI, βI] then lifts to a PL inequality for g = f ○ π with dimension-free constant κ−3/2 , and the descent lemma and the projection combine in three lines to give the rate (Section 4). This bypasses generalized geodesics and convexity along them, which earlier analyses required [AGS08; CMRS20; ACGS21]: their role is played by straight chords of the flat square-root space. Two byproducts follow. First, the fixed-point algorithm of Brahmachari, Rubboli, and Tomamichel [BRT25] is exactly unit-step RGD on the submanifold of invariant matrices, which we prove is totally geodesic (Appendix E). Second, the projected iteration is singular-value-clipped Euclidean GD upstairs (Appendix D), under which the barycenter update is an averaging of Procrustes-aligned square roots—the form in which the unit-step descent lemma is a variance identity (Appendix B). The following theorem summarizes our main convergence guarantee. Theorem 4 (Unified Linear Convergence of Projected BW-GD). Let f ∈ {f1 , f2 } be the objective function with unique optimum S⋆ ∈ [αI, βI] and κ = β/α. For any feasible initial point S0 ∈ [αI, βI] (additionally invariant, S0 ∈ P, when f = f2 ), the iterates {St }t≥0 of Projected BW-GD with unit step size satisfy f (St ) − f (S⋆ ) ≤ (1 −

1

t

) (f (S0 ) − f (S⋆ )). κ3/2

(6)

For the barycenter this improves to (1 − κ′−3/2 ) on the refined interval (Remark 11). Proof. See Section 4. Refined constants.

For the barycenter the condition number can be sharpened. Define

n √ 2 α′ ∶= [ ∑ wi λmin (Ri )]

n

and

i=1

β ′ ∶= λmax ( ∑ wi Ri ),

(7)

i=1

which satisfy α ≤ α′ ≤ β ′ ≤ β and still enclose the optimum, S⋆ ∈ [α′ I, β ′ I] (floor by [ACGS21, Thm. 6], ceiling by [BJL19, Thm. 9]). Projecting onto [α′ I, β ′ I] yields the identical guarantee with κ′ ∶= β ′ /α′ ≤ κ (Remark 11, proved in Appendix C); this requires additional ingredients, since the Ri need not lie in the refined interval. Since β ′ ≤ ∑ni=1 wi λmax (Ri ), our κ′ is no larger than the refined condition number of [ACGS21, Thm. 3], so the rates compare favorably instance-wise. We work with [αI, βI] in the main text, which is simpler and matches the conventions of the literature. Techniques. The analysis splits into an RGD part and a projection part. For the RGD part, the ingredients—the unit-step descent inequality [BJL19; BRT25] and the passage from Euclidean strong convexity to a PL inequality [CMRS20; Fan+24]—are available in the literature; we route them through the total-manifold identification, under which the Riemannian algorithm becomes a Euclidean gradient method, and the classical PL mechanism [Pol63; KNS16] drives the proof. The PL inequality must come from Euclidean strong convexity, for the non-negative curvature of the BW manifold renders the objective geodesically nonconvex [ACGS21, App. B.2]; the descent at the full unit step, by contrast, comes from that same curvature, expressed as a variance inequality on the square-root manifold (Proposition 14). Curvature thus takes away geodesic convexity with one hand and returns unit-step descent with the other. For the projection part, the one-sided argument of [ACGS21, Prop. 3] 7

rests on convexity, which fails for the floor; we establish the two-sided statement with tools from quantum information (Lemma 6). Since the objective is built from squared BW distances to references inside [αI, βI]—each fixed by the projection—non-expansiveness ensures the projection does not undo the progress of the RGD step, while confining the iterate to the set on which the PL inequality holds. Organization. Section 2 collects preliminaries. Section 3 develops the algorithm: the square-root-manifold identification (Section 3.1), the Projection Lemma (Section 3.2), and the algorithm itself (Section 3.3). Section 4 proves the convergence theorem and the a posteriori certificate, and Section 5 presents numerical experiments. The appendices contain the exit instance (Appendix A), the variance-inequality proof of the descent lemma (Appendix B), the refined-interval guarantee (Appendix C), the total-manifold implementation (Appendix D), and the identification of the algorithm of Brahmachari, Rubboli, and Tomamichel [BRT25] as BW-GD on the totally geodesic invariant submanifold (Appendix E).

2

Mathematical Preliminaries

For a positive integer n, [n] ∶= {1, . . . , n}. For d ≥ 2, L(d) denotes the d × d complex matrices and GL(d) the invertible ones. We equip L(d) with the (real) Hilbert–Schmidt (HS, a.k.a. Frobenius) inner product ⟨A, B⟩ ∶= RTr[A† B], where A† is the conjugate transpose.6 The real vector space of Hermitian matrices is Hr(d) ∶= {A ∈ L(d) ∶ A = A† }; within Hr(d) lies the closed cone PSD(d) of positive semidefinite matrices, and its open subset PD(d) of positive definite matrices PSD(d) = {P ∈ Hr(d) ∶ λi (P ) ≥ 0},

and

PD(d) = {P ∈ Hr(d) ∶ λi (P ) > 0},

(8)

where λmin (P ) ≡ λ1 (P ) ≤ . . . ≤ λd (P ) ≡ λmax (P ) are the eigenvalues of P ordered nondecreasingly. U(d) ∶= {U ∈ L(d) ∶ U † U = U U † = Id } denotes the unitary matrices, where Id ≡ I denotes the d × d identity matrix. For A, B ∈ Hr(d) we write A > B (A ≥ B) when A − B is positive (semi)definite. For A, B ∈ PD(d), define the matrix geometric mean as 1 2

A#B ∶= A 2 (A− 2 BA− 2 ) A 2 , 1

1

1

1

(9)

equivalently the unique T ∈ PD(d) solving the Riccati equation T A−1 T = B; it is symmetric, √ 1 A#B = B#A, operator monotone in each argument [Bha09], and satisfies A#(cI) = c A 2 for c > 0. The Bures–Wasserstein manifold. PD(d) is an open subset of the real vector space Hr(d) and hence a smooth manifold, with tangent space TP Pd(d) ≅ Hr(d) at every P . We equip each tangent space with the BW metric tensor 1 1 BW gBW = Tr[XLP (Y )] = Tr[LP (X)Y ], P (X, Y ) ≡ ⟨X, Y ⟩P 2 2

(10)

where LP (X) is the unique Hermitian solution of the Lyapunov equation X = LP (X)P + P LP (X) [BJL19; Bha09]. This metric induces the squared BW distance B(P, Q) defined in 6

In the settings common in machine learning and statistics the matrices are real (e.g. covariance matrices), in which case A† reduces to the transpose A⊺ and the real part in the inner product is redundant. All results of this article hold verbatim in the real case.

8

the introduction; we refer to [BJL19; MMP18] for proofs of the facts collected here. The Riemannian exponential map is ExpBW P [X] ∶= [I + LP (X)] P [I + LP (X)] ;

dom(ExpBW P ) = {X ∶ I + LP (X) ∈ PD(d)}. (11) A more illustrative construction of the BW manifold realizes it as a quotient of the manifold of invertible matrices. Being an open subset of L(d), the set GL(d) is a smooth manifold, which we equip with the Frobenius inner product and distance ∥A − B∥F . Consider on GL(d) the unitary equivalence A ∼ B ⇐⇒ A = BU for some U ∈ U(d), together with the square map π ∶ GL(d) → PD(d), π(A) ∶= AA† , (12) which satisfies π(A) = π(B) if and only if A ∼ B; thus PD(d) ≡ GL(d)/U(d). The equivalence class of P ∈ PD(d), its fibre, is its set of square roots π −1 [P ] = {A ∈ GL(d) ∶ AA† = P } = {P 2 U ∶ U ∈ U(d)}, 1

(13)

1

with P 2 the principal square root.7 We extend the above notation to sets as well: for S ⊆ PD(d), we define π −1 [S] ∶= {A ∶ π[A] ∈ S} ⊆ GL(d). We call GL(d) the total (or square-root) manifold and PD(d) the base manifold. Since right multiplication by a unitary preserves the Frobenius distance, the distance between fibres is well defined, and it is precisely the BW distance: B(P, Q) =

1

min

A∈π −1 [P ], B∈π −1 [Q]

1

∥A − B∥2F = min ∥P 2 U − Q 2 ∥2F ;

(14)

U ∈U(d)

a pair (A, B) attaining the minimum is an optimal pairing for (P, Q). Expanding the right-hand side and using maxU ∈U(d) ReTr[M U ] = ∥M ∥1 recovers the closed form B(P, Q) = Tr[P + Q] − 2F(P, Q). The construction also recovers the BW metric: it is (up to scaling) the unique metric on PD(d) compatible with the Frobenius metric on GL(d) under π [BJL19]. For full-rank P, Q (always the case in this work), the optimal unitary in Eq. (14) 1 1 is U = Pol(P 2 Q 2 ), the unitary polar factor (see [BJL19], or [AF25, Lem. C.1]); moreover 1 1 1 1 (P 2 U, Q 2 ) and (P 2 , Q 2 U † ) are optimal pairings for (P, Q). Further properties of the BW distance. Write BR (⋅) ≡ B(R, ⋅) for the squared BW distance to a fixed R ∈ PD(d). Its Euclidean gradient over PD(d) is ∇BR (S) = I − S −1 #R,

(15)

so the gradient of the barycenter functional is ∇f1 (S) = 12 (I − ∑ni=1 wi [S −1 #Ri ]) [ACGS21; BJL19]. The squared BW distance is strongly convex (SC), and in a form that √does not λ

(R)

min require the reference and the argument to share an interval: for every b > 0, BR is 2b 3/2 SC on {S ∈ PD(d) ∶ S ≤ bI} [BJL18; BRT25], that is, √ 1 λmin (R) BR (Q) ≥ BR (S) + ⟨∇BR (S), Q − S⟩ + ∥Q − S∥2F (16) 3/2 2 2b for all Q, S ≤ bI (Lemma 16 records the Hessian form we use in Appendix C). Specializing to R ∈ [αI, βI] and b = β gives λmin (R) ≥ α, so both f1 and f2 are µ-SC on the corresponding [αI, βI], where we set once and for all 1

1

α2 1 ⎛ α2 ⎞ = 3/2 , µ ∶= 3/2 2 ⎝ 2β ⎠ 4β

so that

7

4µα = κ−3/2 .

(17)

The quantum information audience would recognize these square roots to be related to the purifications of P under the vec isomorphism.

9

The BW distance satisfies the data processing inequality (DPI) [Uhl11; Wat18]: for any quantum channel (completely positive trace-preserving map) Λ, B(P, Q) ≥ B(Λ(P ), Λ(Q)).

(18)

We will use the DPI for pinching channels [Tom13]: given mutually orthogonal projectors E ∶= {Ek }k∈[m] with Ek El = δkl Ek and ∑k Ek = Id , the associated pinching channel is ΛE (X) ∶= ∑k∈[m] Ek XEk . Preliminaries for the invariant matrix projection problem. Let U = {Ui }ni=1 be a projective unitary representation. Its commutant is the unital subalgebra comm(U) ∶= {X ∈ L(d) ∶ [X, Ui ] = 0 for all i} [Wat18], a vector subspace of L(d). Define the invariant sets P ∶= PD(d) ∩ comm(U), H ∶= Hr(d) ∩ comm(U), and G ∶= GL(d) ∩ comm(U). The twirl Φ(X) ∶= n1 ∑ni=1 Ui XUi† is a unital quantum channel and is the Hilbert–Schmidt orthogonal projection onto comm(U): Φ(X) ∈ comm(U)

and

⟨Φ(X), X − Φ(X)⟩ = 0

for all X ∈ L(d).

(19)

An important geometric distinction: for invariant P ∈ P, the tangent space to the full manifold remains all of Hr(d), whereas the tangent space to the invariant submanifold is exactly H: TP PD(d) ≅ Hr(d) and TP P ≅ H for all P ∈ P. (20)

3

The Algorithm: Descent

Projected Bures–Wasserstein Gradient

This section develops the two components of the algorithm and assembles them. Section 3.1 recalls BW gradient descent and its identification with Euclidean GD on the square-root manifold; Section 3.2 proves the Projection Lemma; Section 3.3 states the algorithms.

3.1

BW gradient descent and the square-root manifold

Let f ∶ PD(d) → R be differentiable with Euclidean gradient ∇f (P ) ∈ Hr(d), defined by Df (P )[X] = ⟨∇f (P ), X⟩. The BW (Riemannian) gradient ∇f (P ) is the dual of the = Df (P )[X] = ⟨∇f (P ), X⟩ for all derivative with respect to the BW metric, ⟨∇f (P ), X⟩BW P X ∈ TP Pd(d). A comparison readily yields a relation between the two gradients: 1 ∇f (P ) = LP (∇f (P )) 2

and ∇f (P ) = 2L−1 P (∇f (P )) = 2 (P ∇f (P ) + ∇f (P )P ) ,

(21)

using L−1 P (X) = XP + P X. Riemannian gradient descent with step size η > 0 is therefore P ↦ ExpBW P [−η ∇f (P )] = [I − 2η∇f (P )] P [I − 2η∇f (P )] ,

(22)

defined whenever the step lies in the domain of the exponential map, i.e. I − 2η∇f (P ) > 0.

(23)

For both of our objectives (with twirled gradient for f2 ) the unit step η = 1 always satisfies this: I − 2∇f1 (S) = ∑ni=1 wi [S −1 #Ri ] and I − 2Φ(∇f2 (S)) = Φ(S −1 #R) are positive definite, being a weighted sum, respectively a twirl, of geometric means of positive definite matrices.

10

It is instructive to view this update in the square-root manifold. Let g ∶= f ○ π denote the lifted function, with dom(g) = GL(d); by the chain rule ∇g(A) = 2∇f (AA† )A. For any A ∈ π −1 [P ], a Euclidean GD step of g satisfies A ↦ A − η∇g(A) = (I − 2η∇f (P ))A,

(24)

whose image under π is exactly the RGD update above. Thus BW gradient descent on PD(d) is Euclidean gradient descent on GL(d), equivalently Burer–Monteiro factored GD [BM03; BM05]; this observation appears in [ZHVZ25, Corollary 5.2] and [MLR23], and we exploit it throughout. The next result expresses the BW length of a gradient step through the lifted gradient; the underlying reason it is a one-line computation is that the fibres of π are orbits of an isometry group, so a straight line upstairs that starts at a well-chosen square root remains optimally paired with its starting fibre. For a more formal version of the above statement, see [BJL19, Theorems 3 and 4]. Lemma 5 (Step length via the lifted gradient). Let S ∈ PD(d) and A ∈ π −1 [S]. For any X ∈ Hr(d) with M ∶= I − X > 0, define A′ ∶= M A and S ′ ∶= π(A′ ) = M SM . Then B(S, S ′ ) = ∥A − A′ ∥2F = Tr[X 2 S],

(25)

independently of the choice of A ∈ π −1 [S]. In particular, let f ∶ PD(d) → R be differentiable, g ∶= f ○ π, and let η > 0 satisfy Eq. (23) at S. Taking X = 2η∇f (S) gives A′ = A − η∇g(A) and S ′ = ExpBW S [−η ∇f (S)], whence B(S, S ′ ) = ∥A − A′ ∥2F = η 2 ∥∇g(A)∥2F .

(26)

Proof. Since A − A′ = XA, we have ∥A − A′ ∥2F = Tr[A† X 2 A] = Tr[X 2 S], independent of the † choice of A. Moreover ∥A − A′ ∥2F = Tr[S + S ′ ] − 2R⟨A, A′ ⟩, using S = AA† and S ′ = A′ A′ , so ′ ′ ′ it suffices to show R⟨A, A ⟩ = F(S, S ), i.e. that (A, A ) is an optimal pairing. Indeed, ⟨A, A′ ⟩ = Tr[A† M A] = Tr[M S] = F(S, S ′ ),

(27)

since M = S −1 #S ′ (via the matrix Riccati equation [Bha09]) and F(P, Q) = Tr[P ⋅ (P −1 #Q)] for P, Q ∈ PD(d). For the particular case, the chain rule gives ∇g(A) = 2∇f (S)A, so A′ = M A = A − η∇g(A) with M = I − 2η∇f (S) > 0 by Eq. (23). Consequently, S ′ = 2 2 2 2 2 M SM = ExpBW S [−η ∇f (S)], and Tr[X S] = 4η Tr[∇f (S) S] = η ∥∇g(A)∥F , again by the chain rule. Thus the BW distance between consecutive BW-GD iterates equals the Euclidean length of the lifted gradient for any square root. This lets us monitor progress and certify output purely through gradient norms (Corollary 12).

3.2

Eigenvalue clipping is the non-expansive BW projection

We now prove our central technical contribution—the Projection Lemma: clipping the eigenvalues of a positive definite matrix to [α, β] yields the closed-form BW-metric projection onto the compact set [αI, βI], and this projection is non-expansive (1-Lipschitz). The feasible set here is a spectral one—constraining the whole spectrum to an interval, rather than the trace, the norm, or a distance to a reference—and it is this kind of constraint that arises whenever the conditioning of an iterate must be controlled. We begin by formally defining the eigenvalue-clipping operation. Throughout 0 < α ≤ β < ∞, so that clipping maps PD(d)

11

into itself. We define clipα,β ∶ R → [α, β] as clipα,β (x) ∶= min{max{α, x}, β}, extended to Hermitian matrices spectrally: d

clipα,β (H) ∶= ∑ clipα,β (λi )vi vi† ,

(28)

i=1

where H = ∑di=1 λi vi vi† is the eigendecomposition of H ∈ Hr(d). For non-normal matrices we extend the operation through the singular value decomposition, denoted clipsv α,β . Before the formal proof, it is worth seeing why the first statement should be true. Aligning toward the eigenbasis of P can only help: pinching in that basis fixes P , preserves [αI, βI], and never increases BW distances, so the minimizer may as well commute with P —and once it does, the matrix problem degenerates into d independent scalar problems, each solved by clipping. The formalization of the above statement into a proof is done by using the data processing inequality and strict convexity of the squared BW distance. For the non-expansive property, we use the Alberti–Uhlmann characterization of BW distance and operator monotonicity of the matrix geometric mean. Lemma 6 (Projection Lemma). For any positive definite matrix P , the unique BW projection onto the geodesically convex set [αI, βI] (for 0 < α ≤ β < ∞) is given by eigenvalueclipping: (29) clipα,β (P ) = argmin B(P, Q). Q∈[αI,βI]

Moreover, the projection clipα,β is non-expansive (1-Lipschitz) with respect to the BW distance. For all P, Q ∈ PD(d), B(clipα,β (P ), clipα,β (Q)) ≤ B(P, Q).

(30)

Proof. We first prove that the projection is given by clipping. To this end, observe that the squared BW distance B(P, ⋅) is strictly convex over positive definite matrices [BJL18], and [αI, βI] is convex and compact, so the projection exists and is unique. We claim the projection Q⋆ commutes with P . Let Λ be the pinching channel defined by the eigenprojectors of P , so Λ(P ) = P . For any Q ∈ [αI, βI] we have Λ(Q) ∈ [αI, βI] (Λ is positive and unital, hence order-preserving), and the DPI gives B(P, Λ(Q)) = B(Λ(P ), Λ(Q)) ≤ B(P, Q). Applied to Q⋆ , this shows Λ(Q⋆ ) is also a minimizer; uniqueness forces Λ(Q⋆ ) = Q⋆ , i.e. [Q⋆ , P ] = 0. The problem therefore reduces, in a common eigenbasis, to the scalar problem d

min

Q∈[αI,βI], [Q,P ]=0

1

1

2

B(P, Q) = min ∑ (λi2 − ωi2 ) , ωi ∈[α,β] i=1

(31)

with (λi )i∈[d] the eigenvalues of P , which is minimized by clipping each λi to [α, β], resulting in clipα,β (P ). We now prove the projection is 1-Lipschitz. Begin with the Alberti–Uhlmann variational characterization [Alb83; BJL19]: for P, Q > 0, B(P, Q) = max (Tr[P + Q] − Tr[P X] − Tr[QX −1 ]) = max (⟨P, I − X⟩ + ⟨Q, I − X −1 ⟩) , (32) X>0

X>0

with optimizer X = P −1 #Q. Write P̂ ≡ clipα,β (P ), Q̂ ≡ clipα,β (Q), and let Y ∶= P̂ −1 #Q̂ be the optimizer for the pair (P̂ , Q̂). Since Y > 0 is a feasible point, we have Tr[P + Q] − Tr[P Y ] − Tr[QY −1 ] ≤ B(P, Q). It therefore suffices to prove the following: ?

B(P̂ , Q̂) = ⟨P̂ , I − Y ⟩ + ⟨Q̂, I − Y −1 ⟩ ≤ ⟨P, I − Y ⟩ + ⟨Q, I − Y −1 ⟩ ≤ B(P, Q), 12

(33)

?

or equivalently we want to show that ⟨P̂ − P, Y − I⟩ + ⟨Q̂ − Q, Y −1 − I⟩ ≥ 0. Via the spectral decomposition P = ∑i λi vi vi† , define the gap matrices ∆βP = ∑ (λi − β)vi vi† ≥ 0,

∆αP = ∑ (α − λi )vi vi† ≥ 0,

(34)

λi >β

λi <α

so that P̂ − P = ∆αP − ∆βP , and define ∆αQ , ∆βQ analogously. Then the inequality in question (of Eq. (33)) becomes ?

Tr[∆αP (Y − I)] + Tr[∆βP (I − Y )] + Tr[∆αQ (Y −1 − I)] + Tr[∆βQ (I − Y −1 )] ≥ 0.

(35)

We now show each term is nonnegative via operator monotonicity of the geometric mean: R1 ≥ S1 and R2 ≥ S2 ⇒ R1 #R2 ≥ S1 #S2 [Bha09]. For the first term: Q̂ ≥ αI implies 1 1 Y = P̂ −1 #Q̂ ≥ P̂ −1 #(αI) = α 2 P̂ − 2 . If vi is an eigenvector of P with eigenvalue λi ≤ α, then vi is an eigenvector of P̂ with eigenvalue α, whence √ 1 1 ⟨vi , (Y − I)vi ⟩ ≥ α 2 ⟨vi , P̂ − 2 vi ⟩ − 1 = α/α − 1 = 0, (36) and therefore Tr[∆αP (Y − I)] = ∑i∶λi <α (α − λi )⟨vi , (Y − I)vi ⟩ ≥ 0. The second term is nonnegative by the same argument starting from Q̂ ≤ βI. For the last two, observe Y −1 = (P̂ −1 #Q̂)−1 = P̂ #Q̂−1 , and apply the same argument with the roles of (P, Q) exchanged. The same operation lifts into the exact Frobenius projection in the total manifold: √ √ [A] = clipsv α, β

argmin

B∈π −1 [[αI,βI]]

∥A − B∥F ,

(37)

see Corollary 18 in Appendix D. Unlike its analog in the base BW manifold, the totalmanifold Frobenius projection—whose feasible set π −1 [[αI, βI]] is non-convex—can be expansive (see Remark 19). Remark 7 (Extension to semidefinite matrices). Lemma 6 is stated for positive definite matrices; both assertions extend to all P, Q ∈ PSD(d) by continuity. The map clipα,β and the squared BW distance are jointly continuous, so the non-expansiveness inequality passes to the limit directly. For the projection property, take Pε ∶= P + εI for ε ↘ 0: minimality of clipα,β (Pε ), i.e. B(Pε , clipα,β (Pε )) ≤ B(Pε , Q) for every Q ∈ [αI, βI], survives the limit ε ↘ 0, so clipα,β (P )—which lifts the zero eigenvalues of P to α—remains a BW projection of P onto [αI, βI]. We caution that strict convexity of B(P, ⋅) can fail for singular P along directions supported on ker P , so uniqueness of the projection can fail in the semidefinite case; we never need it, as the algorithm only projects positive definite iterates.

3.3

Projected Bures–Wasserstein gradient descent

A single iteration of Projected BW-GD composes the unconstrained gradient step with the projection. Throughout we take η = 1; the unconstrained unit step is then exactly the fixed-point map of [ÁDCM16] (respectively [BRT25]), which we denote by K: St+1 ∶= clipα,β [K(St )],

K(St ) ∶= ExpBW St [−∇f (St )] = [I − 2∇f (St )] St [I − 2∇f (St )]. (38)

Both instantiations are collected in Algorithm 1; the descent and trapping properties that make the composition work are established in Section 4.

13

Algorithm 1 Projected BW Gradient Descent (barycenter and invariant matrix projection) Require: Barycenter (f = f1 ): ensemble R = (Ri )ni=1 and weight vector w; invariant projection (f = f2 ): R ∈ PD(d) and unitary representation U = {Ui }ni=1 ; iteration count T. ⎧ ⎪ ⎪( mini λmin (Ri ), maxi λmax (Ri )) ∶ f = f1 1: (α, β) ← ⎨ ⎪ ∶ f = f2 ⎪ ⎩(λmin (R), λmax (R)) α+β 2: Initialize S0 ∈ [αI, βI], additionally S0 ∈ P when f = f2 (e.g., S0 ← 2 I) 3: for t = 0, 1, . . . , T − 1 do √ 2 ⎧ 1 1⎞ ⎪ − 12 ⎛ n −1 ⎪ ⎪ 2 2 ⎪ St St 2 ∶ f = f1 ∑i=1 wi St Ri St ⎪ ⎪ ⎪ ⎝ ⎠ ⎪ 4: K(St ) ← ⎨ ▷ Step 1: Unit-sized BW-GD √ 2 ⎪ 1 1 ⎞⎞ 1 1 ⎛ ⎪ ⎛ − − ⎪ ⎪ ⎪ St2 RSt2 St 2 St 2 Φ ∶ f = f2 ⎪ ⎪ ⎪ ⎝ ⎝ ⎠ ⎠ ⎩ Step. 5: St+1 ← clipα,β [K(St )] ▷ Step 2: Projection onto [αI, βI] via Clipping. 6: end for 7: return ST Section 4 says how large the iteration count T must be, and Corollary 12 gives an a posteriori stopping rule that dispenses with fixing T in advance. The projection is free. We note that the projection, via eigenvalue clipping, is essentially free for both problems—it requires nothing more than what the unprojected variant already computes. Indeed the clipping operation, which constructs St+1 from K(St ), requires an ±1/2 eigendecomposition of K(St ). However, the very next iteration opens by forming St+1 in order to evaluate K(St+1 ), which also requires the eigendecomposition of St+1 . Since St+1 = clipα,β [K(St )] is obtained from K(St ) by acting on its eigenvalues alone, the two share an eigenbasis, and the single decomposition serves both: the projection merely does, one step early, what the unprojected iteration would perform at the start of the very next step. The clip itself then costs O(d), a comparison per eigenvalue. Counting spectral decompositions, √ 1/2

1/2

each of cost O(d3 ): an iteration performs n of them for the geometric means St Ri St and one for the new iterate, whether or not it is projected—so both variants cost n + 1 per step, and the clip contributes only the O(d) comparisons. Hence the iteration counts reported in Section 5 may be read as running times. For clarity of exposition, Algorithm 1 is stated without this reuse; the implementation accompanying the article realizes it. The iteration on the total manifold. The base–total equivalence also lets the whole iteration run upstairs: from any A0 ∈ π −1 [S0 ], singular-value-clipped Euclidean GD of g reproduces the projected iteration exactly, returning to PD(d) only upon convergence (Appendix D). We record this not for its computational content but for its structure: upstairs, the update is an averaging of Procrustes-aligned square roots—the form in which the unitstep descent lemma is a variance identity and Euclidean strong convexity becomes the PL inequality that drives Section 4.

14

4

Dimension-Independent Linear Convergence at Unit StepSize

We may now prove the unified convergence theorem, and to this end, we clarify the notation. Throughout this section, f ∈ {f1 , f2 }, with (α, β) read off the inputs as in Algorithm 1 and ⎧ ⎪ ⎪∇g ˆ ∇g ∶= ⎨ ⎪ ⎪ ⎩Φ ○ ∇g

∶ f = f1 ∶ f = f2 .

(39)

Moreover, the domain of interest is [αI, βI] (and π −1 [[αI, βI]]) when f = f1 and [αI, βI]∩P (and π −1 [[αI, βI]]∩G) when f = f2 . In both cases, f is µ-SC on [αI, βI] with µ as in Eq. (17). The optimum S⋆ exists, is unique, and lies in the corresponding domain of interest: existence and uniqueness are classical for the barycenter [AC11] and follow from strong convexity on the compact feasible region for the invariant projection, while the spectral inclusion S⋆ ∈ [αI, βI] is established in [ACGS21; BRT25]. The first ingredient is a PL inequality, obtained by transferring strong convexity through the lift of Section 3.1, a fact previously noted in [Fan+24, Lemma 3.4]; we state it in a form covering both of our problems at once, with the twirl Φ taken to be the identity map for the barycenter problem (U = {I}, so comm(U) = L(d), P = PD(d), and G = GL(d)). Lemma 8 (PL inequality for the lifted problem). Let 0 < a ≤ b < ∞, U be a projective unitary representation with twirl Φ, and f be µ0 -strongly convex on [aI, bI] with S⋆ = argminS∈P∩[aI,bI] f (S). Then g ∶= f ○ π satisfies, for every A ∈ π −1 [[aI, bI]] ∩ G, g(A) − g(A⋆ ) ≤

1 ∥Φ(∇g(A))∥2F , 2(4µ0 a)

(40)

where A⋆ ∈ π −1 [S⋆ ] ∩ G. In particular, for the barycenter functional f1 the twirl is trivial (U = {I}), and g1 ∶= f1 ○ π satisfies g1 (A) − g1 (A⋆ ) ≤ 2(4µ1 0 a) ∥∇g1 (A)∥2F on π −1 [[aI, bI]]. Proof. Let S ∶= π(A) ∈ P ∩ [aI, bI]. For any P ∈ P ∩ [aI, bI], strong convexity gives f (P ) ≥ f (S) + ⟨∇f (S), P − S⟩ +

µ0 ∥P − S∥2F . 2

(41)

Since P − S ∈ H and Φ is the self-adjoint idempotent HS-projection onto comm(U), we may replace ∇f (S) by G ∶= Φ(∇f (S)) in the inner product. Writing Z ∶= P − S and completing the square, µ0 µ0 G 2 1 1 ⟨G, Z⟩ + ∥Z∥2F = ∥Z + ∥ − ∥G∥2F ≥ − ∥G∥2F . (42) 2 2 µ0 F 2µ0 2µ0 Substituting into Eq. (41) and choosing P = S⋆ , f (S) − f (S⋆ ) ≤

1 ∥Φ(∇f (S))∥2F . 2µ0

(43)

To lift, note Φ(∇g(A)) = 2Φ(∇f (S))A, using Φ(XA) = Φ(X)A for A ∈ G; hence ∥Φ(∇g(A))∥2F = 4Tr(Φ(∇f (S))2 S) ≥ 4a∥Φ(∇f (S))∥2F , since S ≥ aI. Combining with Eq. (43) and g(A) = f (S), g(A⋆ ) = f (S⋆ ) gives the claim.

15

(44)

The special case of Lemma 8 for f = f2 appears in [BRT25, Eq. 58]. We broadly follow the proof idea used there and in [KNS16]. The lemma as stated is strictly more general: it covers both problems at once and applies to any strongly convex f , with the setting of [BRT25] recovered on choosing f = f2 . In Appendix E we explain that the twirled gradient, rather than the gradient, appears because (i) in the base-manifold perspective, it is the BW-orthogonal projection of ∇f2 (P ) onto the invariant tangent subspace H ≡ TP P ⊆ TP Pd(d) and (ii) in the total-manifold perspective, it is the (Hilbert–Schmidt) orthogonal projection of the gradient ∇g(A) to the invariant subspace comm(U) (of L(d)). The second ingredient is the descent half of the argument, supplied by the following proposition; combined with the PL inequality above, it yields the linear rate. Proposition 9 (Projected descent lemma). Let f ∈ {f1 , f2 } and g ∶= f ○ π be the lifted functional. Let the current square-root iterate be ⎧ ⎪ ⎪π −1 [[αI, βI]] At ∈ ⎨ −1 ⎪ ⎪ ⎩π [[αI, βI]] ∩ G

∶ f = f1 , ∶ f = f2

(45)

ˆ √ √ [At − ∇g(A St = π(At ), and set At+1 = clipsv t )] and St+1 = π(At+1 ) = clipα,β [K(St )]. α, β Then: (i) (Feasibility.) St+1 ∈ [αI, βI] and At+1 ∈ π −1 [[αI, βI]]; when f = f2 , additionally St+1 ∈ P and At+1 ∈ G. (ii) (Descent.) f (St+1 ) ≤ f (St ) − 12 B(St , K(St )), equivalently, upstairs, 1 2 ˆ g(At+1 ) ≤ g(At ) − ∥∇g(A t )∥F . 2

(46)

Proof. Assume At satisfies Eq. (45) and St ∶= π(At ). We first verify that feasibility is preserved. Membership St+1 ∈ [αI, βI] holds by construction of the clipping map. For f2 , invariance is also preserved. Upstairs, At ∈ G Ô⇒ At − Φ(∇g2 (At )) ∈ G.

(47)

Recall that G ∶= comm(U)∩GL(d). The inclusion (of the unclipped next iterate) in comm(U) follows since both terms lie in comm(U ) (Φ maps into comm(U)) and comm(U) is a vector subspace of L(d). For invertibility, note that Φ(∇g2 (At )) = 2Φ(∇f2 (St ))At , using Φ(XA) = Φ(X)A for A ∈ G; hence At − Φ(∇g2 (At )) = [I − 2Φ(∇f2 (St ))]At , and I − 2Φ(∇f2 (St )) = Φ(St−1 #R) > 0 (see Eqs. (96) and (97)) while At is invertible by assumption. On the base manifold, the unconstrained step K(St ) = Φ(St−1 #R) St Φ(St−1 #R) is a product of elements of comm(U ), and comm(U ) is an algebra; hence K(St ) ∈ P. Since clipα,β acts spectrally, we conclude St+1 ∈ [αI, βI] ∩ P. Finally, to see that At+1 ∈ G, write A′t ∶= At − Φ(∇g2 (At )) ∈ G and note that singular value clipping preserves the polar factor (it modifies only the singular values), so that At+1 = Pol(A′t ) ⋅ clip√α,√β (∣A′t ∣) (48) = (A′t ∣A′t ∣−1 ) ⋅ clip√α,√β (∣A′t ∣) ∈ comm(U), where the last inclusion follows from the facts that every factor in the product is an ele′ 12 ment of comm(U)—for A′t by the previous step; for ∣A′t ∣ = (A′† t At ) and its clipped version because, for Hermitian Z ∈ comm(U) and any function h applied spectrally, Ui h(Z) Ui† = h(Ui ZUi† ) = h(Z), so h(Z) ∈ comm(U); and for ∣A′t ∣−1 since the inverse of an invertible element of comm(U) lies in comm(U )—and that comm(U ) is an algebra and hence it is closed 16

√ √ under products. As the clipped singular values lie in [ α, β], At+1 is invertible; hence At+1 ∈ comm(U) ∩ GL(d) = G. Because the reference matrices—the ensemble (Ri )i∈[n] for f1 , the target R for f2 —lie in [αI, βI] and are fixed by clipα,β , the non-expansiveness of the BW projection guarantees that the objective does not increase upon projecting the unconstrained step K(St ): f (St+1 ) = f (clipα,β [K(St )]) ≤ f (K(St )).

(49)

Furthermore, the objective value at the unconstrained step satisfies the descent inequality f (K(St )) ≤ f (St ) − 12 B(St , K(St )), established for f1 in [BJL19, Theorem 10] (an alternative proof of which is provided as Proposition 14 in the Appendix) and for f2 in [BRT25, Eq. 37]. Combining these two inequalities yields the descent lemma for the projected sequence in the base manifold: 1 f (St+1 ) ≤ f (St ) − B(St , K(St )). (50) 2 To establish the equivalence in the square-root manifold, recall that f ○ π = g, yielding f (St ) = g(At ); moreover, since eigenvalue clipping on the base manifold is equivalent √ √ = clip to singular value clipping on the square-root manifold, π ○ clipsv α,β ○ π (Corolα, β lary 18), we have π(At+1 ) = St+1 and hence f (St+1 ) = g(At+1 ). It remains to establish that 2 ˆ B(St , K(St )) = ∥∇g(A t )∥F , which follows on applying Lemma 5 with ⎧ ⎪ −1 ⎪2∇f1 (St ) ˆ X = ∇g(At )At = ⎨ ⎪ ⎪ ⎩2Φ(∇f2 (St ))

∶ f = f1 . ∶ f = f2

(51)

The choice of X is admissible for Lemma 5: it is Hermitian, and ⎧ ⎪ ⎪∑n wi [St−1 #Ri ] I − X = ⎨ i=1−1 ⎪ ⎪ ⎩Φ(St #R)

∶ f = f1 ∶ f = f2

(52)

is positive definite in both cases, being a weighted sum (respectively, the twirl) of geometric ˆ means of positive definite matrices. Since ∇g(A t ) = XAt , the unconstrained Euclidean step −1 ˆ satisfies At − ∇g(At ) = (I − X)At ∈ π [K(St )], and Eq. (25) yields 2 ˆ B(St , K(St )) = ∥XAt ∥2F = ∥∇g(A t )∥F .

(53)

Combining the above, the descent inequality is equivalently expressed as g(At+1 ) ≤ g(At ) − 1 ˆ 2 2 ∥∇g(At )∥F . This concludes the proof. We now prove that the projected BW-GD algorithm achieves dimension-independent linear convergence at unit step-size. Theorem 10 (Dimension-independent linear convergence). Let f ∈ {f1 , f2 } and S0 ∈ [αI, βI] (additionally, S0 ∈ P if f = f2 ). The iterates of Algorithm 1 satisfy f (St ) − f (S⋆ ) ≤ (1 −

1

t

) (f (S0 ) − f (S⋆ )). κ3/2

(54)

Proof. The proof is obtained through a combination of Proposition 9 (the projected descent lemma) and Lemma 8 (PL inequality from strong convexity). Both are available at every step, as we now check by induction. Feasibility propagates along the iteration by Proposition 9(i): At ∈ π −1 [St ] and St ∈ [αI, βI]

⇒

At+1 ∈ π −1 [St+1 ] and St+1 ∈ [αI, βI], 17

(55)

and, when f = f2 , so does invariance: At , St ∈ comm(U ) ⇒ At+1 , St+1 ∈ comm(U ). The hypotheses hold at t = 0, since S0 ∈ [αI, βI] and A0 ∈ π −1 [S0 ] by assumption; and for f = f2 1

we may take A0 = S02 , which lies in comm(U) because S0 ∈ P and comm(U ) is closed under the functional calculus of its elements (see the proof of Proposition 9). Induction therefore gives St ∈ [αI, βI] and At ∈ π −1 [[αI, βI]] for every t ≥ 0, together with At , St ∈ comm(U) when f = f2 . Thus we may invoke both the projected descent lemma and the PL inequality at every step. The descent lemma gives 1 2 ˆ g(At+1 ) ≤ g(At ) − ∥∇g(A t )∥F 2

⇐⇒

1 2 ˆ εt+1 ≤ εt − ∥∇g(A t )∥F , 2

(56)

where the equivalence is obtained by subtracting g(A⋆ ) = f (S⋆ ) from both sides and defining εt ∶= g(At ) − g(A⋆ ). Since f is µ-SC over [αI, βI], we invoke the PL inequality (Lemma 8, with [aI, bI] = [αI, βI], µ0 = µ and a = α) 1 2 ˆ 4µαεt ≤ ∥∇g(A t )∥F . 2

(57) 3/2

Combining the above two expressions, and using the relation 4µα = ( αβ )

= κ−3/2 , yields

εt+1 ≤ (1 − κ−3/2 )εt . Unrolling the inequality over t iterations we get εt ≤ (1 − κ−3/2 )t ε0 . Substituting εt = f (St ) − f (S⋆ ) and ε0 = f (S0 ) − f (S⋆ ), we get the required result, which is manifestly dimension-independent: f (St ) − f (S⋆ ) ≤ (1 −

1

t

3 ) (f (S0 ) − f (S⋆ )).

(58)

κ2

This concludes the proof. Remark 11 (Sharper constants on the refined interval). For the barycenter problem the constants can be tightened: projecting onto [α′ I, β ′ I]—the refined bounds of Eq. (7), which also enclose the optimum [ACGS21; BJL19]—yields the verbatim analogue of Theorem 10 in the refined condition number, improving the rate to (1 − κ′−3/2 ) with κ′ ∶= β ′ /α′ ≤ κ and the iteration count of Eq. (59) to κ′3/2 log(1/ϵ). The proof of convergence with projection onto the refined interval does not carry over automatically for the following reasons. The projection step of Proposition 9 relied on each Ri being contained in [αI, βI], which fails on the refined interval [α′ I, β ′ I]. What is needed instead is that clipping to the refined interval does not increase f1 , which we prove as Proposition 15. We nonetheless state our results on [αI, βI], whose constants are read off the individual inputs and match the conventions of the literature (the two coincide in the worst case); since our κ′ is no larger than the refined condition number of [ACGS21], the improvement in Table 1 is meaningful instance-wise: κ′3/2 against O(κ′5/2 ). The analysis of the refined bounds is deferred to Appendix C.

4.1

Choosing the iteration count −3/2

From f (St ) − f (S⋆ ) ≤ (1 − κ−3/2 )t (f (S0 ) − f (S⋆ )) ≤ e−tκ

(f (S0 ) − f (S⋆ )), it suffices to take

T ≥ κ3/2 log(1/ϵ)

(59)

to reach relative ϵ accuracy in function value, f (ST ) − f (S⋆ ) ≤ ϵ (f (S0 ) − f (S⋆ )), matching the convention of Table 1. This a priori count is conservative in practice; the following a posteriori alternative is sharper. 18

Corollary 12 (A posteriori certificate). Let f ∈ {f1 , f2 } and let (St )t≥0 be generated by Algorithm 1 from S0 ∈ [αI, βI] (additionally S0 ∈ P when f = f2 ). For every t ≥ 0, f (St ) − f (S⋆ ) ≤

κ3/2 B (St , K(St )) . 2

(60)

Consequently, for an absolute tolerance δ > 0 one may terminate once B(St , K(St )) ≤ 2κ−3/2 δ, at which point f (St ) − f (S⋆ ) ≤ δ; taking δ = ϵ (f (S0 ) − f (S⋆ )) recovers the relative accuracy of Eq. (59). 2 ˆ Proof. Substitute the step-length identity B(St , K(St )) = ∥∇g(A t )∥F of Proposition 9 into the PL inequality of Lemma 8 and use 4µα = κ−3/2 .

The certificate quantity is the (squared) length of the step the iteration has just computed; in practice the resulting count is far smaller than Eq. (59).8

5

Numerical Experiments (a) varying κ (d = 15, n = 10)

[f(St) − f(S ⋆ )] / [f(S0) − f(S ⋆ )]

10

(b) varying d (n = 10, κ = 1000)

0

10−2

10−4 unit step

10−6

small step

κ 03/2 bound

κ = 101 κ = 10

κ = 103

10

unit step

small step

κ 03/2 bound

103

104

105

d=4 d = 16 d = 64

2

−8

10−10

10−12

100

101

102

103

104

105

iteration t

100

101

102

iteration t

Figure 2: Unit-step (blue) against small-step (orange) RGD on pinned ensembles (every Ri of condition number exactly κ), 12 instances per band (translucent) with opaque medians; darker shades for larger parameter values, gaps normalized by their initial value. (a) Varying κ at d = 15: median iterations to the numerical floor range over 6 to 16 for the unit step, against 175 to 16288 for the small step. (b) Varying d at κ = 103 : the unit-step counts are essentially flat over a sixteenfold increase in d, as Theorem 10 asserts; the small-step counts track the conditioning κ′ , not the dimension (see text). Dotted curves show the a priori guarantee of Corollary 17 at each setting’s median κ′ ; the clip never activates, so the projected and unprojected unit steps coincide. Note the logarithmic iteration axis.

We now discuss the numerical experiments accompanying our theoretical results. The scripts and notebook that generate every figure and number reported here are available at the accompanying repository [Afh26]. First, we compare our Projected RGD algorithm against the standard (unprojected) RGD algorithm and the small-step RGD algorithm of [ACGS21], run at the step size η = α′ /(2β ′ ), no smaller than its prescription α′ /(2 ∑ni=1 wi λmax (Ri )) and Across the experiments of Section 5, the certificate (Corollary 12) overestimates the true gap (f (St ) − f (S⋆ )) by roughly its κ3/2 /2 prefactor. However, since the convergence is linear, this multiplicative slack costs only a few additional iterations, and the certified stopping time improves on the a priori count of Eq. (59) by orders of magnitude. 8

19

hence conservative in its favor. On the instance of Example 13 (shown in Fig. 1(b)), we ran projected against unprojected RGD: the projection is active on exactly the two exit steps, and the two trajectories are indistinguishable in convergence rate thereafter. On the same instance the small step, η = α′ /(2β ′ ) = 1.33 × 10−2 , remains feasible throughout but needs 544 iterations to reach a 10−6 relative gap and 1108 to reach 10−12 , over a window in which the unit step converges to machine precision. Unit-step RGD reaches the same relative gaps of 10−6 and 10−12 at iterations 3 and 6 respectively (projected and unprojected alike), a factor of about 180 fewer iterations. Numerics also suggest that unless the ensemble and the initial iterate are adversarially designed, all iterates remain in [αI, βI]: in every subsequent (non-adversarial) experiment the clipping never activates, so the projected and unprojected unit-step RGDs coincide, and the projection acts as a dormant safeguard that converts the numerically observed dimension-independent linear convergence into a guarantee. The primary comparison, in terms of performance, is then between the (projected) unitstep and the small-step RGD. We demonstrate this in Fig. 2(a), which shows RGD exhibits a far superior performance at unit step than at small step: at κ = 103 the unit step reaches the numerical floor in 16 iterations against 16288 for the small step, three orders of magnitude apart. This empirical superiority has been observed previously [ACGS21; BRT25], but a theoretical backing was lacking. Our key contribution is to ground this empirical performance in theory, at the cost of a projection that is (i) essentially free and (ii) almost never needed in practice. We then empirically study the rate of convergence for RGD at unit- and smallstep sizes. As seen in Fig. 2(a), the observed contraction per iteration is far faster than the guaranteed rate (1−κ−3/2 )—and remains so under the refined rate (1−κ′−3/2 ) of Remark 11— so the iteration counts fall well short of the a priori count κ3/2 log(1/ϵ) of Eq. (59), suggesting room for improvement in the constants of Theorem 10. Fixing κ and grading instead by the dimension, Fig. 2(b) isolates the dimensionindependence of the guarantee. For the sweep to be meaningful the effective difficulty must not itself vary with d: with i.i.d. spectra, mini λmin (Ri ) is a minimum of d draws and collapses onto α as d grows, so the refined constants—and with them the small-step size and the sharpened bound—would drift for reasons unrelated to the rate; the pinned ensembles of Fig. 2 remove this drift, leaving only a mild, saturating decrease of β ′ as the eigenbases decorrelate. The unit-step bands for d ∈ {4, 16, 64} then nearly coincide, as do the small-step bands and the sharpened bounds, the residual spread of all three tracking κ′ rather than d: the iteration count is governed by the conditioning, not the dimension, as Theorem 10 asserts.

6

Conclusion

We introduced and analyzed Projected BW-GD: unit-step RGD composed, at every iteration, with the BW projection onto the well-conditioned set, which we showed is eigenvalue clipping (Lemma 6). For both the BW barycenter and the invariant matrix projection problem, the algorithm converges linearly at the dimension-independent rate (1 − κ−3/2 )—the first such guarantee at the unit step size used in practice, and a polynomial improvement, from O(κ′5/2 ) to κ′3/2 , in the iteration complexity of the best dimension-free guarantee—with an a posteriori certificate permitting termination on the fly. We do not establish the rate for the unprojected fixed-point algorithm; rather, an eigenvalue clip of negligible cost converts it into an algorithm with a worst-case guarantee, while coinciding with it on typical inputs. The analysis runs on the flat square-root manifold, where Euclidean strong convexity becomes a PL inequality and the unit-step descent lemma is a variance identity. This perspective greatly simplifies the technical analysis by avoiding notions such as convexity along generalized geodesics. Both ingredients require well-conditioned iterates—a requirement that Example 13 shows can 20

genuinely fail for the unprojected iteration—and the projection enforces it without undoing the progress of the gradient step. The invariant matrix projection problem inherits this picture in full, since the invariant submanifold is totally geodesic and closed under the clip (Appendix E); the rate thereby extends to the quantum-information problems that the formulation unifies [BRT25, Sec. 7]. Two directions are suggested by the experiments. First, from S0 = α+β 2 I the unprojected trajectory remained in [αI, βI] in all our numerics, and the exits of Example 13 required an adversarial initial point; proving that unprojected unit-step RGD preserves [αI, βI] from this canonical initialization would render the projection unnecessary for a default S0 . Second, the observed rates are far faster than the bound, even after the sharpening of Appendix C: the improvement from κ5/2 to κ3/2 is an improvement of the guarantee, and both guarantees are loose against what the iteration actually does. Whether the generic-case rate admits a sharper exponent, or the worst case is attained only on rare instances, remains open.

Acknowledgments The author thanks Marco Tomamichel and Roberto Rubboli for helpful discussions. This project is supported by the NRF Investigatorship award (NRF-NRFI10-2024-0006) and the CQT Young Researcher Career Development Grant.

Use of Large Language Models The author acknowledges the use of Large Language Models (LLMs) in the preparation of this article. This includes standard usage such as proofreading, restructuring, basic questions, and coding. In terms of meaningful contribution, the latter half of Lemma 6 was proven using Google Gemini 3.1 Pro. The key prompt, paraphrased, was “Show that (two-sided) clipping of eigenvalues to an interval is a contraction with respect to the BW distance”, which was inspired by the existing result [ACGS21, Prop. 3], which showed that clipping eigenvalues from above is a BW contraction. As mentioned in the main text of the article, this is important as existing proof methods were unable to guarantee dimension-independent convergence (at unit step-size) because the iterates could dip below αI, which can now be remedied by clipping since it acts as a contraction. The LLM-generated proof was then reworked into the current form for improved presentation. Given the proof for contraction, it was then straightforward to show that eigenvalue-clipping constituted a BW-projection onto [αI, βI] (former half of Lemma 6). Other technical results obtained with the help of LLMs (primarily Claude Fable) include Proposition 15, which proved that clipping onto the refined interval [α′ I, β ′ I] does not increase the objective value, Corollary 18, which showed that singular-value clipping corresponds to Frobenius projection onto the set of square-roots of the well-conditioned interval, and the counterexample instance of Example 13.9

References [Afh25]

A Afham. “Optimizing and Generalizing Quantum Fidelities”. Thesis. 2025. url: https://opus.lib.uts.edu.au/handle/10453/190570.

[Afh26]

A. Afham. Companion code for “Projected Riemannian Gradient Descent for the Bures–Wasserstein Barycenter”. 2026. url: https://github.com/afhamash/ projected-rgd.

9

See Appendix A for the search strategy.

21

[AF25]

A. Afham and Chris Ferrie. “Riemannian-Geometric Generalizations of Quantum Fidelities and Bures-Wasserstein Distance”. In: Journal of Mathematical Physics 66.8 (2025).

[AKF22]

A. Afham, Richard Kueng, and Chris Ferrie. Quantum Mean States Are Nicer than You Think: Fast Algorithms to Compute States Maximizing Average Fidelity. June 2022. arXiv: 2206.08183 [quant-ph]. url: https://arxiv.org/ abs/2206.08183.

[AT26]

A. Afham and Marco Tomamichel. Projections with Respect to Bures Distance and Fidelity: Closed-Forms and Applications. 2026. arXiv: 2602.14732 [quant-ph]. url: https://arxiv.org/abs/2602.14732.

[AC11]

Martial Agueh and Guillaume Carlier. “Barycenters in the Wasserstein Space”. In: SIAM Journal on Mathematical Analysis 43.2 (2011), pp. 904–924.

[AEO26]

Vitalii Aksenov, Martin Eigel, and Mathias Oster. “Anderson Mixing in Bures Wasserstein Space of Gaussian Measures”. In: arXiv preprint arXiv:2601.22038 (2026).

[Alb83]

Peter M Alberti. “A note on the transition probability over C*-algebras”. In: Letters in Mathematical Physics 7.1 (1983), pp. 25–32.

[ACGS21]

Jason Altschuler, Sinho Chewi, Patrik R Gerber, and Austin Stromme. “Averaging on the Bures-Wasserstein Manifold: Dimension-Free Convergence of Gradient Descent”. In: Advances in Neural Information Processing Systems 34 (2021), pp. 22132–22145.

[ÁDCM16]

Pedro C Álvarez-Esteban, E Del Barrio, JA Cuesta-Albertos, and C Matrán. “A Fixed-Point Approach to Barycenters in Wasserstein Space”. In: Journal of Mathematical Analysis and Applications 441.2 (2016), pp. 744–762.

[AGS08]

Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré. Gradient Flows: In Metric Spaces and in the Space of Probability Measures. Springer Science & Business Media, 2008.

[Bei13]

Salman Beigi. “Sandwiched Rényi divergence satisfies data processing inequality”. In: Journal of Mathematical Physics 54.12 (2013).

[Bha09]

Rajendra Bhatia. “Positive Definite Matrices”. In: Positive Definite Matrices. Princeton university press, 2009. doi: 10.1515/9781400827787.

[BJL19]

Rajendra Bhatia, Tanvi Jain, and Yongdo Lim. “On the Bures–Wasserstein Distance between Positive Definite Matrices”. In: Expositiones Mathematicae 37.2 (2019), pp. 165–191.

[BJL18]

Rajendra Bhatia, Tanvi Jain, and Yongdo Lim. “Strong Convexity of Sandwiched Entropies and Related Optimization Problems”. In: Reviews in Mathematical Physics 30.09 (2018), p. 1850014.

[BRT25]

Shrigyan Brahmachari, Roberto Rubboli, and Marco Tomamichel. “A fixedpoint algorithm for matrix projections with applications in quantum information: S. Brahmachari et al.” In: Mathematical Programming (2025), pp. 1–31.

[BM03]

Samuel Burer and Renato DC Monteiro. “A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization”. In: Mathematical programming 95.2 (2003), pp. 329–357.

[BM05]

Samuel Burer and Renato DC Monteiro. “Local minima and convergence in lowrank semidefinite programming”. In: Mathematical programming 103.3 (2005), pp. 427–444. 22

[CMRS20]

Sinho Chewi, Tyler Maunu, Philippe Rigollet, and Austin J Stromme. “Gradient Descent Algorithms for Bures-Wasserstein Barycenters”. In: Conference on Learning Theory. PMLR, 2020, pp. 1276–1304.

[CNR25]

Sinho Chewi, Jonathan Niles-Weed, and Philippe Rigollet. Statistical optimal transport. Springer, 2025.

[Fan+24]

Junyi Fan et al. “On the Convergence of Projected Bures-Wasserstein Gradient Descent under Euclidean Strong Convexity”. In: Proceedings of the 41st International Conference on Machine Learning. Vol. 235. Proceedings of Machine Learning Research. PMLR, 21–27 Jul 2024, pp. 12832–12857. url: https: //proceedings.mlr.press/v235/fan24b.html.

[GW15]

Manish K Gupta and Mark M Wilde. “Multiplicativity of completely bounded p-norms implies a strong converse for entanglement-assisted capacity”. In: Communications in Mathematical Physics 334.2 (2015), pp. 867–887.

[Joz94]

Richard Jozsa. “Fidelity for Mixed Quantum States”. In: Journal of modern optics 41.12 (1994), pp. 2315–2323.

[KNS16]

Hamed Karimi, Julie Nutini, and Mark Schmidt. “Linear convergence of gradient and proximal-gradient methods under the polyak-łojasiewicz condition”. In: Joint European conference on machine learning and knowledge discovery in databases. Springer. 2016, pp. 795–811.

[KRS09]

Robert Konig, Renato Renner, and Christian Schaffner. “The Operational Meaning of Min- and Max-Entropy”. In: IEEE Transactions on Information Theory 55.9 (2009), pp. 4337–4347. doi: 10.1109/tit.2009.2025545. url: https://doi.org/10.1109%5C%2Ftit.2009.2025545.

[KSS21]

Alexey Kroshnin, Vladimir Spokoiny, and Alexandra Suvorikova. “Statistical Inference for Bures–Wasserstein Barycenters”. In: The Annals of Applied Probability 31.3 (2021), pp. 1264–1298.

[Lam+22]

Marc Lambert et al. “Variational inference via Wasserstein gradient flows”. In: Advances in Neural Information Processing Systems 35 (2022), pp. 14434– 14447.

[LL17]

Thibaut Le Gouic and Jean-Michel Loubes. “Existence and consistency of Wasserstein barycenters”. In: Probability Theory and Related Fields 168.3 (2017), pp. 901–917.

[LPRS22]

Thibaut Le Gouic, Quentin Paris, Philippe Rigollet, and Austin J Stromme. “Fast convergence of empirical barycenters in Alexandrov spaces and the Wasserstein space”. In: Journal of the European Mathematical Society 25.6 (2022), pp. 2229–2250.

[LM08]

Adrian S Lewis and Jérôme Malick. “Alternating projections on manifolds”. In: Mathematics of Operations Research 33.1 (2008), pp. 216–234.

[Liu+17]

CL Liu et al. “A new coherence measure based on fidelity”. In: Quantum Information Processing 16.8 (2017), p. 198.

[MMP18]

Luigi Malago, Luigi Montrucchio, and Giovanni Pistone. “Wasserstein Riemannian Geometry of Positive Definite Matrices”. In: arXiv preprint arXiv:1801.09269 (2018). arXiv: 1801.09269.

[MLR23]

Tyler Maunu, Thibaut Le Gouic, and Philippe Rigollet. “Bures-wasserstein barycenters and low-rank matrix recovery”. In: International Conference on Artificial Intelligence and Statistics. PMLR. 2023, pp. 8183–8210. 23

[ONe66]

Barrett O’Neill. “The fundamental equations of a submersion.” In: Michigan Mathematical Journal 13.4 (1966), pp. 459–469.

[PZ19]

Victor M Panaretos and Yoav Zemel. “Statistical Aspects of Wasserstein Distances”. In: Annual review of statistics and its application 6.1 (2019), pp. 405– 431.

[Pol63]

Boris Teodorovich Polyak. “Gradient methods for minimizing functionals”. In: Zhurnal vychislitel’noi matematiki i matematicheskoi fiziki 3.4 (1963), pp. 643– 653.

[Tom13]

Marco Tomamichel. A Framework for Non-Asymptotic Quantum Information Theory. 2013. arXiv: 1203.2142 [quant-ph]. url: https://arxiv.org/abs/ 1203.2142.

[Tom15]

Marco Tomamichel. Quantum Information Processing with Finite Resources: Mathematical Foundations. Vol. 5. Springer, 2015.

[Uhl76]

Armin Uhlmann. “The “Transition Probability” in the State Space of a StarAlgebra”. In: Reports on Mathematical Physics 9.2 (1976), pp. 273–279.

[Uhl11]

Armin Uhlmann. “Transition Probability (Fidelity) and Its Relatives”. In: Foundations of physics 41.3 (2011), pp. 288–298.

[Wat18]

John Watrous. The Theory of Quantum Information. Cambridge: Cambridge University Press, 2018. isbn: 978-1-107-18056-7. doi: 10.1017/9781316848142.

[WG03]

Tzu-Chieh Wei and Paul M Goldbart. “Geometric measure of entanglement and applications to bipartite and multipartite quantum states”. In: Physical Review A 68.4 (2003), p. 042307.

[WY16]

Andreas Winter and Dong Yang. “Operational resource theory of coherence”. In: Physical review letters 116.12 (2016), p. 120404.

[ZP19]

Yoav Zemel and Victor M. Panaretos. “Fréchet Means and Procrustes Analysis in Wasserstein Space”. In: Bernoulli. Official Journal of the Bernoulli Society for Mathematical Statistics and Probability 25.2 (2019), pp. 932–976. doi: 10. 3150/17-BEJ1009. url: https://doi.org/10.3150/17-BEJ1009.

[ZHVZ25]

Shixin Zheng, Wen Huang, Bart Vandereycken, and Xiangxiong Zhang. “Riemannian optimization using three different metrics for Hermitian PSD fixedrank constraints: S. Zheng et al.” In: Computational Optimization and Applications 91.3 (2025), pp. 1135–1184.

A

Unit-Step RGD Exits the Well-Conditioned Set

We now discuss a barycenter instance, where all the matrices involved—the elements of the ensemble, the initial point, and the barycenter—lie strictly in the interior of a compact interval [αI, βI], yet the trajectory exits the interval (transiently). Fig. 3 displays the trajectory and the formal statement follows. This is the same instance, and the same trajectory, as panel (b) of the headline Fig. 1; here it is shown in the complementary view—the full twoeigenvalue picture in the main panel, the signed deviation in the inset—alongside the formal statement. cos(θ) sin(θ) ) be − sin(θ) cos(θ) the rotation matrix by angle θ, and consider the barycenter instance with weights w = Example 13 (Transient exit of Unit-step RGD). Let d = 2, Qθ ∶= (

24

β λmin (St) − α, symlog

100 10−2

102

eigenvalues of St

−4 10−4 0 −10

-4.43e-03

−10−2 −100 0

10

feasible band [αI, βI]

-2.98e-01 1

2

3

4

5

λmax (St)

6

λmin (St)

1

λmin < α

α 100 0

1

2

3

4

5

6

7

8

iteration t

Figure 3: The trajectory of Example 13. Main panel: both eigenvalues of the unprojected unit-step RGD iterates (log scale), with the feasible band [αI, βI] shaded and the infeasible iterates marked. Inset: the signed deviation λmin (St ) − α on a symmetric-log scale, resolving the two undershoots (−2.98 × 10−1 and −4.43 × 10−3 ) against the limit λmin (S⋆ ) − α ≈ +2.52 × 10−3 : the iterate exits [αI, βI] at t = 1, is still outside at t = 2, and re-enters at t = 3.

(0.12, 0.827, 0.053) and ensemble elements 1.05 0 1.0501 0 1.0502 0 R1 = Qθ1 ( ) Q⊺ , R2 = Qθ2 ( ) Q⊺ , R3 = Qθ3 ( ) Q⊺ , 0 315 θ1 0 1.575 θ2 0 4.41 θ3 (61) where the angles are (θ1 , θ2 , θ3 ) = (0.246, 0.338, 0.253). We choose α = 1.05 and β = 315 (κ = 300), together with the strictly interior initialization S0 = Qθ0 (

1.11 0 ) Q⊺ ∈ interior [αI, βI]; 0 310 θ0

θ0 = 1.473.

(62)

Although the barycenter admits no closed form, its location is certified by the eigenvalue bounds of [ACGS21, Thm. 6]: n

√

2

′

λmin (S⋆ ) ≥ ( ∑ wi λmin (Ri )) =∶ α ,

n

√

2

λmax (S⋆ ) ≤ (∑ wi λmax (Ri )) = 10.75,

i=1

(63)

i=1

where α′ = α + 9.33 × 10−5 , so S⋆ lies strictly in the interior of [αI, βI]. Running the unprojected unit-step RGD iteration from S0 yields, in double precision, ⎧ −2.984 × 10−1 ⎪ ⎪ ⎪ ⎪ λmin (St ) − α = ⎨ −4.434 × 10−3 ⎪ ⎪ ⎪ −3 ⎪ ⎩ +2.117 × 10

t = 1, t = 2, t = 3.

(64)

The trajectory therefore exits [αI, βI] at t = 1 and remains outside for two consecutive iterations—indeed, both excursion iterates lie below even the certified floor α′ of the optimum, so the excursion cannot be attributed to the iterates approaching a near-boundary optimum— before re-entering at t = 3 and converging (numerically, λmin (St ) → α + 2.52 × 10−3 ). We choose n = 3 deliberately as for n = 2 the barycenter is a point on the geodesic joining the targets and admits a closed form. The example extends to every problem size: to any n ≥ 3 by splitting R3 into identical copies whose weights sum to w3 , which leaves 25

the iteration unchanged, and to any d ≥ 2 by padding all matrices with a block c Id−2 , c ∈ [α, β], under which the dynamics act block-diagonally and the violations are unchanged. The instance was found using an LLM-driven search (Claude Fable): a randomized sweep over small configurations (d = 2, n = 3) first located rare single-step exits, and an adversarial refinement of the angles, spectra, and weights then deepened these into the two-consecutivestep excursion above, while keeping the ensemble, the initial point, and the certified location of the optimum strictly interior.

B

A Variance Inequality Proof of the Descent Lemma

The descent lemma used in Section 4 is [BJL19, Theorem 10], proved there with tools from optimal transport. We give a self-contained proof of the same result, phrased in the squareroot manifold, where it reduces to the variance identity for a mean in a Euclidean space. Proposition 14. Let (Ri , wi )i∈[n] be a weighted ensemble and g1 ∶= f1 ○ π. For any S ∈ PD(d), 1 1 1 g1 (S 2 ) ≥ g1 (A′ ) + ∥∇g1 (S 2 )∥2F , (65) 2 1

1

where A′ = S 2 − ∇g1 (S 2 ); consequently 1 f1 (S) ≥ f1 (K(S)) + B(S, K(S)). 2

(66)

Proof. For any (Yi )i∈[n] ⊆ L(d) with mean Y ∶= ∑ni=1 wi Yi and any X ∈ L(d), expanding both sides verifies the variance identity n

n

2

2

2

∑ wi ∥X − Yi ∥F = ∑ wi ∥Y − Yi ∥F + ∥X − Y ∥F .

i=1

(67)

i=1

1

1

1

1

1

1

Apply this with X ≡ S 2 and Yi ≡ Ri2 Vi† , where Vi = Pol(S 2 Ri2 ), so that each (S 2 , Ri2 Vi† ) 1

1

is an optimal pairing: ∥S 2 − Ri2 Vi† ∥2F = B(Ri , S). The left-hand side of Eq. (67) is then 1

1

1

1

1

2f1 (S) = 2g1 (S 2 ). The chain rule ∇g1 (S 2 ) = 2∇f1 (S)S 2 = ∑ni=1 wi [S 2 − Ri2 Vi† ] (the second equality is Proposition 21 applied termwise) can be rearranged as n

1

Y = ∑ wi Ri2 Vi† = S 2 − ∇g1 (S 2 ) = A′ , 1

1

(68)

i=1 1

so ∥X − Y ∥2F = ∥∇g1 (S 2 )∥2F . Finally, Y = A′ ∈ π −1 [K(S)] and Yi ∈ π −1 [Ri ], so the variational characterization of the BW distance gives ∥Y −Yi ∥2F ≥ B(Ri , K(S)), whence ∑ni=1 wi ∥Y −Yi ∥2F ≥ 2f1 (K(S)) = 2g1 (A′ ). Substituting the three expressions into Eq. (67) and dividing by two 1 1 gives the first claim; the second follows since ∥∇g1 (S 2 )∥2F = ∥S 2 − A′ ∥2F ≥ B(S, K(S)). We now see where the non-negative curvature of the manifold yields the descent lemma. The total manifold—equipped with the Frobenius metric—is flat, and there the statement is the Euclidean variance identity for the mean of the optimally paired square-roots. The inequality arises from descending to the base manifold, where the natural (BW) distance is obtained by minimizing the Frobenius distance over the fibres (Eq. (14)). Consequently, the quotient map contracts distances. By O’Neill’s formula for Riemannian submersions, sectional curvature can only increase under the quotient [ONe66], so the BW manifold inherits non-negative curvature from the flat total space; on Wasserstein space at large, the 26

non-negative curvature holds in the Alexandrov sense [AGS08, Thm. 7.3.2]. In the proof, the contraction enters with the favorable sign: the averaged square roots project (down to the base manifold via π) to precisely the next RGD iterate, which is at least as close to each reference as the flat identity accounts for, so the gradient-norm term survives as a genuine per-step decrease. This is an instance of a general phenomenon: non-negative curvature yields geodesic smoothness of the barycenter functional, and with it unit-step descent, on the whole Wasserstein space [CMRS20, Thm. 7].

C

Descent under Projection onto the Refined Interval

√ 2 Throughout this appendix, f1 is the barycenter functional and α′ = [∑ni=1 wi λmin (Ri )] , β ′ ∶= λmax (∑ni=1 wi Ri ) are the refined spectral bounds of Eq. (7); they satisfy α ≤ α′ ≤ β ′ ≤ β, enclose the optimum (α′ I ≤ S⋆ ≤ β ′ I; floor by [ACGS21, Thm. 6], ceiling by [BJL19, Thm. 9]), and satisfy β ′ ≤ ∑ni=1 wi λmax (Ri ), so κ′ is no larger than the refined condition number in which the small-step guarantee of [ACGS21, Thm. 3] is stated. Thus it is natural to ask if clipping onto the tighter interval [α′ I, β ′ I] will still lead to convergence of the Projected BW-GD algorithm, which does not follow from the arguments used in Section 4. This is because the descent argument of Proposition 9 relied on the fact that each ensemble element Ri lies in [αI, βI] and is hence a fixed point of clipα,β , which is what forced f1 (clipα,β [K(S)]) ≤ f1 (K(S)). The refined interval does not contain the ensemble— individual Ri may have eigenvalues on either side of [α′ , β ′ ]—so the fixed-point argument is unavailable. By working with the objective function directly, we can show that clipping to this refined interval decreases the objective value, which allows us to use it in Algorithm 1. Proposition 15 (Refined clipping does not increase the objective). For every S ∈ PD(d), it holds that f1 (clipα′ ,β ′ (S)) ≤ f1 (S). (69) Proof. Let Ŝ ∶= clipα′ ,β ′ (S), and let S = ∑dj=1 λj vj vj† be the eigendecomposition; then Ŝ = †

∑j clipα′ ,β ′ (λj )vj vj shares the eigenbasis, and d

∆ ≡ S − Ŝ = ∑ δj vj vj† ,

δj ∶= λj − clipα′ ,β ′ (λj ).

(70)

j=1

The convexity of f1 on PD(d) implies the following first-order convexity condition: f1 (S) ≥ f1 (Ŝ) + ⟨∇f1 (Ŝ), S − Ŝ⟩.

(71)

Thus it suffices to show ⟨∇f1 (Ŝ), S − Ŝ⟩ ≥ 0. By Eq. (15), ∇f1 (Ŝ) = 21 (I − T ) with T ∶= ∑ni=1 wi Ti and Ti ∶= Ŝ −1 #Ri , where each Ti is the unique positive definite solution of the Riccati equation Ti ŜTi = Ri [Bha09]. Expanding in the shared eigenbasis, ⟨∇f1 (Ŝ), S − Ŝ⟩ =

1 d ∑ δj (1 − ⟨vj , T vj ⟩). 2 j=1

(72)

We will show that each summand is non-negative. Since λj ∈ [α′ , β ′ ] ⇒ δj = 0, we can restrict our attention to the ceiling case: λj ≥ β ′ , and the floor case: λj ≤ α′ . Let us first consider the ceiling case. For j such that δj = λj − β ′ > 0, we have Ŝvj = β ′ vj . For each i, inserting

27

Ŝ 2 Ŝ − 2 and applying Cauchy–Schwarz, 1

1

1 2

⟨vj , Ti vj ⟩ = ⟨Ŝ Ti vj , Ŝ

− 21

1 2

vj ⟩ ≤ ∥Ŝ Ti vj ∥ ∥Ŝ

− 21

vj ∥ =

√

⟨vj , Ti ŜTi vj ⟩ ⋅ ¿ Á ⟨vj , Ri vj ⟩ À =Á , β′

√

⟨vj , Ŝ −1 vj ⟩ (73)

where we used Ti ŜTi = Ri and Ŝ −1 vj = vj /β ′ . Averaging and using concavity of the square root, ¿ n Á ⟨vj , (∑ni=1 wi Ri ) vj ⟩ 1 n √ À ⟨vj , T vj ⟩ = ∑ wi ⟨vj , Ti vj ⟩ ≤ √ ′ ∑ wi ⟨vj , Ri vj ⟩ ≤ Á β′ β i=1 i=1 (74) √ λmax (∑ni=1 wi Ri ) = 1, ≤ β′ the final equality by the definition of β ′ ; hence δj (1 − vj† T vj ) ≥ 0 as required. We now consider the floor case. For j such that δj = λj − α′ < 0, we have Ŝvj = α′ vj . Here the required result follows from operator monotonicity of the geometric mean alone. Since Ri ≥ λmin (Ri )I, operator monotonicity of the geometric mean gives the operator inequality √ 1 (75) Ti = Ŝ −1 #Ri ≥ Ŝ −1 # (λmin (Ri )I) = λmin (Ri )Ŝ − 2 . √ 1 Evaluating the quadratic form at vj , where Ŝ − 2 vj = vj / α′ , and averaging, √ n 1 n √ α′ ⟨vj , T vj ⟩ = ∑ wi ⟨vj , Ti vj ⟩ ≥ √ ∑ wi λmin (Ri ) = √ = 1, (76) α′ i=1 α′ i=1 the final equality by the definition of α′ ; hence δj (1 − ⟨vj , T vj ⟩) ≥ 0 in this case as well. Every term of Eq. (72) is non-negative, so ⟨∇f1 (Ŝ), S − Ŝ⟩ ≥ 0, completing the proof. The second ingredient of the analysis of Section 4 is the strong-convexity modulus feeding the PL inequality, and on the refined interval it cannot be borrowed as stated: the modulus of Eq. (17) is derived for references and iterates sharing the interval [αI, βI], whereas here the iterates live in [α′ I, β ′ I] while individual Ri may lie outside it. However, no new work is needed because the strong convexity of Eq. (16) is already mixed: its constant couples λmin (R) of the reference with an upper bound b on the argument alone, with no common interval required. The following lemma records the specialization to f1 and b = β ′ , producing the modulus µ′ whose combination with the PL inequality yields the κ′ -rate. Lemma 16 (specialization of Eq. (16); see [BJL18, proof of Thm. 1]). The barycenter functional f1 is µ′ -SC on {S ∈ PD(d) ∶ S ≤ β ′ I}, where √

3/2

β′ ≤ κ. α′

(77)

√ n √ 1 n 1 α′ 2 ∇ f1 (S) = ∑ wi ∇ BRi (S) ≥ w Id ≡ µ′ ⋅ Id, λ (R ) ⋅ Id = ∑ i min i 3/2 3/2 ′ ′ 2 i=1 4β 4β i=1

(78)

′

µ ∶=

4β

α′

, ′3/2

so that

4µ′ α′ = (

α′ ) β′

= κ′−3/2 ,

κ′ ∶=

Proof. Let S ≤ β ′ I. Applying Eq. (16) to each BRi with b = β ′ , 2

where Id ∶ L(d) → L(d) is the identity map on the matrix space and the last equality by the definition of α′ . 28

With the two ingredients in place, the analysis of Section 4 runs on the refined interval verbatim. Corollary 17 (Convergence on the refined interval). Let S0 ∈ [α′ I, β ′ I] and let the iterates be generated by unit-step BW-GD projected onto the refined interval, St+1 ∶= clipα′ ,β ′ [K(St )]. Then f1 (St ) − f1 (S⋆ ) ≤ (1 −

1

κ

t

) (f1 (S0 ) − f1 (S⋆ )) , ′3/2

(79)

the verbatim analogue of Theorem 10 in the refined constants, and the iteration count of Eq. (59) improves accordingly to κ′3/2 log(1/ϵ). Proof. Write εt ∶= f1 (St ) − f1 (S⋆ ) and fix At ∈ π −1 [St ]; all iterates (by construction of the clipping step) and S⋆ lie in [α′ I, β ′ I]. Chaining the unit-step descent inequality (Proposition 14) with Proposition 15, and using the step-length identity B(St , K(St )) = ∥∇g1 (At )∥2F (Lemma 5), gives εt+1 ≤ εt − 21 ∥∇g1 (At )∥2F . By Lemma 16, f1 is µ′ -SC on {S ∈ PD(d) ∶ S ≤ β ′ I} ⊇ [α′ I, β ′ I], so Lemma 8 (applied with [aI, bI] = [α′ I, β ′ I] and µ0 = µ′ ) gives κ′−3/2 εt = 4µ′ α′ εt ≤ 12 ∥∇g1 (At )∥2F . Combining the two and unrolling yields the rate; the iteration count follows as in Section 4.1. The canonical initialization α+β 2 I of Algorithm 1 may lie outside the refined interval; one ′ ′ ′ ′ clips it onto [α I, β I] first, or initializes at α +β 2 I. We note that the statement is specific to the barycenter problem: for the invariant projection problem (a single reference R), the refined constants reduce to α′ = λmin (R) = α and β ′ = λmax (R) = β, and nothing is gained. Moreover, in the worst case over ensembles, α′ = α and β ′ = β are attained, so the refinement improves the constants on favorable instances while leaving the worst-case rate of Theorem 10 unchanged; on every instance, the count κ′3/2 log(1/ϵ) is expressed in a condition number no larger than that of the O(κ′5/2 log(1/ϵ)) small-step guarantee of [ACGS21], which it therefore dominates.

D

The Projected Iteration on the Total Manifold

The base–total equivalence yields a closed-form implementation of both algorithms of Section 3.3, which we now describe, beginning by transporting the entire iteration—projection included—to the total manifold. Corollary 18 (Singular-value clipping is Frobenius projection). For every A ∈ GL(d), √ √ is the unique Frobenius projection onto the square-roots singular-value clipping clipsv α, β of the well-conditioned set: √ √ (A) = clipsv α, β

argmin

B∈π −1 [[αI,βI]]

∥A − B∥F .

(80)

Moreover, the two projections are intertwined by the quotient map: √ √ = clip π ○ clipsv α,β ○ π. α, β

(81)

Proof. Results of this type are well known for spectral sets: the Frobenius projection onto a set defined by symmetric constraints on the eigenvalues or singular values reduces to projecting the spectrum itself onto the corresponding set of scalars [LM08]. Rather than 29

specializing the general theory, we give a short proof, deriving uniqueness from the Projection Lemma through the quotient structure. Recall that λi (BB † ) = σi (B)2 for any B ∈ L(d). Fix a singular value decomposition A = U ΣV † with σ1 (A) ≥ ⋯ ≥ σd (A) > 0, and write √ √ (A) = U clip√ √ (Σ) V † , which is feasible. Â ≡ clipsv α, β α, β We appeal to the fibre geometry of the total manifold. Let S ∶= π(A) and Ŝ ∶= π(Â) = U clip√α,√β (Σ)2 U † = clipα,β (S)—which is precisely the claimed intertwining identity. More1

1

over, (A, Â) is an optimal pairing for (S, Ŝ): the two commute, so F(S, Ŝ) = Tr[S 2 Ŝ 2 ] = ∑di=1 σi (A) clip√α,√β (σi (A)) = R⟨A, Â⟩, whence ∥A − Â∥2F = B(S, Ŝ). Now let B be any feasible point, with Q ∶= π(B) ∈ [αI, βI]. Then ∥A − B∥2F ≥

min

B ′ ∈π −1 [Q]

∥A − B ′ ∥2F = B(S, Q) ≥ B(S, Ŝ) = ∥A − Â∥2F .

(82)

The first inequality holds simply because B lies in the fibre of Q. The equality that follows is the quotient formula Eq. (14): the Frobenius distance is invariant under right multiplication by unitaries, so the minimum over both fibres is attained already from the fixed representative A. The second inequality is the Projection Lemma (Lemma 6), as Ŝ = clipα,β (S); the final equality is the optimal pairing established above. In particular, since  is feasible, the chain shows that it is a minimizer. If B is a minimizer, then ∥A − B∥2F = ∥A − Â∥2F , so every step in Eq. (82) holds with equality. Equality in the Projection Lemma step forces Q = Ŝ, since the BW projection is unique (Lemma 6); thus B lies in the fibre π −1 [Ŝ], and equality in the first step says B is a nearest point of this fibre to A. But the nearest point of a fibre 1 1 to a full-rank A is unique—it is Ŝ 2 Pol(Ŝ 2 A), from the optimal unitary of Eq. (14) (see also Proposition 21 below)—and  is one such point; hence B = Â. Remark 19 (Projection on the total manifold can be expansive). Unlike with BW projection in the base manifold, the Frobenius projection of Corollary 18 is not non-expansive: already √ for d = 1 and α = 1, the points ±ε project to ±1, an expansion by the factor 1/ε. The √ culprit is the floor—{B ∶ σmin (B) ≥ α} excludes a neighborhood of the singular cone and is non-convex, whereas the ceiling alone is the operator-norm ball, whose projection (the upper clip) is non-expansive by convexity. This does not contradict Lemma 6: the expansion occurs along the fibre direction (±ε lie in the same fibre), which the quotient metric collapses. Moving to the BW distance is precisely what restores non-expansiveness. Theorem 20 (Projected BW-GD is projected Euclidean GD). Let f ∶ PD(d) → R be differentiable, g ∶= f ○ π, and (ηt )t≥0 a sequence of step sizes. Let (St )t≥0 be generated by Projected BW-GD from S0 ∈ [αI, βI], St+1 = clipα,β [ [I − 2ηt ∇f (St )] St [I − 2ηt ∇f (St )] ],

(83)

and (At )t≥0 by singular-value clipped Euclidean GD from any A0 ∈ π −1 [S0 ], √ √ [At − ηt ∇g(At )] . At+1 = clipsv α, β

(84)

Then π[At ] = St for every t ≥ 0. Consequently one may propagate (At )t≥0 alone, for 2 −3/2 ˆ f ∈ {f1 , f2 }, halting once ∥∇g(A δ certifies f (St ) − f (S⋆ ) ≤ δ (Corollary 12), t )∥F ≤ 2κ and return St = π(At ) only upon termination. Proof. Induct on t; the case t = 0 holds by choice of A0 . Suppose π[At ] = St . By Lemma 5 the unprojected step A′t ∶= At − ηt ∇g(At ) satisfies π[A′t ] = [I − 2ηt ∇f (St )]St [I − 2ηt ∇f (St )], the argument of the clip in Eq. (83). Singular-value clipping upstairs implements eigenvalue 30

√ √ = clip clipping downstairs by the intertwining identity π ○ clipsv α,β ○ π of Corollary 18. α, β Hence π[At+1 ] = St+1 .

For the invariant projection problem (f = f2 ), both recursions are run with the twirled ˆ in place of ∇g in Eq. (84)—and the gradients—Φ(∇f2 ) in place of ∇f in Eq. (83) and ∇g theorem holds verbatim: the proof uses only that the step has the form covered by Lemma 5, which accepts any Hermitian X, and that clipping intertwines with π, neither of which is affected by the twirl. This perspective does not lead to a computational advantage: evaluating ∇g(At ) = 2∇f (St )At still passes through St = At A†t and its square roots. Its value is structural, in two respects. First, Theorem 20 exhibits the algorithm as projected Euclidean gradient descent over the non-convex spectral set π −1 [[αI, βI]], the flat-space description under which the PL machinery of Section 4 applies. Second, the gradient step acquires a closed form with a geometric meaning—it replaces the current square root by the weighted average of the optimally paired square roots of the references—which is exactly the form in which the unit-step descent lemma becomes the variance identity of Appendix B; the following proposition formalizes it. Write Pol(X) for the unitary polar factor of an invertible X, so 1 X = Pol(X) ∣X∣ with ∣X∣ ≡ (X † X) 2 . Proposition 21. Let R ∈ PD(d), f2 = 12 B(R, ⋅), and g2 = f2 ○ π. For every A ∈ GL(d), 1

1

A − ∇g2 (A) = R 2 U = argmin ∥A − B∥2F ;

U ≡ Pol(R 2 A).

B∈π −1 [R]

(85)

1

Proof. We first establish the second equality. It suffices to show (A, R 2 U ) is an optimal 1 1 pairing for (S, R), S ≡ AA† , i.e. ⟨R 2 U, A⟩ = F(R, S). From U = Pol(R 2 A), √ √ 1 1 1 1 1 R 2 A = U ∣R 2 A∣ ⇒ U † R 2 A = A† RA ⇒ ⟨R 2 U, A⟩ = Tr[U † R 2 A] = Tr [ A† RA] = F(S, R). (86) Now consider the first equality. We claim ∇g2 (A) = 2∇f2 (S)A = [I − R#S −1 ]A = A − R 2 U, 1

(87)

for which we verify the identity R#S −1 = R 2 V S − 2 with V ≡ Pol(R 2 S 2 ). Since 1 1 1 1 1 1 1 1 1 1 1 (R 2 S 2 )† (R 2 S 2 ) = S 2 RS 2 , the polar decomposition reads R 2 S 2 = V (S 2 RS 2 ) 2 , so 1 1 1 1 1 (S 2 RS 2 ) 2 = V † R 2 S 2 and 1

1

1

1

S −1 #R = S − 2 (S 2 RS 2 ) 2 S − 2 = S − 2 V † R 2 , 1

1 2

which equals its adjoint R V S Pol(X)Ũ ,

− 21

1

1

1

1

1

1

(88) 1 2

(as it is Hermitian). Writing A = S Ũ and using Pol(X Ũ ) =

(R#S −1 )A = R 2 V S − 2 S 2 Ũ = R 2 Pol(R 2 S 2 Ũ ) = R 2 Pol(R 2 A) = R 2 U. 1

1

1

1

1

1

1

1

1

(89)

The proposition makes the total-manifold iteration fully explicit, and with it the reading promised in Section 3.3: the update averages optimally paired square roots of the ensemble elements. Applying Proposition 21 termwise to f1 , and composing with the twirl in the invariant case, the algorithms become n

1

1

√ √ [ ∑ wi R 2 Pol (R 2 At ) ], At+1 = clipsv i i α, β i=1

1 1 √ √ [Φ(R 2 Pol (R 2 At ) )], At+1 = clipsv α, β

31

(90)

for the barycenter and invariant problems respectively, the latter initialized at some A0 ∈ G. The barycenter iteration thus averages the square roots of the ensemble elements that form ˆ an optimal pairing with the current iterate. Writing A′t ≡ At − ∇g(A t ) for the argument of the clip—the unprojected step—the stopping rule of Theorem 20 reads ∥At − A′t ∥2F ≤ 2δκ−3/2 , which is equivalent to requiring the step-length be short.

E

The BRT Fixed-Point Algorithm as BW Gradient Descent on the Invariant Submanifold

In this appendix, we show that the fixed-point algorithm of Brahmachari, Rubboli, and Tomamichel [BRT25] is Bures–Wasserstein gradient descent restricted to the invariant submanifold, contextualizing it within Riemannian optimization. Beyond the contextualization, the payoff is in showing that the algorithm navigates the flat geometry of the invariant square-root subspace, which enables the seamless application of PL inequalities in Section 4. The optimization problem and unconstrained gradient. Given a unitary group U and a non-invariant R ∈ PD(d), the objective is to Bures-project R onto P ≡ comm(U) ∩ PD(d): 1 1 f2 (S) = B(R, S) and ∇f2 (S) = (I − S −1 #R). (91) 2 2 In general the Riemannian gradient ∇f2 (S) ∉ TS P ≅ H; to remain on the invariant submanifold, it must be orthogonally projected, with respect to the BW metric, onto H. The following lemma shows this projection is the twirl. Lemma 22 (Twirling is Projection). For all P ∈ P and X ∈ TP PD(d) ≅ Hr(d), the twirl Φ(X) is the orthogonal projection of X onto TP P with respect to ⟨⋅, ⋅⟩BW P : Φ(X) ∈ TP P ≅ H

and

⟨Φ(X), X − Φ(X)⟩BW = 0. P

(92)

Proof. Φ(X) ∈ H by definition. Applying Φ to the Lyapunov equation X = P LP (X) + LP (X)P , and using that P commutes with every element of U , yields Φ(X) = P Φ(LP (X)) + Φ(LP (X))P , hence Φ(LP (X)) = LP (Φ(X)). Consequently 1 1 = Tr[LP (Φ(X))(X−Φ(X))] = Tr[Φ(LP (X))(X−Φ(X))] = 0, (93) ⟨Φ(X), X−Φ(X)⟩BW P 2 2 since Φ is idempotent and self-adjoint with respect to the HS inner product. The invariant submanifold is totally geodesic. Lemma 22 identifies the Riemannian gradient of f2 restricted to P as the twirled gradient, but RGD on the submanifold requires the exponential map of P with the induced metric. It turns out one may simply use ExpBW S , because P is a totally geodesic submanifold. Note that it does not suffice to verify that the endpoint of the step lies in P: a curve may leave the submanifold and return to it, in which case it would not be a genuine geodesic of the submanifold. Recall that a Riemannian submanifold N ⊆ M is totally geodesic if every geodesic of M starting tangent to N remains in N throughout its interval of existence.

32

Lemma 23 (The invariant submanifold is totally geodesic). P is a totally geodesic submanifold of the BW manifold; consequently BW ExpP P = ExpP ∣T P ,

for all P ∈ P,

P

(94)

wherever either side is defined. BW Proof. Let P ∈ P and X ∈ H ∩ dom(ExpBW P ); recall X ∈ dom(ExpP ) if and only if I + LP (X) > 0 (Eq. (11)). For any t with I + tLP (X) > 0 (in particular all t ∈ [0, 1], since there I + tLP (X) = (1 − t)I + t[I + LP (X)] > 0) we have 1

1

2 2 γ(t) ≡ ExpBW P [tX] = π [(I + tLP (X)) P ] = π [(I + tLP (Φ(X))) P ] 1

1

= π [(I + tΦ(LP (X))) P 2 ] = π [[Φ(I + tLP (X))] P 2 ]

(95)

where we used, in order: Φ(X) = X for invariant X; Φ ○ LP = LP ○ Φ for invariant P (proof of Lemma 22); unitality of Φ; and closure of the subalgebra comm(U) under products. Combined with the domain characterization, γ(t) ∈ PD(d) ∩ comm(U) = P for every t in the interval of existence, which is precisely the statement that P is totally geodesic. Consequently the exponential maps agree on TP P. Deriving the submanifold update rule. By Lemmas 22 and 23, RGD on P with step size η is S ↦ [I − 2ηΦ(∇f2 (S))]S[I − 2ηΦ(∇f2 (S))], (96) with projected gradient Φ(∇f2 (S)) = 21 (I −Φ(S −1 #R)). Indeed, the step along the projected 1

2 Riemannian gradient is ExpBW S [−η Φ(∇f2 (S))] = π [[I − ηLP (Φ(∇f2 (S)))] S ], and since

LP ○ Φ = Φ ○ LP for invariant S (proof of Lemma 22) while ∇f2 = 21 LP (∇f2 ), we have LP (Φ(∇f2 (S))) = Φ(LP (∇f2 (S))) = 2Φ(∇f2 (S)), which yields the displayed update. At η = 1 the inner factor becomes Φ(S −1 #R), recovering the fixed-point map of Brahmachari, Rubboli, and Tomamichel [BRT25]: S ↦ Φ(S −1 #R) S Φ(S −1 #R) ≡ K(S).

(97)

Equivalence in the square-root manifold. The update gains further intuition upstairs: it is projected Euclidean GD on the invariant subspace G = GL(d) ∩ comm(U). Initializing at 1 1 1 an invariant square root S 2 ∈ G, the twirled lifted gradient is Φ(∇g2 (S 2 )) = 2Φ(∇f2 (S))S 2 1 (as S 2 commutes with U), and the projected Euclidean step A′ = S 2 − ηΦ(∇g2 (S 2 )) = [I − 2ηΦ(∇f2 (S))] S 2 1

1

1

maps under π exactly to Eq. (96).

33

(98)

Related documents

• Hessian Schatten-Norm Regularization for Linear Inverse Problems#235364
• Wasserstein Contraction of Coordinate Ascent Variational Inference#238608
• Covariance Shrinkage via Stochastic Interpolation#266177
• A Riemannian Inertial Adaptive Gradient Method with Momentum Restart#1117709
• On the Convergence of Stochastic Low-Rank Adaptation#1116380