Adaptive Bregman Proximal Stochastic Gradient with a Stabilized Barzilai--Borwein Step Size

Adaptive Bregman Proximal Stochastic Gradient with a Stabilized Barzilai–Borwein Step Size Chenhan Jin1 Shengze Xu1 Binghui Xie1 Kaiwen Zhou1 Fan JIA2 James Cheng1 Tieyong Zeng3,4∗ 1

The Chinese University of Hong Kong 3 Beijing Normal-Hong Kong Baptist University

2 4

University of Utah Guangzhou Nanfang College

arXiv:2608.12009v1 [math.OC] 12 Aug 2026

Abstract Bregman proximal stochastic gradient (BPSG) methods bring variance-reduced composite optimization to objectives whose geometry is poorly captured by Euclidean smoothness. Their performance, however, remains sensitive to the step size: raw stochastic curvature estimates can fluctuate sharply, whereas line searches add repeated proximal evaluations. We introduce Ada-BPSG, a line-search-free BPSG method that couples the SAGA gradient table with a stabilized Barzilai–Borwein (BB) candidate. A mediant aggregates incremental secant information so that nearly singular local ratios receive little weight, and an explicit safeguard translates the resulting curvature estimate into the bounded step-size sequence required for convergence. This design yields a direct analytical chain from relative smoothness and component-wise variance control to convergence in finite-dimensional normed spaces. We prove an O(n/K) ergodic rate for convex objectives, a restarted linear rate under relative quadratic growth, and an O(1/K) bound for a Bregman proximal residual in the nonconvex setting. On logistic regression and sparse nonnegative matrix factorization, Ada-BPSG combines low objective values with substantially less sensitivity to the initial step size than standard variance-reduced baselines, while avoiding line search.

1

Introduction

We study finite-sum composite optimization on a finite-dimensional normed vector space E: n

min F (x) := f (x) + h(x), x∈E

f (x) :=

1X fi (x), n

(1.1)

i=1

where [n] := {1, . . . , n} indexes the component losses fi : E → R, and h : E → (−∞, +∞] is a proper convex regularizer whose proximal subproblem is tractable. This model includes regularized empirical risk minimization, constrained learning, and structured estimation. When h = 0 and E = Rd , stochastic gradient descent (SGD) samples ik ∈ [n] and updates xk+1 = xk − ηk ∇fik (xk ) with step size ηk > 0. For a nonsmooth h, the Euclidean proximal stochastic-gradient update replaces the gradient step by   1 2 e ∥x − xk ∥2 + h(x) , (1.2) xk+1 = arg min ⟨∇k , x⟩ + x 2ηk e k is an unbiased estimator of the full gradient ∇f (xk ) conditional on the iteration where ∇ history. With a constant step size, the persistent variance of a plain sampled gradient can prevent the iterates from settling near a solution. Two complementary responses are therefore ∗

Corresponding author: tieyongzeng@bnbu.edu.cn

1

central to stochastic optimization: variance-reduced (VR) estimators, which make the gradient noise vanish, and adaptive step sizes, which respond to the local scale of the objective. The first family includes SAG, SAGA, SVRG, and SARAH [15, 23, 33, 36]; the second includes the Barzilai–Borwein (BB) rule, Polyak step sizes, AdaGrad, and Adam [5, 17, 26, 35]. Their combination has produced adaptive VR methods such as SVRG-BB, SAG-BB, AdaSVRG, AI-SARAH, ADASPIDER, AdProxGD, and 2AdaSPIDER [16, 21, 25, 31, 37, 40]. The Euclidean model behind (1.2) can nevertheless be a poor match for modern objectives. A globally Lipschitz Euclidean gradient is unavailable in applications such as neural-network training, quantum tomography, and matrix or tensor factorization [14, 20, 27]. Bregman geometry replaces the squared Euclidean distance by a divergence generated by a problem-adapted kernel ψ:   1 e xk+1 = arg min ⟨∇k , x⟩ + Dψ (x, xk ) + h(x) . (1.3) x ηk Section 3.1 formally defines the Bregman divergence Dψ . The kernel encodes the geometry of the domain, while relative smoothness replaces a global Euclidean Lipschitz-gradient condition [7, 8, 29]. This geometric extension sharpens the step-size problem. Bregman divergences are generally asymmetric, so Euclidean polarization identities no longer drive the descent proof, and stochastic secant ratios may vary markedly across sampled components. Diminishing schedules trade this instability for slow progress, while line searches add repeated objective or proximal evaluations. Our goal is to retain stochastic curvature adaptation without either cost.

1.1

Related work and positioning

Bregman projection, forward–backward, mirror-prox, accelerated Bregman proximal-gradient, and relative-smooth first-order methods establish the deterministic geometric foundation of this work [6–9, 11, 19, 29, 41]. In Euclidean geometry, optimized proximal-gradient schemes have also been applied to structured inverse problems, including ℓ1 - and total-variation-regularized non-blind image deblurring [44]. Their stochastic counterparts incorporate sampled gradients, primal–dual splitting, or variance reduction [34, 38, 42, 43]. In particular, BPSG-SAGA and BPSG-SARAH extend Bregman proximal optimization to variance-reduced estimators, and extrapolated BPSG uses line search to control the resulting updates [42, 43]. Adaptive stochastic methods approach the same stability question through the step size. AdaGrad and Adam accumulate gradient statistics [17, 26, 45]; universal and relative-smoothness schemes estimate local regularity in mirror-descent, primal–dual, or extragradient frameworks [1–4, 13, 24, 28]; and BB methods extract scalar curvature from secant pairs [5]. Stabilized BB and BB-VR variants improve this estimate in stochastic Euclidean problems [10, 30, 40, 46, 47]. Related powered and double-adaptive methods use BB or SPIDER-style estimators in Euclidean finite-sum and minimax settings [21, 22, 46]. Ada-BPSG connects these lines through one mechanism: the SAGA table supplies incremental secant pairs, a stabilized BB rule converts them into a curvature-sensitive candidate, and clipping converts that candidate into a provably stable Bregman proximal step. Table 1 places it among the closest Bregman proximal methods. This leads to the central question of the paper: Can a line-search-free BPSG method adapt to stochastic curvature while preserving convergence in non-Euclidean composite optimization?

1.2

Contributions

We answer this question with Ada-BPSG. Its contributions are threefold. • A stable adaptive step from SAGA information.

2

Table 1: Positioning of Ada-BPSG among Bregman proximal (variance-reduced) methods. All operate on a distance-generating kernel ψ. “LS-free” marks methods needing no line search; “Needs L̄, M ” lists the smoothness constants required to set the step or safeguard. Rates are ergodic orders in the iteration count K at fixed batch size; “lin. (RQG)” is the restarted linear rate under relative quadratic growth. A dash marks a regime the cited analysis does not cover: the stochastic Bregman baselines [42, 43] are established only in the nonconvex setting, so the convex ergodic and restarted-linear guarantees are, to our knowledge, new to Ada-BPSG. Crucially, our analysis is carried out in a general finite-dimensional normed space—using only the dual norm, with no Euclidean polarization or norm decomposition—so the matching O(n/K) order is obtained under strictly weaker geometric assumptions rather than by improving the rate. Over BPSG-SAGA, Ada-BPSG adds only two scalar accumulations per refreshed component. Method BPG [8] BPSG-SAGA [42, 43] BPSGE [43] Ada-BPSG (ours)

Oracle

Step-size rule

LS-free

Needs L̄, M

Convex

Nonconvex

full SAGA SAGA/SARAH SAGA

fixed 1/L̄ fixed line search safeguarded BB

Yes Yes No Yes

L̄ L̄, M L̄, M L̄, M

O(1/K) – – O(n/K); lin. (RQG)

O(1/K) O(1/K) O(1/K) O(1/K)

We aggregate recent component-wise secant pairs by a mediant, equivalently a denominatorweighted mean of local BB ratios. The weighting suppresses ratios associated with nearly singular curvature estimates. Clipping and a monotone update then turn this raw candidate into a bounded, line-search-free step-size sequence. For linear-prediction losses, the same quantities can be maintained with an O(n) scalar table instead of the general O(nd) SAGA table. • A convergence analysis aligned with the algorithm. Relative smoothness controls Bregman descent, local component smoothness controls the SAGA table variance, and the safeguard closes the resulting Lyapunov recursion. The entire argument is carried out in a general finite-dimensional normed space, relying only on the dual norm, dual-norm Young inequalities, and the three-point Bregman inequality: it never invokes the Euclidean polarization identity or the norm decomposability that Euclidean adaptive variance-reduced analyses rest on. This chain yields an O(n/K) ergodic rate for convex objectives, a restarted linear rate under relative quadratic growth, and an O(1/K) bound on the expected squared Bregman proximal residual for nonconvex objectives. Matching these Euclidean-SAGA orders under the weaker non-Hilbert geometry—rather than improving them—is the technical content of the analysis. The convex ergodic and restarted-linear guarantees moreover have no counterpart in the SAGA/SARAH-based stochastic Bregman methods we build on, which are analyzed only in the nonconvex regime [42, 43]. • Experiments that test the design mechanism. A curvature diagnostic shows why arithmetic averages of stochastic BB ratios can spike and how the mediant stabilizes the candidate. Across four logistic-regression datasets, Ada-BPSG reaches low objective gaps in fewer effective passes and remains substantially less sensitive to initialization than standard VR baselines. On a simplex-constrained Poisson inverse problem under the entropy kernel—an instance meeting every hypothesis of the analysis, in which variance reduction is genuinely needed—it exhibits the predicted convex rate in genuinely non-Euclidean geometry and, by adapting past the conservative worst-case safeguard, improves on every fixed-step baseline held to its a-priori guaranteed step by more than two orders of magnitude—a gain that carries over unchanged to a regionpooled Poisson unmixing of the real Samson hyperspectral scene; on sparse nonnegative matrix factorization it attains the lowest objective curves among the compared Bregman methods while avoiding the line search used by the accelerated baseline.

3

2

Algorithm

Ada-BPSG adds one curvature-adaptive scalar to the BPSG-SAGA update. The construction proceeds from the classical BB secant ratio, identifies why an arithmetic average of stochastic ratios remains unstable, and replaces that average by a safeguarded mediant computed from information already stored by SAGA.

2.1

Barzilai–Borwein step size

The BB method is a scalar quasi-Newton scheme: it approximates the local Hessian by ηk−1 I and chooses ηk to fit a secant equation. In the Euclidean setting, define the iterate and gradient differences sk := xk − xk−1 and yk := ∇f (xk ) − ∇f (xk−1 ). Minimizing either ∥ηk−1 sk − yk ∥2 or ∥sk − ηk yk ∥2 gives the two classical rules [5] ηkBB1 =

sTk sk sTk yk

ηkBB2 =

and

sTk yk . ykT yk

(2.1)

The first and second expressions correspond to the two residual problems, respectively. Both estimate inverse curvature along the latest displacement; our construction follows the BB2 orientation because SAGA directly supplies component-gradient differences.

2.2

Why stochastic BB ratios need stabilization

A BB ratio uses only one displacement and is therefore sensitive to sampled curvature. In the BB2 rule, a small ykT yk can produce an excessively large step, especially near a solution or along a low-curvature component. Denominator regularization and explicit bounds alleviate this behavior [10, 18, 30, 47], but an arithmetic mean over recent ratios can still be dominated by a single nearly singular denominator. Figure 1 isolates this failure mode on the ijcnn1 dataset [12]. For the regularized logistic component λ fi (x) = log(1 + exp(−bi aTi x)) + ∥x∥22 , 2 ai ∈ Rd is a feature vector and bi ∈ {−1, +1} is its label. The Hessian at epoch k is n

Hk =

1X pi (xk )ai aTi + λI, n

pi (x) :=

i=1

exp(−bi aTi x) , [1 + exp(−bi aTi x)]2

so pi (x) records the sample-dependent curvature weight. Even as these weights settle, the vanilla and arithmetic-mean stabilized BB rules in Figure 1b spike to 103 –104 . In contrast, the mediant candidate stays within the smallest and smoothest range among the plotted rules. The diagnostic motivates a division of labor that runs through the paper: the mediant stabilizes the empirical curvature estimate, and the safeguard enforces the bounds used by the convergence analysis.

2.3

Stabilized Barzilai–Borwein candidate

For positive denominators, the mediant of ratios satisfies P X bi ai a Pi i = P . i bi j bj bi i

Hence the mediant is a denominator-weighted arithmetic mean: a nearly singular denominator contributes little weight instead of dominating the aggregate.

4

stable

0.25 10

4

0.20 10 Step Size

0.15

0.10

0.05

10

10

0.00

10

3

BB1 BB2 BB1_STA BB2_STA OURS

2

1

0

1 2 3 4 5 7 9 11 13 15 17 19 21 23 25

5

Number of effective passes

10 15 20 Number of effective passes

(a)

25

(b)

Figure 1: Stochastic curvature and BB step sizes on ijcnn1, averaged over four runs. Left: distribution of the component curvature weights pi (xk ). Right: step-size trajectories over 30 effective passes; the mediant-based rule suppresses the spikes produced by vanilla and arithmetic-mean BB variants.

To apply this idea to SAGA, let φti denote the most recent point stored for component i at iteration t, and let Wk collect the component updates since the previous step-size refresh. Each record (i, t) ∈ Wk supplies the secant pair si,t := φt+1 − φti , i

t yi,t := ∇fi (φt+1 i ) − ∇fi (φi ).

Algorithm 1 maintains the two window sums X ξk := |⟨yi,t , si,t ⟩|,

δk :=

(i,t)∈Wk

X

∥yi,t ∥2∗ ,

(i,t)∈Wk

where ⟨·, ·⟩ is the duality pairing and ∥ · ∥∗ is the dual norm. The raw adaptive candidate is  ξk 1 , ξk > 0, ηbkBB = α δk + θξk (2.2)  ηk−1 , ξk = 0. The absolute value makes the numerator nonnegative when a nonconvex component produces a negative secant pairing. The buffer θ > 0 prevents denominator degeneracy and bounds the raw candidate by 1/(αθ); α > 0 controls its scale. In Euclidean space, the construction reduces to an aggregated BB2-type ratio. The safeguard maps this curvature estimate to qk = clip[ηmin ,ηmax ] (b ηkBB ) := min{ηmax , max{ηmin , ηbkBB }},

(2.3)

where 0 < ηmin ⩽ ηmax . Algorithm 1 then uses the monotone update ηk = max{ηk−1 , qk }, which preserves the interval and supplies the monotonicity needed by the convex Bregman telescope. Thus the raw ratio governs adaptation within the interval, while the safeguard supplies the stability invariant used in Section 3.

5

2.4

Main method

Algorithm 1 integrates the candidate into BPSG-SAGA. A batch Bk refreshes the corresponding e k , and the same table table entries, the standard SAGA correction forms the unbiased estimator ∇ differences update ξk and δk . These accumulators are reset every m iterations after refreshing the step size, so adaptation requires only scalar reductions beyond the gradient-table operations. The parameters have separate, explicit roles: m is the refresh frequency, θ regularizes the denominator, α rescales the raw ratio, and [ηmin , ηmax ] defines the stability range. The algorithm adapts within this range instead of following a prescribed diminishing schedule. When E = Rd and ψ(x) = 12 ∥x∥22 , the Bregman update becomes the Euclidean proximal update (1.2); we call this specialization SAGA-BB in the experiments.

2.5

Efficient implementation

The general SAGA table stores one d-dimensional gradient per component and therefore costs O(nd) memory. Linear-prediction losses have the form fi (x) = Ψi (⟨ai , x⟩), where ai ∈ Rd is a data vector and Ψi : R → R is a scalar loss [48]. In this case, storing only ⟨ai , φki ⟩ recovers both the gradient-table correction and the secant sums in (2.2). The implementation in Appendix A.2 consequently reduces the table memory to O(n) without changing the Ada-BPSG update. This reduction follows from the factorization ∇fi (φki ) = Ψ′i (⟨ai , φki ⟩)ai , where Ψ′i denotes the derivative of the scalar loss. When component i is refreshed, its old and new scalar predictions determine both the SAGA gradient correction and its contributions to ξk and δk ; the feature norm ∥ai ∥∗ can be precomputed. Ada-BPSG therefore needs neither a separate curvature table nor a history of past iterates. The same reuse preserves the computational profile of the underlying BPSG-SAGA iteration. For a mini-batch of size b, step-size adaptation adds two scalar accumulations per refreshed component and a constant-cost safeguard whenever the step size is updated. It introduces no extra component-gradient evaluations or proximal subproblems. This is the implementation counterpart of the method’s central design: the stochastic correction and the curvature estimate use the same local information, while the safeguard converts that information into the stable step-size sequence required by the analysis.

6

Algorithm 1 Ada-BPSG 1: Input: max iterations K, batch-size b ∈ {1, . . . , n}, initial point x0 , kernel ψ

initial step-size η0 ∈ [ηmin , ηmax ], update frequency m ⩾ 1 BB scaling α > 0, BB buffer θ > 0, safeguard bounds 0 < ηmin ⩽ ηmax 2: Initialize: P φ0i = x0 , gradient table {∇fi (φ0i )}ni=1 n 1 µ0 = n i=1 ∇fi (φ0i ), ξ0 = 0, and δ0 = 0 3: for k = 0, 1, . . . , K − 1 do 4: if k > 0 then 5: ηk = ηk−1 6: end if 7: if k > 0 ( and k mod m = 0 then ξk 1 , ξk > 0, 8: ηbk = α δk +θξk ηk−1 , ξk = 0, 9: ηk = max{ηk−1 , clip[ηmin ,ηmax ] (b ηk )} 10: ξk = 0, δk = 0 11: end if 12: Randomly pick a batch Bk uniformly without replacement 13: For i ∈ Bk , set φk+1 = xk and update ∇fi (φk+1 ) i i k+1 14: For i ∈ / Bk , set φi = φki
k+1 ek = 1 P 15: ∇ ) − ∇fi (φki ) + µk i (φi i∈Bk ∇f b o n e k , x⟩ + 1 Dψ (x, xk ) 16: xk+1 = arg minx h(x) + ⟨∇ 17: 18: 19: 20:

ηk k+1 k Compute yi = ∇fi (φi ) − ∇fi (φi ) for i ∈ Bk k+1 k Compute si = Pφi − φi for i ∈ Bk

ξk+1 = ξk + i∈Bk |⟨yi , si ⟩| P δk+1 = δk + i∈Bk ∥yi ∥2∗
P ) − ∇fi (φki ) 21: µk+1 = µk + n1 i∈Bk ∇fi (φk+1 i 22: end for 23: return xK

3

Convergence Analysis

The analysis follows the same chain as the algorithm. Relative smoothness converts the Bregman proximal update into descent, component smoothness controls the SAGA estimation error, and the safeguard keeps the adaptive step within the range that balances these two terms. We first define the geometry and the step-size invariant, then establish variance bounds for the convex and nonconvex regimes, and finally derive the convergence rates.

3.1

Preliminaries

Normed-space notation. Let E be a finite-dimensional real normed vector space with norm ∥ · ∥, and let E∗ be its dual space, the space of continuous linear functionals on E. For y ∈ E∗ and x ∈ E, ⟨y, x⟩ denotes their duality pairing, and ∥y∥∗ := max{⟨y, x⟩ : ∥x∥ ⩽ 1}. This expression defines the dual norm ∥ · ∥∗ . The filtration Fk := σ(B0 , . . . , Bk−1 ) contains the mini-batches sampled before iteration k, and Ek [·] := E[· | Fk ] denotes conditional expectation given this history. Since E is finite-dimensional, it is automatically complete, hence Banach; closed bounded subsets are compact, all norms are equivalent, and E∗ has the same dimension. Every existence and boundedness fact used below—solvability of the proximal subproblem and existence of x∗ ∈ D—follows from dim E < ∞ together with the strong convexity of ψ, so no completeness or sequential-compactness assumption beyond finite-dimensionality is invoked. In fact the convergence results are non-asymptotic, finite-horizon bounds (Theorems 3.6–3.11) and use no 7

limit-point or subsequence argument at all. The generality relative to the Euclidean setting is therefore geometric—the norm need not derive from an inner product, so no polarization identity or norm decomposition is available—rather than topological. Bregman geometry. A distance-generating function (DGF), or Bregman kernel, is a continuously differentiable function ψ : E → R that is 1-strongly convex with respect to ∥ · ∥: 1 ψ(x) − ψ(y) − ⟨∇ψ(y), x − y⟩ ⩾ ∥x − y∥2 . 2 Its Bregman divergence is the first-order linearization error Dψ (x, y) := ψ(x) − ψ(y) − ⟨∇ψ(y), x − y⟩. Strong convexity gives Dψ (x, y) ⩾ 21 ∥x − y∥2 , which connects Bregman descent to the norm differences used in the SAGA variance bounds. Given a direction G ∈ E∗ and step size η > 0, define the Bregman proximal mapping by   1 Proxψ (x, G) := arg min ⟨G, u⟩ + D (u, x) + h(u) ψ h,η u∈E η whenever the minimizer exists. Assumption 3.4 ensures that this mapping is single-valued along the algorithmic trajectory. Definition 3.1 (Smooth adaptability [8, 29, 43]). Let X ⊂ E be nonempty, convex, and open, and let f : X → (−∞, +∞] be proper, lower semicontinuous, and continuously differentiable on X . Given a DGF ψ, the pair (f, ψ) is (L̄, L)-smooth adaptable on X if constants L̄ > 0 and ¯ L ⩾ 0 satisfy, for every x, y ∈ X , ¯ −LDψ (x, y) ⩽ f (x) − f (y) − ⟨∇f (y), x − y⟩ ⩽ L̄Dψ (x, y). ¯ The upper inequality is the relative-smoothness bound used for descent; the lower inequality controls negative curvature. The convergence proofs below use the upper constant L̄ on a deterministic region D ⊂ X . Example 3.2 (A genuinely non-Euclidean instance). Let E = Rd carry the norm ∥ · ∥ = ∥ · ∥1 , so that the dual of the probability simplex P norm is ∥ · ∥∗ = ∥ · ∥∞ , and let h = ι∆ be the indicator P ∆ = {x ⩾ 0 : i xi = 1}. The negative-entropy kernel ψ(x) = i xi log xi is 1-strongly convex with respectP to ∥ · ∥1 by Pinsker’s inequality, and its divergence is the Kullback–Leibler divergence Dψ (x, y) = P i xi log(xi /yi ), which is asymmetric and admits no polarization identity. For a finite sum f = n1 i fi with bounded Hessians on ∆, the pair (f, ψ) is relatively smooth; the SAGA e k and the accumulators ξk = P |⟨yi,t , si,t ⟩| and δk = P ∥yi,t ∥2 are then measured correction ∇ ∞ in ℓ∞ , and the proximal step is the KL projection onto ∆, i.e. a softmax renormalization. Every constant in Theorems 3.6–3.11 is expressed through the ℓ1 /ℓ∞ pairing, with σb ⩾ 1 typically strict; the Euclidean specialization ψ = 12 ∥ · ∥22 used in Section 4 is the degenerate case ∥ · ∥ = ∥ · ∥∗ = ∥ · ∥2 , σb = 1. Standard Euclidean-underlying stochastic Bregman analyses, which measure gradients and variance in ℓ2 , do not cover this instance directly.

3.2

Safeguarded step-size framework

The candidate ηbkBB , its clipped value qk , and the monotone update ηk = max{ηk−1 , qk } are defined in (2.2)–(2.3). Because they use only records collected before Bk is sampled, ηk is Fk -measurable. The update maintains the deterministic invariant ηmin ⩽ ηk ⩽ ηmax ,

ηk ⩾ ηk−1

for all k.

The lower bound sets the scale in the convex rate, the upper bound absorbs relative smoothness and gradient-estimation error, and monotonicity telescopes the Bregman distance in the convex proof. The nonconvex argument uses only the upper bound. 8

Remark 3.3 (Cost of monotonicity). The nondecreasing rule ηk = max{ηk−1 , qk } is precisely what makes ak = 1/ηk nonincreasing and thus lets the Bregman distance telescope in Theorem 3.6. Its price is that the step may rise but never fall: starting from ηmin , it ramps toward the clipped curvature estimate and then holds. Overshoot is controlled not by retreat but by the fixed cap ηmax , so the adaptation is confined to how fast and how far the step climbs within [ηmin , ηmax ]. The nonconvex analysis (Theorem 3.11) does not use monotonicity and relies only on the upper bound, so there the step is free to vary in both directions inside the interval.

3.3

Assumptions

Assumption 3.4. Items 1–5 apply to all results; Item 6 is used only in the convex analysis. 1. Geometry and prox. The function ψ is continuously differentiable and 1-strongly convex with respect to ∥ · ∥. The function h : E → (−∞, +∞] is proper, closed, and convex. For every x ∈ D, G ∈ E∗ , and η ∈ [ηmin , ηmax ], the Bregman proximal subproblem has a unique solution. 2. Deterministic bounded region. There exists a deterministic bounded convex set D ⊂ X such that xk , φki ∈ D almost surely for all k and i. In the convex results, an optimal solution x∗ also belongs to D. 3. Sampling. Conditional on Fk , the mini-batch is sampled uniformly from [n] with size b and without replacement. Write I1 , . . . , Ib for the sampled indices. A deterministic mini-batch second-moment constant σb ⩾ 1 satisfies, for any deterministic v1 , . . . , vn ∈ E∗ 1 P with v̄ = n i vi ,  2 b n X 1 σb X E (vIj − v̄)  ⩽ ∥vi − v̄∥2∗ . b bn j=1

∗

i=1

In the Euclidean norm, one may take σb = 1. Conditional on Fk , the table vectors are fixed, so this inequality applies directly to the SAGA estimator. 4. Relative smoothness. The pair (f, ψ) is (L̄, L)-smooth adaptable on D; only the upper ¯ bound with L̄ is used in the convergence proofs. This is a relative smoothness condition for the aggregate differentiable term f with respect to the chosen kernel ψ. 5. Component local smoothness. Each ∇fi is Lipschitz on D with deterministic constant M , i.e., ∥∇fi (x) − ∇fi (y)∥∗ ⩽ M ∥x − y∥, x, y ∈ D, i ∈ [n]. 6. Convex self-bounding condition. For the convex results, each fi is convex on D and satisfies ∥∇fi (x) − ∇fi (y)∥2∗ ⩽ 2M (fi (x) − fi (y) − ⟨∇fi (y), x − y⟩) ,

x, y ∈ D.

Role of the assumptions. The deterministic set D localizes the component smoothness needed for variance reduction. A concrete verifiable sufficient condition is that dom h is compact and the kernel ψ is continuous on dom h: the Bregman proximal subproblem then minimizes a lower-semicontinuous function over a compact feasible set, so every iterate xk and every table point φki stays in D := dom h almost surely, and M is the finite Lipschitz constant of the component gradients on that compact set. This makes Item 2 an assumption on the problem data rather than on the trajectory. The sparse-NMF experiment of Section 4 does not satisfy this out of the box, because the objective is invariant under the rescaling (U, V ) 7→ (tU, t−1 V ) and its sublevel sets are therefore unbounded; appending an explicit box 0 ≤ U, V ≤ U to the 9

cardinality and nonnegativity constraints restores a compact D while leaving the optimal value unchanged for U large enough, at the cost of enforcing the box in the proximal step. Fixing D before the run keeps M deterministic. This local component condition serves a different purpose from relative smoothness: L̄ controls descent of the aggregate f in Bregman geometry, whereas M controls differences between stored component gradients. For convex components, Item 6 converts those gradient differences into the linearization errors used by the objective-gap Lyapunov function. Convex Euclidean M -smooth functions satisfy this self-bounding inequality automatically. The nonconvex analysis instead tracks distances between the current iterate and the SAGA table points, and therefore uses only Items 1–5.

3.4

Variance bounds

For the convex analysis, let x∗ be an optimal solution in D and choose ρ∗ ∈ ∂h(x∗ ) satisfying the first-order condition ∇f (x∗ ) + ρ∗ = 0. Define ∆i (x) := fi (x) − fi (x∗ ) − ⟨∇fi (x∗ ), x − x∗ ⟩,

H(x) := h(x) − h(x∗ ) − ⟨ρ∗ , x − x∗ ⟩.

Here ∆i (x)Pis the linearization error of component fi , H(x) is the corresponding error of h, and ∆(x) := n1 i ∆i (x). Convexity makes all three quantities nonnegative; the optimality condition also gives F (x) − F (x∗ ) = ∆(x) + H(x). This identity links the SAGA variance to the objective gap. Lemma 3.5 (Convex SAGA variance). Under Assumption 3.4, the SAGA estimator in Algorithm 1 is conditionally unbiased, e k ] = ∇f (xk ) Ek [∇ and its conditional variance in the convex setting satisfies " # n h i 16σ M X 1 b e k − ∇f (xk )∥2∗ ⩽ Ek ∥∇ ∆(xk ) + ∆i (φki ) . b n i=1

3.5

Convex analysis

Lemma 3.5 makes the stochastic error proportional to the current and stored objective gaps. Combining this estimate with the three-point Bregman proximal inequality yields a Lyapunov decrease in which the safeguard absorbs the variance term. Theorem 3.6. Suppose Assumption 3.4 holds in the convex setting, and the safeguarded step-sizes satisfy 1 ηmin ⩽ ηk ⩽ ηmax , ηk ⩾ ηk−1 , ηmax ⩽ . L̄ + 32σb M/b P Let x̄K = K1 K k=1 xk . Then    2 n 1 1 ∗ ∗ ∗ E[F (x̄K )] − F (x ) ⩽ + (F (x0 ) − F (x )) + Dψ (x , x0 ) . K 4b 2 ηmin Theorem 3.6 therefore gives an O(n/K) ergodic rate for fixed b. Its dependence on ηmin records the cost of the lower safeguard, while the admissible ηmax couples relative smoothness and SAGA variance. The BB candidate remains free to adapt anywhere inside this interval.

10

Remark 3.7 (What the guarantee does and does not claim). The admissible cap ηmax ≤ (L̄ + 32σb M/b)−1 encodes exactly the smoothness constants L̄ and M that a well-tuned fixed step would use. Consequently Theorem 3.6 does not improve on the O(n/K) rate of fixed-step BPSGSAGA; it certifies that adapting within [ηmin , ηmax ] does not degrade that rate. Ada-BPSG is thus adaptive within a safe interval whose upper end still depends on problem constants, not free of them; the empirical payoff (Section 4) is robustness to the initial step size rather than a faster −1 worst-case rate. The leading constant is moreover governed by ηmin Dψ (x∗ , x0 ), so a conservative ηmin enlarges the bound: the BB candidate removes the need to guess where in [ηmin , ηmax ] the step should sit, but not the need to set the interval from L̄ and M . The analysis is nonetheless conducted entirely in the general normed space E: it uses only the dual norm, dual-norm Young inequalities, and the three-point identity of Lemma A.1, and at no point invokes the Euclidean polarization identity 2⟨a, b⟩ = ∥a∥2 + ∥b∥2 − ∥a − b∥2 or the norm decomposability on which the Euclidean adaptive-BB analyses rely. The price of dropping them is explicit rather than hidden: the mini-batch second-moment constant obeys σb ⩾ 1 with equality only in the Euclidean norm (Assumption 3.4, Item 3), Pand the variance P bound of Lemma 3.5 carries a factor 4 where the Hilbert variance identity i ∥ai − ā∥2 = i ∥ai ∥2 − n∥ā∥2 would give 1. Recovering the Euclidean O(n/K) order in this non-Hilbert geometry, rather than improving it, is where the difficulty lies, so the contribution is added generality at no loss of rate. Relative quadratic growth means that the objective gap controls the Bregman distance to a solution. Under this condition, restarting the averaged convex guarantee turns a constant-factor reduction per round into geometric convergence. Theorem 3.8. In addition to the assumptions of Theorem 3.6, suppose F satisfies the relative quadratic growth condition F (x) − F (x∗ ) ⩾ µDψ (x∗ , x), Let Ab :=

x ∈ D.

n 1 + . 4b 2

If each restart round reinitializes the SAGA table at x(t) , starts or continues with a step-size in [ηmin , ηmax ] satisfying the monotonicity condition within the round, and runs Theorem 3.6 for   1 K ⩾ 4 Ab + µηmin iterations and sets x(t+1) to the averaged iterate of round t, then after τ restart rounds,   E[F (x(τ ) )] − F (x∗ ) ⩽ 2−τ F (x(0) ) − F (x∗ ) . The same condition converts the objective contraction into a Bregman distance bound:    2−τ  F (x(0) ) − F (x∗ ) . E Dψ (x∗ , x(τ ) ) ⩽ µ If the kernel is 1-strongly convex, this further bounds the squared norm distance by E∥x(τ ) − x∗ ∥2 ⩽ 2E[Dψ (x∗ , x(τ ) )] whenever the displayed Bregman distance is evaluated on D.

3.6

Nonconvex analysis

We next allow f to be nonconvex and assume that F has a finite lower bound Finf . Because an objective gap alone does not characterize stationarity, we use a Bregman proximal gradient mapping. 11

Definition 3.9 (Bregman proximal gradient mapping). For η > 0, define Gη (x) =

x − Proxψ h,η (x, ∇f (x)) η

.

Let x ek+1 = Proxψ h,ηk (xk , ∇f (xk )). The full-gradient proximal point x ek+1 removes the SAGA estimation error from the stationarity measure. We normalize its displacement by the fixed upper safeguard: Rk :=

xk − x ek+1 . ηmax

This normalization yields a common scale across adaptive iterations. The residual Rk vanishes exactly when the full-gradient Bregman proximal step leaves xk unchanged. Lemma 3.10 (Nonconvex SAGA variance). Under Items 1–5 of Assumption 3.4, n h i 2 X e k − ∇f (xk )∥2∗ ⩽ 4σb M ∥xk − φki ∥2 . Ek ∥∇ bn i=1

Theorem 3.11. Suppose the common conditions in Assumption 3.4 hold, f may be nonconvex, and F is bounded below by Finf > −∞. Let Γb :=

8σb M 2 n(2n − b) . b3

If Γb ηmax + L̄ ⩽

1 ηmax

,

then for R drawn independently and uniformly from {0, . . . , K − 1},   2(F (x0 ) − Finf ) E ∥RR ∥2 ⩽ . Kηmax Theorem 3.11 gives an O(1/K) bound on the expected squared residual. The quantity Γb captures how mini-batch sampling and table refresh interact with local component smoothness, while the upper safeguard balances this variance term against relative smoothness. Because 3/2 /(nM )), −1 forces η Γb = O(M 2 n2 /b3 ), the stability condition Γb ηmax + L̄ ≤ ηmax max = O(b so the residual bound 2(F (x0 ) − Finf )/(Kηmax ) inherits the O(n)-type gradient complexity √ characteristic of SAGA rather than the improved O( n) dependence available to SARAH- or SPIDER-type recursive estimators. We nonetheless build on SAGA because its table supplies the component secant pairs (si,t , yi,t ) directly, which is what the stabilized candidate consumes; porting the same candidate to a memory-light recursive estimator is left to future work. Section 4 examines how the stabilized candidate behaves with practically chosen safeguards and kernels.

4

Numerical Experiments

The experiments probe the mechanisms developed above: whether the adaptive SAGA step improves gradient efficiency, whether stabilization reduces sensitivity to the initial scale, whether the convex O(n/K) guarantee is realized inside a genuinely non-Euclidean instance that meets every hypothesis of Assumption 3.4, whether the same behavior transfers from synthetic to real measured data, and whether the design remains effective on a harder non-Euclidean problem that falls outside the theory. We use E = Rd throughout and vary the distance-generating function ψ with the problem geometry. 12

4.1

Regularized binary classification

We first choose the Euclidean kernel ψ(x) = 12 ∥x∥22 and solve ℓ2 -regularized logistic regression, n

1X λ F (x) = log(1 + exp(−bi aTi x)) + ∥x∥22 , n 2 i=1

λ=

1 , n

where ai ∈ Rd and bi ∈ {−1, +1} are the features and label of sample i. Under this kernel, Ada-BPSG specializes to SAGA-BB. We compare it with SAGA, SVRG, SVRG-BB, loopless SVRG, and SARAH on mushrooms, ijcnn1, w8a, and covtype from LibSVM [12]. Every method uses the same tuning grid: the initial step size lies in {10−3 , 10−2 , 10−1 , 1, 10, 100} and the batch size lies in {1, 8, 16, 64}. For methods with an epoch length or refresh interval, we set it to m = n/b. Figure 2 connects the curvature diagnostic in Figure 1 to optimization performance. SAGABB reaches the smallest plotted objective gaps with the fewest effective passes on all four datasets. The gain is especially pronounced on mushrooms, w8a, and covtype, where its curve separates early from both SAGA and the SVRG-family baselines. Thus the stabilized secant information does more than prevent isolated step-size spikes: it improves the gradient efficiency of the underlying SAGA update. The initial-scale experiment in Figure 3 fixes b = 8 and measures the final gradient norm after 120 effective passes. The vanilla VR methods vary by several orders of magnitude across the grid and can become unstable at its extremes. SVRG-BB reduces this dependence on some datasets but remains uneven on w8a. SAGA-BB maintains a low gradient norm across the entire grid on every dataset. This behavior matches the algorithmic design: the initial value governs only the iterations before the first refresh, after which the mediant candidate adapts within the safeguard interval. Appendix B verifies the same pattern across additional batch sizes and with Huber loss.

4.2

Simplex-constrained Poisson inverse problem

The second experiment moves to a non-Euclidean kernel while staying strictly inside the theory, and is chosen so that variance reduction is genuinely needed. We recover mixing weights on the probability simplex from count observations, n  1 X min F (x) = ⟨ai , x⟩ − bi log⟨ai , x⟩ , x∈∆ n

n o P x = 1 , ∆ = x ∈ Rd : x ⩾ 0, j j

(4.1)

i=1

where the rows ai > 0 are known nonnegative dictionary atoms and bi ∼ Poisson(⟨ai , x⋆ ⟩) are counts generated by a planted simplex vector x⋆ . Each fi (x) = Ψi (⟨ai , x⟩) with Ψi (u) = u−bi log u is a convex linear-prediction loss, so the O(n) scalar-table implementation of Appendix A.2 applies verbatim, and h = ι∆ is the indicator of the simplex. Following Example 3.2 P we equip d E = R with the ℓ1 norm, its dual ℓ∞ , and the negative-entropy kernel ψ(x) = j xj log xj , which is 1-strongly convex with respect to ∥ · ∥1 ; the Bregman proximal step of Algorithm 1 then has the closed form of an entropic mirror update, e k ]j ) xk,j exp(−ηk [∇ xk+1,j = P , e l xk,l exp(−ηk [∇k ]l ) and the accumulators ξk , δk are formed in the ℓ∞ dual norm. Because the atoms are strictly positive, the prediction obeys ⟨ai , x⟩ ⩾ p := mini minj aij > 0 on the simplex, which keeps the loss away from the singularity of log and makes every constant in Assumption 3.4 finite: ∆ is compact and convex, so the deterministic region D = ∆ is available a priori; each Ψi is convex with Ψ′′i (u) = bi /u2 , giving the component-gradient bound M = (maxi bi /p 2 ) maxi ∥ai ∥2∞ 13

ijcnn1

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−4

10−6

10−8

10−10

10−12

−2

10−4

10−6

10−8

10−10

0

20

40

60

80

100

10−12

120

0

10

20

30

40

50

60

70

10−2

10−4

10−6

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−8

10−12

80

0

Number of effective passes

(a)

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

100

10−10

Number of effective passes

covtype

102

100

Objective minus best

10

Objective minus best

Objective minus best

10

−2

w8a

102 SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

100

Objective minus best

mushrooms 100

20

10−2

10−4

10−6

10−8

10−10

40

60

80

100

10−12

120

0

Number of effective passes

(b)

20

40

60

80

100

120

Number of effective passes

(c)

(d)

Figure 2: Gradient efficiency for regularized logistic regression. The horizontal axis counts componentgradient evaluations normalized by n, and the vertical axis is the objective gap to the best value attained in the corresponding experiment. Lines and shaded regions show the mean and standard deviation over five runs. All methods use their best configuration from the common tuning grid. mushrooms

ijcnn1

w8a

covtype

100

10−2

10−2

10−6

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−8

10−10

10

−3

10

−2

Gradient norm

Gradient norm

10−4

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−6

10−8

10

10

10

1

10

2

10−7

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless

10−11

10−12

0

10−5

10−4

10−6

10−8

10−9

10−10

−1

Gradient norm

10−3

10−4

Gradient norm

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless

100 10−1

10−2

10

−3

10

−2

10

−1

10

0

10

1

10

2

10−3

10−2

10−10

10−12 10−1

100

101

102

10−3

10−2

10−1

100

Initial step size

Initial step size

Initial step size

Initial step size

(a)

(b)

(c)

(d)

101

102

Figure 3: Sensitivity to the initial step size after 120 effective passes with batch size b = 8. Each part reports the final gradient norm; finite runs within the displayed range are shown.

of Items 5–6 and a finite relative-smoothness scale L̄ against the entropy kernel. Crucially, unlike a quadratic data term, the Poisson loss has no global Euclidean Lipschitz gradient, so constant-step SGD retains a variance floor while a variance-reduced step can drive the gap arbitrarily low: this is exactly the regime the method targets, and the convex rate of Theorem 3.6 is expected to hold in it. The geometry is non-Hilbert throughout, with σb > 1 in general. We use a controlled synthetic instance with d = 500, n = 5000, atoms ai with entries drawn uniformly in [0.05, 1.05], a planted interior signal x⋆ , and Poisson counts bi ; the batch size is b = 8 and the refresh interval is m = n/b. Because the problem is convex, we obtain the reference value F ⋆ from an independent high-accuracy full-batch entropic mirror descent with backtracking, so that no plotted stochastic method is scored against its own best iterate. We compare Ada-BPSG against deterministic BPG, against BPSG with SGD, SAGA, and SARAH estimators, and against the extrapolated line-search variant BPSGE. The decisive feature of this instance is that its worst-case component constant M is dominated by the log singularity and is therefore far larger than the curvature the iterates actually visit: here L̄ ≈ 18 but M ≈ 2.2 × 103 , whereas the local component smoothness at typical predictions is only Mloc ≈ 18. Figure 4 (a)–(b) therefore adopt an a-priori protocol: every fixed-step baseline uses the largest step its own analysis guarantees from these quantities—(L̄ + M )−1 for SGD, (L̄ + 32M/b)−1 for SAGA, and (L̄ + M/b)−1 for the recursive SARAH estimator—while Ada-BPSG adapts within a safeguard ηmax = (L̄ + Mloc /b)−1 built from the local constant, started from the lower safeguard ηmin . This local cap exceeds the worst-case safeguard (L̄ + 32σb M/b)−1 of Theorems 3.6–3.11, so it is a heuristic relaxation rather than a theory-faithful setting; to make the gap explicit we also plot “Ada-BPSG (theory cap)”, the same rule confined to the analysis cap with the global M . The complementary question—what a baseline achieves if its step is instead hand-tuned per instance—is answered separately in Figure 4 (c). Figure 4 reads along the same three axes as the classification study, now in Bregman geometry. In Figure 4 (a) the gap of Ada-BPSG decays steadily to 1.3 × 10−5 along the O(n/K) reference slope, exhibiting the convex-rate behavior predicted by Theorem 3.6 in a bona fide non-Euclidean instance. The fixed-step baselines, held to the steps their own theory licenses, remain near 10−2 : because the worst-case M is inflated by the log singularity, the guaranteed steps are two to

14

BPSG-SGD BPSG-SAGA BPSG-SARAH BPSGE Ada-BPSG (theory cap) Ada-BPSG (ours)

10 2

10 3 F(xK) F

10 3 F(xK) F

10 2

BPG BPSG-SGD BPSG-SAGA BPSG-SARAH BPSGE Ada-BPSG (theory cap) Ada-BPSG (ours) O(n/K) ref.

10 4

10 4 final F F

10 2

10 4

10 6 10 8 10 10

10 5

0

20

40 60 effective passes

10 5

80

(a) objective gap vs. effective passes

0.0

0.5

1.0 1.5 2.0 wall-clock time (s)

10 12

2.5

(b) objective gap vs. wall-clock time

BPSG-SAGA BPSG-SGD Ada-BPSG (ours)

10 3

10 2

10 1 100 initial step size

101

102

(c) sensitivity to the initial step

Figure 4: Simplex-constrained Poisson inverse problem under the entropy kernel. (a) Objective gap F (xK ) − F ⋆ against effective passes, with the O(n/K) reference slope of Theorem 3.6 overlaid; each fixed-step baseline uses the a-priori step guaranteed by its own analysis, and “Ada-BPSG (theory cap)” is the adaptive rule confined to the worst-case safeguard, whereas “Ada-BPSG (ours)” relaxes the cap to the local component constant. (b) The same gap against wall-clock time. (c) Final objective gap after a fixed budget as the initial step sweeps the grid {10−3 , . . . , 102 }: for the fixed-step methods the swept value is the step, while for Ada-BPSG it is only the starting point of the adaptive rule. Lines and shaded regions show the mean and standard deviation over five runs at batch size b = 8.

four orders of magnitude below the useful curvature scale, and the methods barely move in 80 passes. The “theory cap” run of Ada-BPSG, confined to the same worst-case safeguard, sits with them—so on this instance the analysis-faithful step is itself conservative, and the speedup comes precisely from adapting beyond it. This is the intended behavior of the stabilized rule: its secant candidate climbs from ηmin to the local safeguard, improving on the best a-priori baseline by more than two orders of magnitude while reading only the SAGA information it already maintains, and it does so as a heuristic relaxation of the worst-case cap, as is standard for BB-type steps. Figure 4 (b) shows the same separation in wall-clock time: Ada-BPSG reaches 10−5 in a few seconds, each step-size refresh costing only two scalar accumulations and no line search. Figure 4 (c) explains what tuning would buy the baselines and why adaptivity still wins in practice: the final gap of Ada-BPSG is essentially flat across the entire initial-step grid, whereas fixed-step SAGA is sharply V-shaped—it can reach a much lower gap, but only in a narrow band around a single step several orders of magnitude above the guaranteed cap, and it diverges beyond it—while SGD is pinned at its variance floor. Ada-BPSG thus recovers the accuracy of a hand-tuned step without any tuning and without the attendant fragility. Together Figure 4 exhibits the predicted O(n/K) decay and the practical value of adapting past the worst-case safeguard, complementing the harder out-of-theory test that follows.

4.3

Real-data confirmation: hyperspectral unmixing

The synthetic study isolates the mechanism on controlled data; we now verify that it survives on measured observations, with the same estimator, safeguards, and hyperparameters. We use the Samson scene [49], 95 × 95 pixels over L = 156 contiguous bands with three reference materials (soil, tree, water). To assemble a large finite sum on the simplex, we pool a spatially contiguous, spectrally homogeneous patch R under a shared-abundance model: all pixels in R are taken to share one abundance x ∈ ∆d over a dictionary D ∈ RL×d of measured spectra. Each material + contributes an endmember bundle of K = 6 representative spectra—cluster centers of the pixels it dominates, the standard construction for spectral variability [39]—so d = 18, and summing the weights within a bundle recovers a physical material abundance. Every pixel–band pair (p, ℓ) with p ∈ R is one component, with atom ai = Dℓ,: ≥ 0 and count bi = Ypℓ ≥ 0 the measured radiance, giving n  1 X min ⟨ai , x⟩ − bi log⟨ai , x⟩ , n = |R| · L, x∈∆ n i=1

15

the same generalized Kullback–Leibler (Poisson) fidelity as (4.1), now with real atoms and radiances. We select the 10×10 patch that is simultaneously most mixed and most homogeneous— mean abundance (0.37, 0.53, 0.10), a genuinely interior target—which yields n = 15 600. Exactly as in the synthetic case, the worst-case component constant is inflated by measured reflectances that dip to p ≈ 8.6 × 10−3 in some bands: here L̄ ≈ 68 and M ≈ 5.6 × 103 , whereas the local component smoothness is only Mloc ≈ 23. We keep the batch size b = 8, the a-priori baseline protocol, and the two Ada-BPSG safeguards of Section 4.2 unchanged, and again obtain F ⋆ from an independent full-batch mirror-descent solver. Figure 5 reproduces the synthetic behavior on real data. Held to the steps their own analysis follows the worst-case M , the fixed baselines stall between 10−2 and 10−3 (BPG 8.8 × 10−3 , SGD 3.6 × 10−3 , SAGA 7.1 × 10−3 , SARAH 1.4 × 10−3 ), and the “theory cap” run of Ada-BPSG sits among them at 7.1 × 10−3 : the analysis-faithful step is once more conservative on the instance. Relaxing the safeguard to the local constant, Ada-BPSG drives the gap to 2.2 × 10−6 —more than two orders of magnitude below the best a-priori baseline—tracking the O(n/K) reference slope and reaching 10−5 in a few seconds of wall-clock (Figure 5 (b)). Figure 5 (c) is transparent about the role of tuning: with its step hand-swept, fixed SAGA reaches near machine precision over a broad band of large steps, below Ada-BPSG, so on this instance a well-chosen constant step is very forgiving once found; the adaptive rule instead stays flat at ≈ 2 × 10−6 across five orders of magnitude of initial step, recovering a low gap with neither a per-instance search nor a line search, from the SAGA information it already maintains. The mechanism the synthetic experiment isolates thus governs a real unmixing problem under the identical configuration.

4.4

Sparse nonnegative matrix factorization

The third experiment tests the non-Euclidean update on sparse nonnegative matrix factorization (NMF). Given M ∈ RI1 ×I2 and a target rank r < min{I1 , I2 }, we solve 1 min Φ(U, V ) := ∥M − U V ∥2F , U,V 2 subject to U, V ⩾ 0, ∥U:j ∥0 ⩽ s1 ,

∥Vj: ∥0 ⩽ s2

(4.2) (j ∈ [r]).

I1 ×r 2 where U ∈ R+ , V ∈ Rr×I , U:j is column j of U , Vj: is row j of V , and ∥ · ∥0 counts nonzero + entries. Equivalently, the constraints form the extended-valued term h, whose indicator is zero on the feasible set and +∞ outside it. We use the quartic kernel

ψ(U, V ) = 3ψ1 (U, V ) + ∥M∥F ψ2 (U, V ), with  ψ1 (U, V ) :=

∥U ∥2F + ∥V ∥2F 2

2 ,

ψ2 (U, V ) :=

∥U ∥2F + ∥V ∥2F , 2

which matches the relative-smooth geometry of matrix factorization [32]. Prior results permit the relative-smoothness scale L̄ ⩾ 1 for this kernel [8, 43]; we set L̄ = 1. Because the cardinality constraints make h nonconvex, this experiment evaluates the algorithmic mechanism beyond the convex-regularizer setting of Assumption 3.4. On the ORL dataset of 400 normalized face images, we set r = 25, sample 5% of the components per iteration, and test four sparsity levels. We compare Ada-BPSG with BPG and with BPSG/BPSGE using SGD, SAGA, or SARAH estimators. BPSGE is the extrapolated linesearch variant. The smoothness scale M is estimated by the power method used in Wang et al. [43], and the initial step size is selected from the same grid as in the classification experiments. Each method runs for 200 epochs. 16

10 4 10 5

BPSG-SGD BPSG-SAGA BPSG-SARAH BPSGE Ada-BPSG (theory cap) Ada-BPSG (ours)

10 3 F(xK) F

10 3 F(xK) F

10 2

BPG BPSG-SGD BPSG-SAGA BPSG-SARAH BPSGE Ada-BPSG (theory cap) Ada-BPSG (ours) O(n/K) ref.

20

40 60 effective passes

80

(a) objective gap vs. effective passes

BPSG-SAGA BPSG-SGD Ada-BPSG (ours)

10 4 10 6

10 4 10 5

0

10 2

final F F

10 2

0

2

4 6 wall-clock time (s)

(b) objective gap vs. wall-clock time

10 8 10 10 10 12 10 4

10 3

10 2 10 1 initial step size

100

101

(c) sensitivity to the initial step

Figure 5: Region-pooled Poisson unmixing on the Samson hyperspectral scene, read along the same axes as Figure 4. Fixed-step baselines use their a-priori guaranteed steps; “Ada-BPSG (theory cap)” is confined to the worst-case safeguard, “Ada-BPSG (ours)” relaxes it to the local component constant. (c) sweeps the initial step—the step itself for the fixed methods, only the starting point for Ada-BPSG. Mean and standard deviation over five runs at b = 8.

(a) s1 = I1 /2 and s2 = (b) s1 = I1 /4 and s2 = (c) s1 = I1 /6 and s2 = (d) s1 = I1 /8 and s2 = I2 /2 I2 /4 I2 /6 I2 /8 Figure 6: Sparse NMF on ORL under four cardinality levels. Ada-BPSG attains the lowest objective trajectory in each regime, including against the extrapolated BPSGE methods that use line search.

Figure 6 completes the empirical chain. Ada-BPSG descends faster than the stochastic BPSG variants and remains slightly below the two extrapolated BPSGE curves throughout the four sparsity regimes. Together with the Euclidean results, this shows that the stabilized adaptive step transfers from ordinary proximal SAGA to a problem-adapted Bregman geometry without introducing line-search evaluations.

5

Conclusion

Ada-BPSG turns the gradient information already maintained by SAGA into a stable, curvatureadaptive Bregman proximal step. Its mediant aggregation downweights unreliable local BB ratios, while clipping and monotonicity supply the deterministic invariant needed for convergence. This separation keeps the method line-search-free and makes the algorithmic and theoretical roles of the step-size components explicit. The analysis follows a single chain across regimes: relative smoothness yields Bregman descent, local component smoothness controls the SAGA table error, and the safeguard balances the two. The resulting guarantees comprise an O(n/K) convex rate, a restarted linear rate under relative quadratic growth, and an O(1/K) bound on the expected squared proximal residual in the nonconvex setting. The experiments reflect the same mechanism. Ada-BPSG improves gradient efficiency and initial-step robustness in Euclidean logistic regression, exhibits the predicted convex rate on a simplex-constrained Poisson inverse problem under the entropy kernel—a non-Euclidean instance meeting every hypothesis of the analysis, where it adapts past the conservative worst-case safeguard to outrun every fixed-step baseline held to its guaranteed step, with the same advantage confirmed on real hyperspectral unmixing data—and achieves the lowest objective trajectories among the compared Bregman methods for sparse NMF, without the line search used by the extrapolated baselines.

17

Limitations and directions. The current theory analyzes the clipped sequence on a deterministic bounded region with a strongly convex kernel; the admissible safeguard depends on L̄ and the local component constant M . General finite-sum models also retain the O(nd) SAGA table, although linear-prediction losses admit the O(n) implementation in Appendix A.2. Extending the same stabilized candidate to memory-light SVRG, SARAH, or SPIDER estimators, and deriving safeguards under weaker kernel and localization conditions, are natural next steps.

18

References [1] Kimon Antonakopoulos and Panayotis Mertikopoulos. Adaptive first-order methods revisited: Convex optimization without lipschitz requirements. In Advances in Neural Information Processing Systems, volume 34, 2021. [2] Kimon Antonakopoulos, E. Veronica Belmega, and Panayotis Mertikopoulos. An adaptive mirror-prox method for variational inequalities with singular operators. In Advances in Neural Information Processing Systems, volume 32, 2019. [3] Kimon Antonakopoulos, Dong Quan Vu, Volkan Cevher, Kfir Y. Levy, and Panayotis Mertikopoulos. A universal black-box optimization method with almost dimension-free convergence rate guarantees. In Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 694–726, 2022. [4] Francis Bach and Kfir Y. Levy. A universal algorithm for variational inequalities adaptive to smoothness and noise. In Proceedings of the Thirty-Second Conference on Learning Theory, volume 99 of Proceedings of Machine Learning Research, pages 164–194. PMLR, 2019. [5] Jonathan Barzilai and Jonathan M Borwein. Two-point step size gradient methods. IMA journal of numerical analysis, 8(1):141–148, 1988. [6] Heinz H. Bauschke and Jonathan M. Borwein. Legendre functions and the method of random bregman projections. Journal of Convex Analysis, 4(1):27–67, 1997. [7] Heinz H Bauschke, Jérôme Bolte, and Marc Teboulle. A descent lemma beyond lipschitz gradient continuity: first-order methods revisited and applications. Mathematics of Operations Research, 42(2):330–348, 2017. [8] Jérôme Bolte, Shoham Sabach, Marc Teboulle, and Yakov Vaisbourd. First order methods beyond convexity and lipschitz gradient continuity with applications to quadratic inverse problems. SIAM Journal on Optimization, 28(3):2131–2151, 2018. [9] Minh N. Bui and Patrick L. Combettes. Bregman forward-backward operator splitting. Set-Valued and Variational Analysis, 29(3):583–603, 2021. [10] Oleg Burdakov, Yu-Hong Dai, and Na Huang. Stabilized barzilai-borwein method. arXiv preprint arXiv:1907.06409, 2019. [11] Dan Butnariu and Alfredo N. Iusem. Totally Convex Functions for Fixed Points Computation and Infinite Dimensional Optimization. Kluwer Academic Publishers, 2000. [12] Chih-Chung Chang and Chih-Jen Lin. Libsvm: a library for support vector machines. ACM transactions on intelligent systems and technology (TIST), 2(3):1–27, 2011. [13] Michael B. Cohen, Aaron Sidford, and Kevin Tian. Relative lipschitzness in extragradient methods and a direct recipe for acceleration. In 12th Innovations in Theoretical Computer Science Conference, volume 185 of LIPIcs, pages 62:1–62:18, 2021. doi: 10.4230/LIPIcs. ITCS.2021.62. [14] Pierre Comon, Gene Golub, Lek-Heng Lim, and Bernard Mourrain. Symmetric tensors and symmetric tensor rank. SIAM Journal on Matrix Analysis and Applications, 30(3): 1254–1279, 2008. [15] Aaron Defazio, Francis R. Bach, and Simon Lacoste-Julien. SAGA: A fast incremental gradient method with support for non-strongly convex composite objectives. In Advances in Neural Information Processing Systems, pages 1646–1654, 2014. 19

[16] Benjamin Dubois-Taine, Sharan Vaswani, Reza Babanezhad, Mark Schmidt, and Simon Lacoste-Julien. SVRG meets adagrad: Painless variance reduction. CoRR, abs/2102.09645, 2021. [17] John C. Duchi, Elad Hazan, and Yoram Singer. Adaptive subgradient methods for online learning and stochastic optimization. J. Mach. Learn. Res., 12:2121–2159, 2011. URL http://dl.acm.org/citation.cfm?id=2021068. [18] Roger Fletcher. On the barzilai-borwein method. In Optimization and control with applications, pages 235–256. Springer, 2005. [19] Filip Hanzely, Peter Richtárik, and Lin Xiao. Accelerated bregman proximal gradient methods for relatively smooth convex optimization. Computational Optimization and Applications, 79(2):405–440, 2021. doi: 10.1007/s10589-021-00273-8. [20] Marzieh Hasannasab, Johannes Hertrich, Sebastian Neumayer, Gerlind Plonka, Simon Setzer, and Gabriele Steidl. Parseval proximal neural networks. Journal of Fourier Analysis and Applications, 26:1–31, 2020. [21] Feihu Huang and Yuning Luo. Faster double adaptive gradient methods. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 39, pages 27018–27026, 2025. [22] Feihu Huang, Xidong Wu, and Heng Huang. Efficient mirror descent ascent methods for nonsmooth minimax problems. Advances in Neural Information Processing Systems, 34: 10431–10443, 2021. [23] Rie Johnson and Tong Zhang. Accelerating stochastic gradient descent using predictive variance reduction. In Advances in Neural Information Processing Systems, pages 315–323, 2013. [24] Ali Kavis, Kfir Y. Levy, Francis Bach, and Volkan Cevher. UniXGrad: A universal, adaptive algorithm with optimal guarantees for constrained optimization. In Advances in Neural Information Processing Systems, volume 32, 2019. [25] Ali Kavis, Stratis Skoulakis, Kimon Antonakopoulos, Leello Tadesse Dadi, and Volkan Cevher. Adaptive stochastic variance reduction for non-convex finite-sum minimization. Advances in Neural Information Processing Systems, 35:23524–23538, 2022. [26] Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In Yoshua Bengio and Yann LeCun, editors, 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015. URL http://arxiv.org/abs/1412.6980. [27] Sen Kuang, Benwei Wu, and Shuang Cong. A fusion measurement approach to improve quantum state tomography efficiency and accuracy. IEEE Transactions on Instrumentation and Measurement, 69(6):3049–3060, 2019. [28] Kfir Y. Levy, Alp Yurtsever, and Volkan Cevher. Online adaptive methods, universality and acceleration. In Advances in Neural Information Processing Systems, volume 31, 2018. [29] Haihao Lu, Robert M. Freund, and Yurii Nesterov. Relatively smooth convex optimization by first-order methods, and applications. SIAM Journal on Optimization, 28(1):333–354, 2018. doi: 10.1137/16M1099546. [30] Ke Ma, Jinshan Zeng, Jiechao Xiong, Qianqian Xu, Xiaochun Cao, Wei Liu, and Yuan Yao. Stochastic non-convex ordinal embedding with stabilized barzilai-borwein step size. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 32, 2018. 20

[31] Yura Malitsky and Konstantin Mishchenko. Adaptive proximal gradient method for convex optimization. Advances in Neural Information Processing Systems, 37:100670–100697, 2024. [32] Mahesh Chandra Mukkamala and Peter Ochs. Beyond alternating updates for matrix factorization with inertial bregman proximal gradient algorithms. Advances in Neural Information Processing Systems, 32, 2019. [33] Lam M. Nguyen, Jie Liu, Katya Scheinberg, and Martin Takác. SARAH: A novel method for machine learning problems using stochastic recursive gradient. In Proceedings of the 34th International Conference on Machine Learning, 2017. [34] Van Dung Nguyen and Bàng Công Vũ. A stochastic variance reduction algorithm with bregman distances for structured composite problems. Optimization, 72(6):1463–1484, 2023. doi: 10.1080/02331934.2022.2027413. [35] Boris T Polyak. Some methods of speeding up the convergence of iteration methods. Ussr computational mathematics and mathematical physics, 4(5):1–17, 1964. [36] Nicolas Le Roux, Mark Schmidt, and Francis R. Bach. A stochastic gradient method with an exponential convergence rate for finite training sets. In Advances in Neural Information Processing Systems, pages 2672–2680, 2012. [37] Zheng Shi, Nicolas Loizou, Peter Richtárik, and Martin Takác. AI-SARAH: adaptive and implicit stochastic recursive gradient methods. CoRR, abs/2102.09700, 2021. [38] Antonio Silveti-Falls, Cesare Molinari, and Jalal Fadili. A stochastic bregman primal-dual splitting algorithm for composite optimization. arXiv preprint arXiv:2112.11928, 2021. [39] Ben Somers, Maciel Zortea, Antonio Plaza, and Gregory P. Asner. Automated extraction of image-based endmember bundles for improved spectral unmixing. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 5(2):396–408, 2012. doi: 10.1109/JSTARS.2011.2181340. [40] Conghui Tan, Shiqian Ma, Yu-Hong Dai, and Yuqiu Qian. Barzilai-borwein step size for stochastic gradient descent. In Advances in Neural Information Processing Systems, pages 685–693, 2016. [41] Quang Van Nguyen. Forward-backward splitting with bregman distances. Vietnam Journal of Mathematics, 45(3):519–539, 2017. [42] Qingsong Wang and Deren Han. A bregman stochastic method for nonconvex nonsmooth problem beyond global lipschitz gradient continuity. Optimization Methods and Software, 38(5):914–946, 2023. [43] Qingsong Wang, Zehui Liu, Chunfeng Cui, and Deren Han. A bregman proximal stochastic gradient method with extrapolation for nonconvex nonsmooth problems. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 38, pages 15580–15588, 2024. [44] Qingsong Wang, Shengze Xu, Xiaojiao Tong, and Tieyong Zeng. An improved optimal proximal gradient algorithm for non-blind image deblurring. arXiv preprint arXiv:2502.07602, 2025. [45] Rachel Ward, Xiaoxia Wu, and Léon Bottou. Adagrad stepsizes: sharp convergence over nonconvex landscapes. In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, Proceedings of the 36th International Conference on Machine Learning, ICML 2019, 9-15 June 2019, Long Beach, California, USA, volume 97 of Proceedings of Machine Learning Research, pages 6677–6686. PMLR, 2019. URL http://proceedings.mlr.press/v97/ ward19a.html. 21

[46] Zhuang Yang. Improved powered stochastic optimization algorithms for large-scale machine learning. Journal of Machine Learning Research, 24(241):1–29, 2023. [47] Danqing Zhou, Shiqian Ma, and Junfeng Yang. Adabb: Adaptive barzilai-borwein method for convex optimization. arXiv preprint arXiv:2401.08024, 2024. [48] Kaiwen Zhou, Qinghua Ding, Fanhua Shang, James Cheng, Danli Li, and Zhi-Quan Luo. Direct acceleration of saga using sampled negative momentum. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 1602–1610, 2019. [49] Feiyun Zhu. Hyperspectral unmixing: Ground truth labeling, datasets, benchmark performances and survey. arXiv preprint arXiv:1708.05125, 2017.

22

A

Proofs and Supplementary Details

A.1

Auxiliary inequalities

Lemma A.1 (Three-point proximal inequality). Let   1 + x = arg min h(x) + ⟨g, x⟩ + Dψ (x, xk ) . x η Then, for any u ∈ E, h(x+ ) − h(u) + ⟨g, x+ − u⟩ ⩽

 1 Dψ (u, xk ) − Dψ (u, x+ ) − Dψ (x+ , xk ) . η

Proof. The optimality condition gives 1 0 ∈ ∂h(x+ ) + g + (∇ψ(x+ ) − ∇ψ(xk )). η Choose ξ + ∈ ∂h(x+ ) satisfying this inclusion. By convexity of h, h(u) ⩾ h(x+ ) + ⟨ξ + , u − x+ ⟩. Substituting ξ + = −g − η −1 (∇ψ(x+ ) − ∇ψ(xk )) and using the Bregman three-point identity gives the result.

A.2

Scalar-table implementation

For a linear-prediction loss fi (x) = Ψi (⟨ai , x⟩), the scalar Φki := ⟨ai , φki ⟩ determines the stored component gradient. Algorithm 2 maintains these scalars together with the aggregate gradient µk , reducing SAGA table storage from O(nd) to O(n) [48]. The two secant accumulators remain scalar, so the adaptive step requires no additional vector table.

A.3

Proof of Lemma 3.5

Let

n

aki := ∇fi (xk ) − ∇fi (φki ),

āk :=

1X k ai . n i=1

1 P

Since µk = n

k i ∇fi (φi ), the estimator can be written as

ek = 1 ∇ b

X

aki + µk .

i∈Bk

The uniform sampling assumption gives e k ] = āk + µk = ∇f (xk ). Ek [∇ Moreover, e k − ∇f (xk ) = 1 ∇ b

X

(aki − āk ).

i∈Bk

By the sampling second-moment bound in Assumption 3.4, e k − ∇f (xk )∥2∗ ⩽ σb Ek ∥∇ bn ⩽ 23

n X

∥aki − āk ∥2∗

i=1 n 4σb X

bn

i=1

∥aki ∥2∗ .

Algorithm 2 Efficient Implementation of Ada-BPSG 1: Input: max iterations K, batch-size b ∈ {1, . . . , n}, initial point x0

initial step-size η0 ∈ [ηmin , ηmax ], update frequency m ⩾ 1 BB scaling α > 0, BB buffer θ > 0, safeguard bounds 0 < ηmin ⩽ ηmax 2: Initialize: P scalar table Φ0i = ⟨ai , x0 ⟩ for i = 1, . . . , n n 1 µ0 = n i=1 ∇Ψi (Φ0i )ai , ξ0 = 0, and δ0 = 0 3: for k = 0, 1, . . . , K − 1 do 4: if k > 0 then 5: ηk = ηk−1 6: end if 7: if k > 0  and k mod m = 0 then ξk 1 , ξk > 0, α δ + θξk 8: ηbk = k  ηk−1 , ξk = 0, 9: ηk = max{ηk−1 , clip[ηmin ,ηmax ] (b ηk )} 10: ξk = 0, δk = 0 11: end if 12: Randomly pick a batch Bk uniformly without replacement. 13: Store Φk+1 = ⟨ai , xk ⟩ for i ∈ Bk i k 14: Set Φk+1 / Bk i P= Φi for i ∈
1 e ) − ∇Ψi (Φki ) ai + µk 15: ∇k = b i∈Bk ∇Ψi (Φk+1 i n o e k , x⟩ + 1 Dψ (x, xk ) 16: xk+1 = arg minx h(x) + ⟨∇ ηk
P k+1 k 17: ξk+1 = ξk + i∈Bk (Φk+1 − Φ ) ∇Ψ ) − ∇Ψi (Φki ) i (Φi i i  
2 P 18: δk+1 = δk + i∈Bk ∇Ψi (Φk+1 ) − ∇Ψi (Φki ) · ∥ai ∥2∗ i
P 19: µk+1 = µk + n1 i∈Bk ∇Ψi (Φk+1 ) − ∇Ψi (Φki ) ai i 20: end for 21: return xK

For the convex setting, the self-bounding condition implies ∥aki ∥2∗ ⩽ 2∥∇fi (xk ) − ∇fi (x∗ )∥2∗ + 2∥∇fi (φki ) − ∇fi (x∗ )∥2∗   ⩽ 4M ∆i (xk ) + ∆i (φki ) . Combining the two displays yields " # n X 1 16σ M b e k − ∇f (xk )∥2∗ ⩽ ∆(xk ) + ∆i (φki ) . Ek ∥∇ b n i=1

A.4

Proof of Theorem 3.6

Let

n

Qi (x) := ∆i (x) + H(x),

Tk :=

11X Qi (φki ), ωn

ω :=

i=1

Since φ0i = x0 and ∇f (x∗ ) + ρ∗ = 0, T0 =

n (F (x0 ) − F (x∗ )) . 4b

The table update gives   b b Ek [Tk+1 ] = 1 − Tk + (F (xk ) − F (x∗ )) . n nω 24

4b . n

e k , we first Using the upper smooth-adaptability of f and Lemma A.1 with u = x∗ and g = ∇ obtain F (xk+1 ) − F (x∗ ) ⩽ f (xk ) − f (x∗ ) − ⟨∇f (xk ), xk − x∗ ⟩ e k , xk+1 − xk ⟩ + ⟨∇f (xk ) − ∇ e k , xk − x∗ ⟩ + ⟨∇f (xk ) − ∇   1 1 Dψ (xk+1 , xk ) + + L̄ − (Dψ (x∗ , xk ) − Dψ (x∗ , xk+1 )) . ηk ηk By convexity of f , the first line is nonpositive. Taking conditional expectation and using e k ] = ∇f (xk ) gives Ek [∇ e k , xk+1 − xk ⟩ Ek [F (xk+1 ) − F (x∗ )] ⩽ Ek ⟨∇f (xk ) − ∇   1 1 Ek Dψ (xk+1 , xk ) + (Dψ (x∗ , xk ) − Ek Dψ (x∗ , xk+1 )) . + L̄ − ηk ηk By Young’s inequality and Dψ (xk+1 , xk ) ⩾ 12 ∥xk+1 − xk ∥2 , e k , xk+1 − xk ⟩ ⩽ Ek ⟨∇f (xk ) − ∇

1 e k − ∇f (xk )∥2 + βEk Dψ (xk+1 , xk ). Ek ∥∇ ∗ 2β

Take β = 32σb M/b. Since ηmax ⩽ (L̄ + β)−1 , the coefficient of Dψ (xk+1 , xk ) is nonpositive. Lemma 3.5 gives n 1 1 11X 2 e Ek ∥∇k − ∇f (xk )∥∗ ⩽ ∆(xk ) + ∆i (φki ). 2β 4 4n i=1 1 P ∗ Because H ⩾ 0, we have ∆(xk ) ⩽ F (xk ) − F (x ) and n i ∆i (φki ) ⩽ ωTk . Combining these bounds with the recursion for Tk+1 , and using b/(nω) = 1/4 and ω/4 = b/n, gives 1 Ek [F (xk+1 ) − F (x∗ ) + Tk+1 ] ⩽ Tk + (F (xk ) − F (x∗ )) 2 1 + (Dψ (x∗ , xk ) − Ek Dψ (x∗ , xk+1 )) . ηk Rearranging, 1 1 Ek [F (xk+1 ) − F (x∗ )] ⩽ Tk − Ek Tk+1 + (F (xk ) − F (x∗ ) − Ek [F (xk+1 ) − F (x∗ )]) 2 2 1 ∗ (Dψ (x , xk ) − Ek Dψ (x∗ , xk+1 )) . + ηk Let ak = 1/ηk . Since the convex theorem assumes ηk is nondecreasing, ak is nonincreasing. Therefore, K−1 X

E [ak (Dψ (x∗ , xk ) − Dψ (x∗ , xk+1 ))] ⩽ a0 Dψ (x∗ , x0 ) ⩽

k=0

1 ηmin

Dψ (x∗ , x0 ).

Summing the preceding one-step inequality and using TK ⩾ 0 yields K−1

1 X 1 1 E[F (xk+1 ) − F (x∗ )] ⩽ T0 + (F (x0 ) − F (x∗ )) + Dψ (x∗ , x0 ). 2 2 ηmin k=0

Jensen’s inequality for x̄K = K1

PK

k=1 xk completes the proof.

25

A.5

Proof of Theorem 3.8

Condition on the history up to the beginning of restart round t. The restart round reinitializes the SAGA table by setting all table points to x(t) and uses a safeguarded step-size sequence satisfying the assumptions of Theorem 3.6 within that round. Therefore the initial Lyapunov term in Theorem 3.6 has the same form with x0 replaced by x(t) . Applying Theorem 3.6 to that round and using the relative quadratic growth condition gives  1 Dψ (x∗ , x(t) ) ⩽ F (x(t) ) − F (x∗ ) . µ n + 12 , Hence, with Ab = 4b    h i 1 2 (t+1) ∗ (t) Ab + F (x(t) ) − F (x∗ ) . E F (x ) − F (x ) | x ⩽ K µηmin

Taking total expectation, the stated lower bound on K makes the contraction factor at most 1/2, and induction gives the claim.

A.6

Proof of Lemma 3.10

Using the notation from the proof of Lemma 3.5, e k − ∇f (xk ) = 1 ∇ b

X

(aki − āk ).

i∈Bk

The sampling second-moment bound and Jensen’s inequality give n

e k − ∇f (xk )∥2∗ ⩽ Ek ∥∇

4σb X k 2 ∥ai ∥∗ . bn i=1

By component Lipschitz continuity on the deterministic set D, ∥aki ∥2∗ = ∥∇fi (xk ) − ∇fi (φki )∥2∗ ⩽ M 2 ∥xk − φki ∥2 . This proves the lemma.

A.7

Proof of Theorem 3.11

Let x ek+1 = Proxψ h,ηk (xk , ∇f (xk )). e k . Using upper smooth-adaptability of f and Lemma A.1 for the stochastic Let ek := ∇f (xk ) − ∇ prox step with u = x ek+1 gives e k , xk+1 − x F (xk+1 ) ⩽ f (xk ) + h(e xk+1 ) + ⟨∇f (xk ), xk+1 − xk ⟩ − ⟨∇ ek+1 ⟩   1 1 1 + L̄ − Dψ (xk+1 , xk ) + Dψ (e xk+1 , xk ) − Dψ (e xk+1 , xk+1 ). ηk ηk ηk Applying Lemma A.1 to the full-gradient prox step with u = xk gives h(e xk+1 ) − h(xk ) + ⟨∇f (xk ), x ek+1 − xk ⟩ ⩽ −

1 1 Dψ (xk , x ek+1 ) − Dψ (e xk+1 , xk ). ηk ηk

Combining the last two displays yields   1 F (xk+1 ) ⩽ F (xk ) + ⟨ek , xk+1 − x ek+1 ⟩ + L̄ − Dψ (xk+1 , xk ) ηk 1 1 − Dψ (xk , x ek+1 ) − Dψ (e xk+1 , xk+1 ). ηk ηk 26

Young’s inequality and the strong convexity of ψ imply ⟨ek , xk+1 − x ek+1 ⟩ ⩽

1 ηk ∥ek ∥2∗ + Dψ (e xk+1 , xk+1 ). 2 ηk

Taking conditional expectation and applying Lemma 3.10, we obtain   1 1 Ek [F (xk+1 )] ⩽ F (xk ) − Dψ (xk , x Ek Dψ (xk+1 , xk ) ek+1 ) + L̄ − ηk ηk n 2σb ηk M 2 X + ∥xk − φki ∥2 . bn i=1

Define

n

1X Sk := ∥xk − φki ∥2 , n

p :=

i=1

b , n

and use the deterministic Lyapunov coefficient c :=

4σb ηmax M 2 n . b2

Let dk := xk+1 − xk and rik := xk − φki . For a fixed sampled batch Bk , the table update gives ( dk , i ∈ Bk , k+1 xk+1 − φi = k dk + ri , i ∈ / Bk . For any γ > 0, Young’s inequality yields   1 k 2 ∥dk + ri ∥ ⩽ 1 + ∥dk ∥2 + (1 + γ)∥rik ∥2 . γ Consequently, 1 X 1 X ∥dk ∥2 + ∥dk + rik ∥2 n n i∈Bk i∈B / k   1−p 1+γ X k 2 ∥dk ∥2 + ⩽ 1+ ∥ri ∥ . γ n

Sk+1 =

i∈B / k

Taking conditional expectation and using Dψ (xk+1 , xk ) ⩾ 12 ∥dk ∥2 and Ek [1{i∈B / k } ] = 1 − p, we obtain   1−p Ek [Sk+1 ] ⩽ 2 1 + Ek Dψ (xk+1 , xk ) + (1 + γ)(1 − p)Sk . (A.1) γ Combining (A.1) with the one-step bound for F (xk+1 ) gives 1 ek+1 ) Ek [F (xk+1 ) + cSk+1 ] ⩽ F (xk ) − Dψ (xk , x ηk    1 1−p + L̄ − + 2c 1 + Ek Dψ (xk+1 , xk ) ηk γ   2σb ηk M 2 + c(1 + γ)(1 − p) + Sk . b Choose γ = p/2. Since (1 + γ)(1 − p) ⩽ 1 − p/2 and ηk ⩽ ηmax , the definition of c gives c(1 + γ)(1 − p) + 27

2σb ηk M 2 ⩽ c. b

Moreover,   2(2 − p) 1−p 2(2n − b) = 2 1+ = . γ p b Therefore,   8σb ηmax M 2 n(2n − b) 1−p ⩽ 2c 1 + = Γb ηmax . γ b3 The condition in Theorem 3.11 implies   1 1−p L̄ − ⩽0 + 2c 1 + ηk γ for every k, because ηk ⩽ ηmax . Hence Ek [F (xk+1 ) + cSk+1 ] ⩽ F (xk ) + cSk −

1 Dψ (xk , x ek+1 ). ηk

Taking total expectation and summing gives, using S0 = 0, SK ⩾ 0, and F (xK ) ⩾ Finf , K−1 X k=0

 E

 1 Dψ (xk , x ek+1 ) ⩽ F (x0 ) − Finf . ηk

Since Dψ (u, v) ⩾ 21 ∥u − v∥2 and ηk ⩽ ηmax , K−1

1 X E∥xk − x ek+1 ∥2 ⩽ F (x0 ) − Finf . 2ηmax k=0

For R sampled independently and uniformly from {0, . . . , K − 1} and RR = (xR − x eR+1 )/ηmax , this yields 2(F (x0 ) − Finf ) E∥RR ∥2 ⩽ . Kηmax

B

Additional experiments

B.1

Additional experiments for logistic loss

Figures 7–9 extend the initial-step experiment in Figure 3 to batch sizes 1, 16, and 64. Ada-BPSG is denoted by SAGA-BB in this Euclidean setting. Across these batch sizes, SAGA-BB preserves the broad low-gradient region observed at b = 8, whereas the nonadaptive baselines remain sensitive to the initial scale. This consistency supports the intended role of the adaptive refresh: it reduces dependence on both the starting step and the number of component updates aggregated at each iteration.

B.2

Additional experiments for Huber loss

We next replace logistic loss by Huber loss in the same ℓ2 -regularized binary-classification problems. The Huber penalty is quadratic for residuals near zero and linear beyond a fixed threshold, combining local smoothness with reduced sensitivity to large residuals. This change tests whether the step-size behavior persists across loss geometries as well as batch sizes. Figure 10 reproduces the gradient-efficiency pattern observed with logistic loss: SAGA-BB reaches the lowest objective gap in fewer effective passes on each dataset. Figures 11–14 test the complementary stability question across four batch sizes; the consistently broad low-gradient region ties the empirical robustness to the stabilized adaptive mechanism rather than to a single logistic-regression configuration. 28

ijcnn1

10−5

10−7

10−9

10−11

10−3

10−1

10−2

101

100

10

10−8

10−4 SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−6

10−8

10−10

10−12

10−3

10−2

10−1

101

100

10−5

10−7

10−9

10−11

10−13

10−12

102

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

−2

10−3

10−6

10−10

SAGA_BB SVRG SAGA SVRG_Loopless SARAH

covtype 10−1

Gradient norm

10−4

Gradient norm

Gradient norm

10−3

w8a 100

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−2

Gradient norm

mushrooms 10−1

102

10−3

10−1

10−2

101

100

102

10−3

10−1

10−2

101

100

Initial step size

Initial step size

Initial step size

Initial step size

(a)

(b)

(c)

(d)

102

Figure 7: Initial-step sensitivity for logistic loss with batch size b = 1. Gradient norms are capped at 1 for readability. mushrooms

ijcnn1

w8a

covtype

100 10−1

10−1

10−2 10−2

10−7

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−9

10−3

10

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

−6

10−8

10

Gradient norm

10−5

Gradient norm

Gradient norm

Gradient norm

10−4

10−4

10−3

−4

10−6

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−10

10 10−1

10−2

101

100

10−8

−12

102

10−3

10−2

10−1

101

100

102

10−3

10−1

10−2

101

100

10−7

10−10

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−13

10−16

102

10−3

10−1

10−2

101

100

Initial step size

Initial step size

Initial step size

Initial step size

(a)

(b)

(c)

(d)

102

Figure 8: Initial-step sensitivity for logistic loss with batch size b = 16. Gradient norms are capped at 1 for readability. mushrooms

ijcnn1

w8a 100

10−2

10

10−2

10−5

10−7

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−9

10−11

10−3

−8

10−4

10−3 10−4 10−5

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−7

10−10

10−6

10−12 101

100

10−1

10−13

10−10

10−1

10−2

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−6

10

Gradient norm

10−4

Gradient norm

Gradient norm

10−3

covtype

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

−1

Gradient norm

10−1

10−16

10−7

102

10−3

10−2

10−1

101

100

102

10−3

10−1

10−2

101

100

102

10−3

10−1

10−2

101

100

Initial step size

Initial step size

Initial step size

Initial step size

(a)

(b)

(c)

(d)

102

Figure 9: Initial-step sensitivity for logistic loss with batch size b = 64. Gradient norms are capped at 1 for readability. w8a SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−4

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−6

0

20

10−4

10−6

10−8

10−10

40

60

80

100

10−12

120

100

10−2

10−2

10−4

10−6

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−8

10−10

0

5

Number of effective passes

10

15

20

25

30

35

10−12

40

0

20

Number of effective passes

(a)

covtype

100

Objective minus best

Objective minus best

Objective minus best

10−2

10−8

ijcnn1

100

10−2

Objective minus best

mushrooms 100

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−4

10−6

10−8

10−10

40

60

80

100

10−12

120

0

20

Number of effective passes

(b)

40

60

80

100

120

Number of effective passes

(c)

(d)

Figure 10: Gradient efficiency on four LibSVM datasets with Huber loss. Lines and shaded regions show the mean and standard deviation over five runs. mushrooms

covtype

10−2

10−3

10−3 10−4 10−5

10−5

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−7

10−9

10−3 SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−5

10

10−6

−7

10−9

10−7

10−1

100

101

102

10−6

10−8

10−12

10−13 10−2

10−4

10−10

10−11

10−3

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

100 10−1

Gradient norm

Gradient norm

10−2

w8a

10−1

Gradient norm

−1

Gradient norm

10

ijcnn1

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

100

10−3

10−2

10−1

100

101

102

10−3

10−2

10−1

100

101

102

10−3

10−2

10−1

100

Initial step size

Initial step size

Initial step size

Initial step size

(a)

(b)

(c)

(d)

101

102

Figure 11: Initial-step sensitivity for Huber loss with batch size b = 1. Gradient norms are capped at 1 for readability.

29

ijcnn1

100

10−2

10−4

Gradient norm

10

10−3 10−4 10−5

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−6 10−7 10−3

10−2

10−1

100

101

w8a

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

covtype

100 10−1

10−2

Gradient norm

10−2

Gradient norm

−1

10−6

10−8

10

10−3

Gradient norm

mushrooms 100

−4

10−6

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−10

10−12

102

10−8

10−3

10−2

10−1

100

101

102

10−3

10−2

10−1

100

101

10−5

10−7 SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−9

10−11

102

10−3

10−2

10−1

100

Initial step size

Initial step size

Initial step size

Initial step size

(a)

(b)

(c)

(d)

101

102

Figure 12: Initial-step sensitivity for Huber loss with batch size b = 8. Gradient norms are capped at 1 for readability.

10

ijcnn1

100

10−5 SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−7

10−3

10−2

10−1

100

101

covtype

100 10−1

10−2

Gradient norm

Gradient norm

Gradient norm

10−4 10−3

w8a

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−2

10−6

10−8

10−3

Gradient norm

mushrooms

−1

10−4

10−6

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−10 10−8

10−12

102

10−3

10−2

10−1

100

101

102

10−3

10−2

10−1

100

101

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−5

10−7

10−9

10−11

102

10−3

10−2

10−1

100

Initial step size

Initial step size

Initial step size

Initial step size

(a)

(b)

(c)

(d)

101

102

Figure 13: Initial-step sensitivity for Huber loss with batch size b = 16. Gradient norms are capped at 1 for readability. mushrooms

ijcnn1

100

w8a

covtype

100

100 10

10

−2

−1

10−1

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−5

10−6 10−3

10−2

10−1

100

101

102

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−6

10−8

10

−3

10−4

10−10

10−5

10−12

10−6

10−3

10−2

10−1

100

101

102

Gradient norm

10−4

Gradient norm

Gradient norm

Gradient norm

10−3

10−2

10−2

10−4 10−2

SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH 10−3

10−2

10−1

100

101

102

10−4

10−6 SAGA_BB SVRG_BB SVRG SAGA SVRG_Loopless SARAH

10−8

10−3

10−2

10−1

100

Initial step size

Initial step size

Initial step size

Initial step size

(a)

(b)

(c)

(d)

101

102

Figure 14: Initial-step sensitivity for Huber loss with batch size b = 64. Gradient norms are capped at 1 for readability.

30

Related documents

• Single-Loop Gradient Algorithms for Pessimistic Bilevel Optimization Problems#998732
• Federated stochastic bilevel optimization with fully first-order gradients#1110713
• S$^3$LDBO: A Snapshot Single-Loop Algorithm for Decentralized Bilevel Optimization#246538
• Convergence Rate Analysis for Monotone Accelerated Proximal Gradient Method#629792
• Convergence Rate Analysis for Monotone Accelerated Proximal Gradient Method#627930