Conceptio › Archive › arXiv CS

arXiv CSopen access

Randomized Subspace Nesterov Accelerated Gradient

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

neural-networks

machine learning, deep learning, neural networks

Randomized Subspace Nesterov Accelerated Gradient Gaku Omiya1,2 , Pierre-Louis Poirion2 , and Akiko Takeda1,2

arXiv:2605.00740v1 [math.OC] 1 May 2026

Department of Mathematical Informatics, The University of Tokyo, Tokyo, Japan 2 Center for Advanced Intelligence Project, RIKEN, Tokyo, Japan

Abstract Randomized-subspace methods reduce the cost of first-order optimization by using only low-dimensional projected-gradient information, a feature that is attractive in forward-mode automatic differentiation and communication-limited settings. While Nesterov acceleration is well understood for full-gradient and coordinate-based methods, obtaining accelerated methods for general subspace sketches that use only projected-gradient information and can improve over full-dimensional Nesterov acceleration in oracle complexity is technically nontrivial. We develop randomized-subspace Nesterov accelerated gradient methods for smooth convex and smooth strongly convex optimization under matrix smoothness and generic sketch moment assumptions. The key technical ingredient is a three-sequence formulation tailored to matrix smoothness, which recovers the corresponding classical Nesterov methods in the full-dimensional case. The resulting theory establishes accelerated oracle-complexity guarantees and makes explicit how matrix smoothness and the sketch distribution enter the complexity. It also provides a unified basis for comparing sketch families and identifying when randomized-subspace acceleration improves over full-dimensional Nesterov acceleration in oracle complexity.

Keywords: Randomized subspace methods; Nesterov acceleration; Convergence analysis; Convex optimization; Matrix smoothness

Introduction

We consider the unconstrained optimization problem: min f (x),

x∈Rd

where f : Rd → R is differentiable, and focus on both the smooth convex and smooth strongly convex settings. Whenever a minimizer exists, we denote by x⋆ an optimal solution and set f ⋆ := f (x⋆ ). First-order methods are the workhorse of large-scale optimization due to their scalability and low per-iteration cost. Among them, acceleration techniques—most notably Nesterov’s accelerated gradient (NAG) method [24, 25]—play a fundamental role. They achievepoptimal convergence rates, improving from O(1/N ) to O(1/N 2 ) in the convex setting, and from L/µ to L/µ dependence in the strongly convex case. As a result, acceleration has become an indispensable component in modern optimization algorithms. High-dimensional optimization and randomized subspace methods. Many machine learning applications give rise to large-scale, high-dimensional optimization problems. To address this challenge, randomized coordinate and block-coordinate descent methods have been extensively studied [22, 26, 33, 39]. At each iteration, these methods update only one coordinate or a small subset of coordinates, thereby reducing the per-iteration cost. A natural generalization is given by randomized subspace methods, which update along a randomly chosen low-dimensional subspace. A representative update, introduced by [17], is xk+1 = xk − αk Pk Pk⊤ ∇f (xk ), 1

Table 1: Comparison of oracle complexity under a directional-derivative oracle model, where each directionalderivative evaluation counts as one oracle call. Here R0 := ∥x0 − x⋆ ∥, ∆0 := f (x0 ) − f (x⋆ ), and the sketch parameters ℓ and ω are defined in Assumption 2.2. The RS-GD bounds of [17] are re-evaluated under the same assumptions as our proposed method; see Appendix H. After instantiating ω and ℓ for standard sketches, the resulting oracle factors are compared in Section 5. Method

Convex

Strongly convex ∆0 GD O dR0 ϵ O dL log ϵ µ ∆0 L 2L Prop. H.1 Prop. H.2 RS-GD [17] O ωrR0 ϵ O ℓr µ log ϵ q q ∆0 NAG [24, 25] O dR0 Lϵ O d L µ log ϵ q √ √ q ∆0 RS-NAG (this work) O ωℓr2 R0 Lϵ Thm. 3.3 O ωℓr2 L Thm. 4.3 µ log ϵ 2L

where Pk ∈ Rd×r is a random sketch matrix with r ≤ d; here, r corresponds to the dimension of the subspace. This framework includes randomized coordinate and block-coordinate updates as special cases, and recovers standard gradient descent when r = d and Pk Pk⊤ = Id . Randomized subspace methods are attractive for several reasons. They provide more flexibility than coordinate-aligned updates, and they are particularly effective under memory and communication constraints. For instance, reverse-mode AD is efficient for full-gradient computation but typically requires storing intermediate quantities, whereas forward-mode AD can compute directional derivatives with lower memory overhead. In a forward-mode implementation, computing the full gradient may require d directional derivatives, whereas subspace methods require only r. Here, Pk Pk⊤ ∇f (x) is understood as Pk (Pk⊤ ∇f (x)), with only Pk⊤ ∇f (x) queried. Similarly, in distributed optimization, transmitting an r-dimensional sketch is significantly cheaper than communicating a full gradient. These advantages have led to growing interest in randomized subspace gradient methods [3, 8, 17, 27–29]. Limitations of existing acceleration methods. Despite this progress, a fundamental gap remains. While acceleration is essential in first-order optimization, its integration into randomized subspace methods is still poorly understood. For coordinate descent methods, accelerated variants have been developed, including Nesterov’s accelerated coordinate descent [26] and subsequent refinements such as APCG [21], APPROX [7], ALPHA [32], and non-uniform sampling schemes [1]. However, these methods rely heavily on coordinate-wise structure and do not extend naturally to general subspace directions. Acceleration has also been studied for compressed gradient descent in distributed optimization [20]. However, under our oracle model, this does not give the desired oracle-complexity advantage over full-dimensional Nesterov acceleration. A detailed comparison with compressed-gradient methods is deferred to Appendix A. A natural attempt is to directly combine Nesterov acceleration with randomized subspace gradients by replacing the full gradient with Pk Pk⊤ ∇f : xk+1 = yk − ηPk Pk⊤ ∇f (yk ),

yk+1 = xk+1 + βk (xk+1 − xk ).

However, accelerated guarantees for this direct two-sequence scheme are not obtained by a straightforward adaptation of the classical analysis. The classical two-sequence Nesterov analysis relies on delicate estimatesequence arguments, which do not carry over directly to randomized subspace updates. In particular, the direct argument closes only under restrictive near-full-dimensional conditions; a precise discussion is deferred to Appendix J. Our approach and contributions.

These observations lead to the following question:

Can one design accelerated methods that use only randomized subspace gradients while achieving improved oracle complexity?

In this paper, we answer this question affirmatively. We propose randomized-subspace Nesterov accelerated gradient (RS-NAG) methods, to our knowledge the first accelerated framework for general randomized subspace gradient methods that can achieve favorable oracle complexity compared with both non-accelerated randomized-subspace methods and standard NAG; see Table 1. Our key technical contribution is a novel three-sequence formulation tailored to matrix smoothness, which combines a sketched descent step with an auxiliary estimate sequence and enables a clean convergence analysis. Our contributions. • Accelerated randomized subspace methods. We propose Nesterov-type randomized subspace methods for both convex and strongly convex optimization, recovering the corresponding classical Nesterov methods in the full-dimensional case. • Oracle complexity under matrix smoothness. We prove convergence and oracle-complexity bounds that capture the interaction between matrix smoothness and the sketch distribution. • Comparison with full-dimensional acceleration. Our bounds identify when randomized subspace acceleration can outperform full-dimensional Nesterov acceleration in oracle complexity. • Unified comparison of sketching strategies. We analyze Haar, coordinate, and Gaussian sketches, revealing their relative convergence bounds and identifying optimal sketch dimensions in terms of oracle complexity. Additional probability guarantees. Beyond the expectation bounds stated in Theorems 3.3 and 4.3, Appendix I also provides uniform-in-time high-probability bounds and almost-sure eventual rates with only mild losses. These results show that the accelerated behavior predicted by the expectation bounds is not merely an average-over-runs phenomenon, but persists with high probability uniformly over time and eventually along almost every run. Notation: Unless stated otherwise, ∥ · ∥ denotes the Euclidean norm for vectors and the corresponding operator norm for matrices induced by the Euclidean norm. Let Id denote the d × d identity matrix. For a matrix A ∈ Rd×d , diag(A) ∈ Rd×d denotes the diagonal matrix whose diagonal entries coincide with those of A. For symmetric matrices A, B ∈ Rd×d , we write A ⪯ B if B − A is positive semidefinite. Throughout the paper, we assume that the ambient dimension satisfies d ≥ 2. For sketch-based methods, we assume 1 ≤ r ≤ d.

Preliminaries

2.1

Oracle complexity

We measure oracle complexity in terms of directional-derivative queries. Under the oracle model considered here, which is compatible with forward-mode AD, a full-dimensional gradient evaluation costs d directionalderivative queries. In contrast, a projected gradient evaluation, which computes a projection of the gradient onto a lower-dimensional subspace, e.g., P ⊤ ∇f (x) for a matrix P ∈ Rd×r with r ≤ d, costs only r such queries. Thus, in this oracle model, the cost scales with r rather than d. This oracle measure is also relevant in distributed settings with communication bottlenecks, where transmitting an r-dimensional sketch requires sending r real numbers, as opposed to d for a full gradient, and therefore also reduces the communication volume from d to r scalars; see Section B for details. Accordingly, our complexity results are intended for forward-mode and/or communication-bottlenecked regimes, rather than as a claim about universal wall-clock speedups across all implementations.

2.2

Problem setting and sketch assumptions

Assumption 2.1 (Matrix smoothness). Let f : Rd → R be differentiable. Assume that there exists a fixed matrix L ∈ Rd×d with L ⪰ 0, L ̸= 0, such that, for all x, y ∈ Rd , 1 f (y) ≤ f (x) + ⟨∇f (x), y − x⟩ + (y − x)⊤ L(y − x). 2

(1)

Define L := ∥L∥. Since L ⪯ LId , (1) implies the standard scalar L-smooth descent lemma. Assumptions of this form have recently been used increasingly in matrix-smooth optimization and compression; see, e.g., Flynn et al. [8], Hanzely and Richtárik [13], Li et al. [19], Maranjyan et al. [23], Wang et al. [37]. Assumption 2.2 (Sketch moment conditions). Let P ∈ Rd×r be a random matrix. Assume that for some constants ω > 0 and ℓ > 0, E[P P ⊤ ] = Id ,

(2)

⊤ 2

(3)

⊤

(4)

E[(P P ) ] ⪯ ωId , ⊤

E[P P LP P ] ⪯ ℓL Id .

We will show in Section 5 that this assumption is satisfied by Haar, coordinate, and Gaussian sketches. In all algorithms and analyses below, the sketch matrices {Pk }k≥0 are assumed to be i.i.d. copies of a random matrix P ∈ Rd×r . Proposition 2.3. Under Assumptions 2.1 and 2.2, any constant ω satisfying (3) satisfies ω ≥ d/r, and any constant ℓ satisfying (4) satisfies ℓ ≥ 1. Consequently, every admissible pair (ω, ℓ) satisfies ℓω ≥ 1. Moreover, for any ω satisfying (3), (4) holds with ℓ = ω. The proof is deferred to Section C. By the last statement of Proposition 2.3, for any admissible ω, (4) also holds with ℓ = ω. Therefore, we assume without loss of generality throughout the rest of the paper that ℓ ≤ ω. The explicit constants derived in Section 5 and summarized in Table 2 satisfy this convention.

RS-NAG for convex problems (RS-NAG-C)

We first consider the convex, not necessarily strongly convex, setting. Assumption 3.1 (Convexity). Let f : Rd → R be differentiable and convex. Assume that f admits a minimizer. Under Assumption 3.1, we propose RS-NAG for convex problems (RS-NAG-C), a randomized-subspace variant of standard NAG for smooth convex optimization [24, 25]. The method is given in Algorithm 1. Compared with standard NAG, RS-NAG-C uses only the sketched gradient Pk Pk⊤ ∇f (yk ) in place of the full gradient ∇f (yk ). The corresponding two-sequence NAG recursion is recalled in Section D, (9)–(12). We prove below that RS-NAG-C reduces to this recursion when r = d and Pk Pk⊤ = Id . Proofs for this section are deferred to Section D. Proposition 3.2 (Convex case: reduction to standard Nesterov). Suppose Assumptions 2.1 and 2.2 and Assumption 3.1 hold. Consider the full-sketch case r = d, and assume that Pk Pk⊤ = Id

for all k ≥ 0,

Algorithm 1 RS-NAG-C Require: x0 ∈ Rd , constants L, ω, ℓ > 0, and an i.i.d. sketch sequence {Pk }k≥0 1: m ← 1/(2Lℓ) 2: A0 ← 0, z0 ← x0 3: for k = 0, 1, 2, . . √ . do m + m2 + 2ωmAk 4: ak+1 ← ω 5: Ak+1 ← Ak + ak+1 Ak ak+1 6: yk ← xk + zk Ak+1 Ak+1 1 7: xk+1 ← yk − Pk Pk⊤ ∇f (yk ) Lℓ 8: zk+1 ← zk − ak+1 Pk Pk⊤ ∇f (yk ) 9: end for Therefore, the proposed randomized-subspace methods in the convex setting can be viewed as generalizations of the standard NAG. We next state the convergence guarantee for RS-NAG-C. Theorem 3.3. Suppose Assumptions 2.1 and 2.2 and Assumption 3.1 hold, and let {xk , yk , zk } be generated by Algorithm 1. Then, for all N ≥ 1, 2

E[f (xN ) − f ⋆ ] ≤ 2Lωℓ

∥x0 − x⋆ ∥ . N2

(5)

In particular, for R0 := ∥x0 − x⋆ ∥, the iteration complexity to guarantee E[f (xN ) − f ⋆ ] ≤ ϵ is ! r Lωℓ , N = O R0 ϵ and since one iteration uses r oracle calls, the oracle complexity is ! r Lωℓr2 #Oracle = rN = O R0 . ϵ Compared with the convex RS-GD √ bound in Table 1, RS-NAG-C improves the accuracy dependence of the oracle complexity from 1/ϵ to 1/ ϵ, matching the acceleration effect of NAG over GD.

RS-NAG for strongly convex problems (RS-NAG-SC)

We next consider the strongly convex setting. Assumption 4.1 (Strong convexity). Let f : Rd → R be differentiable and µ-strongly convex with µ > 0. Under Assumption 4.1, we propose RS-NAG for strongly convex problems (RS-NAG-SC), a randomizedsubspace variant of standard NAG for smooth strongly convex optimization [24, 25]. The method is given in Algorithm 2. Compared with standard NAG, RS-NAG-SC uses only the sketched gradient Pk Pk⊤ ∇f (yk ) in place of the full gradient ∇f (yk ). The corresponding two-sequence NAG recursion is recalled in Section E, (39)–(41). We prove below that RS-NAG-SC reduces to this recursion when r = d and Pk Pk⊤ = Id . Proofs for this section are deferred to Section E. Proposition 4.2 (Strongly convex case: reduction to standard Nesterov). Suppose Assumptions 2.1 and 2.2 and Assumption 4.1 hold. Consider the full-sketch case r = d, and assume that Pk Pk⊤ = Id

for all k ≥ 0, 5

Algorithm 2 RS-NAG-SC d Require: px0 ∈ R , constants L, µ, ω, ℓ > 0, and an i.i.d. sketch sequence {Pk }k≥0 1: θ ← µ/(Lωℓ) 2: z0 ← x0 3: for k = 0, 1, 2, . . . do θ 1 xk + zk 4: yk ← 1+θ 1+θ 1 5: xk+1 ← yk − Pk Pk⊤ ∇f (yk ) Lℓ θ 6: zk+1 ← (1 − θ)zk + θyk − Pk Pk⊤ ∇f (yk ) µ 7: end for

which is the case, for example, for the Haar and block coordinate sketches. Then Assumption 2.2 is satisfied with ω = 1, ℓ = 1. With this choice, the (xk , yk )-sequence generated by Algorithm 2 coincides with classical two-sequence Nesterov accelerated gradient method for strongly convex objectives. Therefore, the proposed randomized-subspace methods in the strongly convex setting can be viewed as generalizations of the standard NAG. We next state the convergence guarantee for RS-NAG-SC. Theorem 4.3. Suppose Assumptions 2.1 and 2.2 and Assumption 4.1 hold, and let {xk , yk , zk } be generated by Algorithm 2. Define ∆0 := f (x0 ) − f ⋆ . Then the parameter θ in Algorithm 2 satisfies θ ∈ (0, 1], and for all N ≥ 0, E[f (xN ) − f ⋆ ] ≤ 2(1 − θ)N ∆0 .

(6)

In particular, the iteration complexity to guarantee E[f (xN ) − f ⋆ ] ≤ ϵ is s ! Lωℓ ∆0 log , N =O µ ϵ and since one iteration uses r oracle calls, the oracle complexity is s ! Lωℓr2 ∆0 #Oracle = rN = O log . µ ϵ Compared with the strongly convex RS-GD bound inp Table 1, RS-NAG-SC improves the conditionnumber dependence of the oracle complexity from L/µ to L/µ, matching the acceleration effect of NAG over GD.

Examples of sketches: Haar, coordinate, and Gaussian

We verify Assumption 2.2 for three standard sketches: Haar, Coordinate, and Gaussian sketches. We then compare the resulting sketch-dependent oracle-complexity factors. Proofs for this section are deferred to Section F. Throughout this section, let L = ∥L∥, where L is the matrix in Assumption 2.1, and define reff :=

tr(L) , L

δdiag :=

∥ diag(L)∥ . L

√ The oracle-complexity bounds in Theorems 3.3 and 4.3 are governed√ by the sketch-dependent factor ωℓr2 . For full-dimensional Nesterov, the corresponding factor is d. Hence ωℓr2 < d means an improvement over full-dimensional Nesterov in our oracle bound. 6

Table 2: Sketch-dependent constants for Haar, coordinate, and Gaussian sketches. Here β = d(d − r)/((d + 2)(d − 1)). Sketch Haar Coordinate Gaussian

√

ℓ

d r d r d+r+1 r

reff d r d r−1 d−r + δdiag r d−1 d−1 r + 1 + reff r d

1−β+β

ωℓr2 r reff d 1−β+β d r r−1 d−r d + δdiag d−1 d−1 p (d + r + 1)(r + 1 + reff )

Proposition 5.1. Let d ≥ 2, 1 ≤ r ≤ d, L ⪰ 0, and L =p∥L∥ > 0. Then Assumption 2.2 holds with the constants in Table 2 for the following sketches: Haar P =p d/r R, where R consists of the first r columns of a Haar-distributed orthogonal matrix; Coordinate P = d/r S, where S consists of r uniformly sampled distinct columns of Id ; and Gaussian Pij ∼ N (0, 1/r) i.i.d. Moreover, 1 reff ≤ δdiag ≤ 1, δdiag ≥ . d d √ We next minimize the sketch-dependent oracle factor ωℓr2 over the sketch dimension r. Here ω and ℓ are the values in Table 2, which generally depend on r. 1 ≤ reff ≤ d,

Proposition 5.2 (Optimal sketch dimension and comparison of sketch constants). √ For each of the Haar, Coordinate, and Gaussian sketches, with ω and ℓ chosen as in Table 2, the factor ωℓr2 is minimized over r ∈ {1, . . . , d} at r = 1. At r = 1, the Haar, Coordinate, and Gaussian factors, denoted by QH , QC , and QG , respectively, are as follows. Values and ranges. r r reff + 2 3 QH = d , d ≤ QH ≤ d, d+2 d+2 √ p QC = d δdiag , d ≤ QC ≤ d, p p QG = (d + 2)(reff + 2), 3(d + 2) ≤ QG ≤ d + 2. Relations between sketch factors. 2 QH , QG = 1 + d

QH ≤

√

3 QC .

Moreover, for L = e1 e⊤ 1, r QH =

3 QC . d+2

Proposition 5.2 shows that the best bound in each sketch family is attained at r = 1. Since the√ fulldimensional Nesterov factor is d, the ranges above show that the factor Q can decrease to about d in favorable cases. Haar and Coordinate are always no worse than full-dimensional Nesterov, and Haar is always better than Gaussian, although √ the two are nearly identical for large d. Compared with Coordinate, Haar is never worse by more than 3, while it can be much better, as shown by the final example.

Numerical Experiments

We evaluate the convex and strongly convex versions of RS-NAG, namely RS-NAG-C and RS-NAG-SC. Unless otherwise stated, each curve is the mean over independent random seeds, and the shaded region denotes mean ± one standard deviation across runs. The seeds determine the Gaussian initialization; for randomized-subspace methods, they also determine the sampled sketch sequence. 7

(a) Convex diagonal

(b) Convex dense

Figure 1: Oracle-axis convergence on the four quadratic problems. The horizontal axis shows the number of oracle calls, and the vertical axis shows the objective gap f (xk ) − f ⋆ on a logarithmic scale.

6.1

Quadratic objectives

We first consider four quadratic objectives f (x) = 12 x⊤ Lx on Rd , with d = 1000 and f ⋆ = 0. These instances are designed to isolate the effects of the effective rank reff and the diagonal quantity δdiag appearing in the sketch-dependent constants. The first two instances are convex but not strongly convex, while the last two are strongly convex. The diagonal instances have small reff and large δdiag , whereas the dense instances have small δdiag . The four matrices L, together with the corresponding values of L, µ, reff , and δdiag , are specified in Appendix G.1. We use oracle budget 10,000. For randomized-subspace methods, we set r = 1 in the main experiments, the theoretically preferred choice. Appendix G.2 reports an r-sweep. We use 10 independent random seeds with x0 ∼ N (0, Id ) for each seed. For randomized-subspace methods, we consider Haar, Block-coordinate, and Gaussian sketches. We compare NAG with RS-NAG-C in the convex setting and RS-NAG-SC in the strongly convex setting, and plot f (xk ) − f ⋆ vs. oracle calls. Discussion. The results are consistent with the theoretical predictions in Proposition 5.2. On the diagonal instances, where reff is small and δdiag is large, Haar and Gaussian sketches outperform both the Blockcoordinate sketch and full-dimensional Nesterov acceleration. On the dense instances, where δdiag is small, the Block-coordinate sketch becomes the most effective, while Haar and Gaussian remain competitive. Thus, the observed oracle-axis behavior reflects the sketch-dependent quantities QH , QG , QC .

6.2

Logistic regression

We next evaluate RS-NAG-SC on ℓ2 -regularized logistic regression. Given binary classification data {(ai , yi )}ni=1 , where ai ∈ Rd and yi ∈ {−1, +1}, we consider f (x) =

n µ 1X ∥x∥22 , log 1 + exp(−yi a⊤ i x) + n i=1 2

µ > 0.

(7)

Since µ > 0, the objective is µ-strongly convex. For the matrix smoothness constant, let A ∈ Rn×d be the data matrix whose i-th row is a⊤ i , and use L=

1 ⊤ A A + µId , 4n

(8)

see Appendix G.3 for the derivation. We evaluate six real-world binary-classification benchmarks: colon-cancer [2], hiva agnostic [11], bioresponse [12], gisette [10], leukemia [9], and duke (Duke breast-cancer) [38]. For each dataset, we set µ = 1/n, use r = 1, and set L = ∥L∥, computed numerically as the largest eigenvalue of the matrix L in (8). We compare GD, NAG-SC, RS-GD, and RS-NAG-SC with Haar, coordinate, and Gaussian sketches. We initialize each run from a Gaussian random vector and plot the objective gap f (xk ) − fref against oracle calls,

(a) colon-cancer

(b) hiva agnostic

(d) gisette

(e) leukemia

(f) duke

Figure 2: Oracle-axis comparison for ℓ2 -regularized logistic regression on six real-world datasets. The horizontal axis shows oracle calls, and the vertical axis shows f (xk ) − fref on a logarithmic scale, where fref is computed by L-BFGS-B. We compare GD, NAG-SC, RS-GD, and RS-NAG-SC with Haar, coordinate, and Gaussian sketches. Each curve is the mean over 3 random seeds, and the shaded region shows one standard deviation. For each dataset, µ = 1/n, r = 1, and the oracle budget is 1,000,000. where fref is computed by L-BFGS-B [4, 41]. The oracle axis can also be read as communication bits under fixed-precision distributed implementations.1 The results are shown in Figure 2, and the corresponding dataset-dependent quantities are summarized in Table 3. Additional implementation details, dataset sources, and reference-solver details are provided in Appendix G. Appendix G.4 further reports experiments on six additional real-world datasets. Discussion. Overall, RS-NAG-SC performs strongly, and the empirical trends align with the datasetdependent quantities in Table 3. For hiva agnostic, bioresponse, and leukemia, QH and QG are smaller than QC , and the Haar and Gaussian sketches indeed perform better. By contrast, for colon-cancer, gisette, and duke, QC is smaller than QH and QG , which is consistent with the relatively strong performance of the coordinate sketch. In particular, for colon-cancer and leukemia, where the corresponding Q-values are relatively close, the empirical performance of the three sketches is also broadly comparable. Overall, these results suggest that the Q values can serve as a useful practical guide when choosing the sketch distribution before running the method.

Conclusion

We introduced randomized-subspace Nesterov accelerated gradient methods for smooth convex and strongly convex optimization under an oracle model, following the same cost-sensitive viewpoint as randomizedsubspace methods in general. The methods use projected gradients, recover standard full-dimensional Nesterov acceleration when r = d, and enjoy oracle-complexity guarantees under matrix smoothness. For three 1 This follows the communication-bit accounting in the experiments of Li et al. [20], where an r-dimensional sparse message is counted as 32r bits. The essential comparison is d versus r transmitted scalars per iteration; with b-bit scalars, these correspond to bd and br bits.

Table 3: Dataset-dependent quantities for the datasets. Here d is the ambient dimension, n is the number of training samples, and QH , QG , QC denote the r = 1 constants defined in Proposition 5.2 for the Haar, Gaussian, and Coordinate sketches, respectively. Dataset colon-cancer hiva agnostic bioresponse gisette leukemia duke

reff

2000 1617 1776 5000 7129 7129

132.1725 92.3197 88.0942 129.9973 215.7402 189.2243

132.3047 92.4338 88.1934 130.0493 215.8007 189.2773

116.1338 226.1505 231.3375 86.1400 237.3443 146.8262

6.7435 3.2773 2.3746 1.3812 4.5306 3.0240

δdiag

0.0034 62 0.0196 4229 0.0170 3751 0.0003 6000 0.0011 38 0.0004 44

standard sketch distributions—Haar, coordinate, and Gaussian sketches—we derived explicit rates, compared the resulting sketch-dependent constants, and identified the theoretically preferred subspace dimension under our oracle model. Future work includes designing sketch distributions beyond the canonical choices considered here. In particular, inspired by non-uniform sampling in accelerated coordinate descent [1], it would be interesting to develop L-aware sketch distributions that lead to faster convergence in oracle calls.

Acknowledgments This project has been partially supported by the Japan Society for the Promotion of Science (JSPS) through JSPS KAKENHI Grant Number JP23H03351 and JST CREST Grant Number JPMJCR24Q2.

References [1] Zeyuan Allen-Zhu, Zheng Qu, Peter Richtárik, and Yang Yuan. Even faster accelerated coordinate descent using non-uniform sampling. In Proceedings of the 33rd International Conference on Machine Learning, volume 48 of Proceedings of Machine Learning Research, pages 1110–1119, 2016. [2] Uri Alon, Naama Barkai, Daniel A. Notterman, Kurt Gish, Suzanne Ybarra, Daniel Mack, and Arnold J. Levine. Broad patterns of gene expression revealed by clustering analysis of tumor and normal colon tissues probed by oligonucleotide arrays. Proceedings of the National Academy of Sciences, 96(12): 6745–6750, 1999. [3] Atılım Güneş Baydin, Barak A. Pearlmutter, Don Syme, Frank Wood, and Philip Torr. Gradients without backpropagation, 2022. arXiv:2202.08587. [4] Richard H. Byrd, Peihuang Lu, Jorge Nocedal, and Ciyou Zhu. A limited memory algorithm for bound constrained optimization. SIAM Journal on Scientific Computing, 16(5):1190–1208, 1995. [5] Coralia Cartis and Lindon Roberts. Scalable subspace methods for derivative-free nonlinear least-squares optimization. Mathematical Programming, 199(1–2):461–524, 2023. [6] Chih-Chung Chang and Chih-Jen Lin. LIBSVM: A library for support vector machines. ACM Transactions on Intelligent Systems and Technology, 2(3):27:1–27:27, 2011. [7] Olivier Fercoq and Peter Richtárik. Accelerated, parallel, and proximal coordinate descent. SIAM Journal on Optimization, 25(4):1997–2023, 2015. [8] Thomas Flynn, Patrick Johnstone, and Shinjae Yoo. Problem-dependent convergence bounds for randomized linear gradient compression, 2024. arXiv:2411.12898.

[9] Todd R. Golub, Donna K. Slonim, Pablo Tamayo, Christine Huard, Michelle Gaasenbeek, Jill P. Mesirov, Hilary Coller, Mignon L. Loh, James R. Downing, Mark A. Caligiuri, Clara D. Bloomfield, and Eric S. Lander. Molecular classification of cancer: Class discovery and class prediction by gene expression monitoring. Science, 286(5439):531–537, 1999. [10] Isabelle Guyon, Steve Gunn, Asa Ben-Hur, and Gideon Dror. Result analysis of the NIPS 2003 feature selection challenge. In Advances in Neural Information Processing Systems, volume 17, pages 545–552, 2005. [11] Isabelle Guyon, Amir Saffari, Gideon Dror, and Gavin C. Cawley. Agnostic learning vs. prior knowledge challenge. In Proceedings of the International Joint Conference on Neural Networks, pages 829–834, 2007. [12] Ben Hamner, dcthompson, and Jorg. Predicting a biological response. Kaggle competition, 2012. URL https://kaggle.com/competitions/bioresponse. [13] Filip Hanzely and Peter Richtárik. Accelerated coordinate descent with arbitrary sampling and best rates for minibatches. In Proceedings of the 22nd International Conference on Artificial Intelligence and Statistics, volume 89 of Proceedings of Machine Learning Research, pages 304–312, 2019. [14] Warren Hare, Lindon Roberts, and Clément W. Royer. Expected decrease for derivative-free algorithms using random subspaces. Mathematics of Computation, 94(351):277–304, 2025. [15] Mingyi Hong, Davood Hajinezhad, and Ming-Min Zhao. Prox-PDA: The proximal primal-dual algorithm for fast distributed nonconvex optimization and learning over networks. In Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pages 1529–1538, 2017. [16] Markelle Kelly, Rachel Longjohn, and Kolby Nottingham. The UCI machine learning repository. https: //archive.ics.uci.edu, 2023. [17] David Kozak, Stephen Becker, Alireza Doostan, and Luis Tenorio. A stochastic subspace approach to gradient-free optimization in high dimensions. Computational Optimization and Applications, 79(2): 339–368, 2021. [18] David Kozak, Cesare Molinari, Lorenzo Rosasco, Luis Tenorio, and Silvia Villa. Zeroth-order optimization with orthogonal random directions. Mathematical Programming, 199(1–2):1179–1219, 2023. [19] Hanmin Li, Avetik Karagulyan, and Peter Richtárik. Det-CGD: Compressed gradient descent with matrix stepsizes for non-convex optimization. In Proceedings of the 12th International Conference on Learning Representations, 2024. [20] Zhize Li, Dmitry Kovalev, Xun Qian, and Peter Richtárik. Acceleration for compressed gradient descent in distributed and federated optimization. In Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 5895–5904, 2020. [21] Qihang Lin, Zhaosong Lu, and Lin Xiao. An accelerated proximal coordinate gradient method. In Advances in Neural Information Processing Systems, volume 27, pages 3059–3067, 2014. [22] Zhaosong Lu and Lin Xiao. On the complexity analysis of randomized block-coordinate descent methods. Mathematical Programming, 152(1–2):615–642, 2015. [23] Artavazd Maranjyan, Mher Safaryan, and Peter Richtárik. Gradskip: Communication-accelerated local gradient methods with better computational complexity. Transactions on Machine Learning Research, 2025. [24] Yurii Nesterov. A method of solving a convex programming problem with convergence rate O(1/k 2 ). Soviet Mathematics Doklady, 27(2):372–376, 1983.

[25] Yurii Nesterov. Introductory Lectures on Convex Optimization: A Basic Course, volume 87 of Applied Optimization. Springer, 2004. [26] Yurii Nesterov. Efficiency of coordinate descent methods on huge-scale optimization problems. SIAM Journal on Optimization, 22(2):341–362, 2012. [27] Yurii Nesterov and Vladimir Spokoiny. Random gradient-free minimization of convex functions. Foundations of Computational Mathematics, 17(2):527–566, 2017. [28] Ryota Nozawa, Pierre-Louis Poirion, and Akiko Takeda. Randomized subspace gradient method for constrained optimization, 2023. arXiv:2307.03335. [29] Gaku Omiya, Pierre-Louis Poirion, and Akiko Takeda. Convergence analysis of randomized subspace normalized SGD under heavy-tailed noise, 2026. arXiv:2601.20399. [30] John C. Platt. Fast training of support vector machines using sequential minimal optimization. In Bernhard Schölkopf, Christopher J. C. Burges, and Alexander J. Smola, editors, Advances in Kernel Methods: Support Vector Learning, pages 185–208. MIT Press, Cambridge, MA, 1999. [31] Danil Prokhorov. Ijcnn 2001 neural network competition. Slide presentation in IJCNN’01, Ford Research Laboratory, 2001. [32] Zheng Qu and Peter Richtárik. Coordinate descent with arbitrary sampling I: Algorithms and complexity. Optimization Methods and Software, 31(5):829–857, 2016. [33] Peter Richtárik and Martin Takáč. Iteration complexity of randomized block-coordinate descent methods for minimizing a composite function. Mathematical Programming, 144(1–2):1–38, 2014. [34] Lindon Roberts and Clément W. Royer. Direct search based on probabilistic descent in reduced spaces. SIAM Journal on Optimization, 33(4):3057–3082, 2023. [35] Mher Safaryan, Filip Hanzely, and Peter Richtárik. Smoothness matrices beat smoothness constants: Better communication compression techniques for distributed optimization. In Advances in Neural Information Processing Systems, volume 34, pages 25688–25702, 2021. [36] Joaquin Vanschoren, Jan N. van Rijn, Bernd Bischl, and Luis Torgo. OpenML: Networked science in machine learning. ACM SIGKDD Explorations Newsletter, 15(2):49–60, 2013. [37] Bokun Wang, Mher Safaryan, and Peter Richtárik. Theoretically better and numerically faster distributed optimization with smoothness-aware quantization techniques. In Advances in Neural Information Processing Systems, volume 35, pages 9841–9852, 2022. [38] Mike West, Carrie Blanchette, Holly Dressman, Erich Huang, Seiichi Ishida, Rainer Spang, Harry Zuzan, John A. Olson, Jeffrey R. Marks, and Joseph R. Nevins. Predicting the clinical status of human breast cancer by using gene expression profiles. Proceedings of the National Academy of Sciences of the United States of America, 98(20):11462–11467, 2001. [39] Stephen J. Wright. Coordinate descent algorithms. Mathematical Programming, 151(1):3–34, 2015. [40] Tao Yang, Xinlei Yi, Junfeng Wu, Ye Yuan, Di Wu, Ziyang Meng, Yiguang Hong, Hong Wang, Zongli Lin, and Karl H. Johansson. A survey of distributed optimization. Annual Reviews in Control, 47: 278–305, 2019. [41] Ciyou Zhu, Richard H. Byrd, Peihuang Lu, and Jorge Nocedal. Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization. ACM Transactions on Mathematical Software, 23(4):550–560, 1997.

Related work

A.1

Other random-subspace derivative-free and zeroth-order methods

Random-subspace ideas have also been studied in related but different derivative-free and zeroth-order optimization settings. Examples include model-based derivative-free methods in random subspaces [5], directsearch methods using random subspaces [34], zeroth-order methods based on orthogonal random directions [18], and expected-decrease analyses for derivative-free algorithms using random subspaces [14]. These works concern derivative-free and zeroth-order optimization, whereas our focus is on Nesterov-type acceleration under a directional-derivative oracle.

A.2

Relation to accelerated compressed-gradient methods

Several works have studied acceleration in compressed-gradient methods for distributed optimization. These works are related to our goal of reducing the amount of first-order information used per iteration, but their models and resulting oracle-complexity implications are different from ours. First, consider the generic-compressor framework of Li et al. [20]. If one takes C(v) = P P ⊤ v as a compressor and applies their algorithm, then their compression parameter becomes ωLi = ω − 1 under our notation, since E[P P ⊤ ] = Id and E[(P P ⊤ )2 ] ⪯ ωId . Accordingly, their bounds yield oracle complexities s ! r ! L L ∆0 and O rω log O R0 rω ϵ µ ϵ in the convex and strongly convex cases, respectively, under our oracle model. By Proposition 2.3, we have rω ≥ d. Thus, this generic-compressor route is not better than standard full-dimensional Nesterov acceleration in our oracle model. In this sense, if P P ⊤ is used only through a generic unbiased-compressor framework, the benefit of randomized-subspace structure does not appear in the resulting oracle bounds. The work of Safaryan et al. [35] is also related, as it studies smoothness-aware compression and accelerated variants under matrix smoothness. Their setting is distributed: each local loss fi is equipped with a local smoothness matrix Li , which is used together with a random diagonal sketch matrix to define a smoothnessaware compression mechanism. As discussed in their limitations section, this approach requires the server 1/2 to store Li for all workers; hence it is not expected to be practical for large d unless the matrices Li have special structure, such as low-rank or diagonal structure. This differs from our setting, where we allow general randomized subspace sketches Pk and assume only a global matrix smoothness condition for the objective f . In particular, our framework does not require local smoothness matrices for individual workers or the storage of such matrices at the server. Thus, while both approaches exploit matrix smoothness, the optimization model and sketching mechanism are different.

Communication perspective on randomized subspace methods

Consider the distributed optimization problem n

min f (x),

x∈Rd

f (x) :=

1X fi (x), n i=1

where n is the number of workers and fi is the local loss on worker i. This is a standard objective form in distributed optimization; see, e.g., [15, 40]. Each worker sends ∇fi (yk ) ∈ Rd to the server, which forms n

∇f (yk ) =

1X ∇fi (yk ). n i=1 13

To reduce the communication cost per round, we instead use a shared sketch of dimension r. At iteration k, the server broadcasts the random seed defining Pk ∈ Rd×r , and each worker i computes and sends si,k := Pk⊤ ∇fi (yk ) ∈ Rr . The server then averages n

s̄k :=

1X si,k n i=1

and reconstructs gk := Pk s̄k . By linearity, n

gk = Pk

1X ⊤ P ∇fi (yk ) n i=1 k

! = Pk Pk⊤ ∇f (yk ).

Hence the iterates coincide with those of the corresponding single-machine version of the proposed method applied to f . Therefore, all convergence guarantees proved for the single-machine method apply also to the distributed shared-sketch implementation. In particular, the corresponding convex and strongly convex convergence rates carry over to the distributed setting. The matrix smoothness assumption is also natural in this setting. Indeed, if each local loss fi satisfies 1 fi (y) ≤ fi (x) + ⟨∇fi (x), y − x⟩ + (y − x)⊤ Li (y − x), 2

Li ⪰ 0,

then averaging over i = 1, . . . , n yields 1 f (y) ≤ f (x) + ⟨∇f (x), y − x⟩ + (y − x)⊤ L(y − x), 2

L :=

1X Li . n i=1

Such local smoothness matrices Li arise naturally in the matrix-smoothness framework for distributed optimization; see [35].

Missing Proofs for Section 2

Proof of Proposition 2.3. Set A := P P ⊤ . Then A is symmetric positive semidefinite and rank(A) ≤ r. Also, by Assumption 2.1, L ⪰ 0, L ̸= 0, and hence L = ∥L∥ > 0. We prove in order the lower bounds ω ≥ d/r and ℓ ≥ 1 for constants satisfying (3) and (4), respectively, the product bound ℓω ≥ 1, and finally that (4) holds with ℓ = ω whenever ω satisfies (3). (i) Lower bound on ω. Let λ1 , . . . , λd ≥ 0 be the eigenvalues of A. Since at least d − r of them are zero, the Cauchy–Schwarz inequality gives d d 2 X 2 X tr(A) = λi ≤ r λ2i = r tr(A2 ). i=1

i=1

Hence tr(A2 ) ≥

(tr(A))2 . r

Taking expectations, we obtain E[tr(A2 )] ≥

1 E (tr(A))2 . r

Now apply Jensen’s inequality to the convex function t 7→ t2 : 2 E (tr(A))2 ≥ E[tr(A)] .

Using the unbiasedness condition (2), E[A] = Id , we obtain E[tr(A)] = tr(E[A]) = tr(Id ) = d. Therefore, E[tr(A2 )] ≥

d2 . r

On the other hand, taking the trace in (3) gives E[tr(A2 )] = tr(E[A2 ]) ≤ tr(ωId ) = ωd. Combining the lower and upper bounds yields ωd ≥

d2 , r

ω≥

d . r

and hence

(ii) Lower bound on ℓ. Let u ∈ Rd be a unit eigenvector corresponding to the largest eigenvalue L of L, that is, ∥u∥ = 1,

Lu = Lu.

Since L ⪰ 0, we may write its spectral decomposition as L = Luu⊤ +

d X

λi vi vi⊤ ,

λi ≥ 0.

i=2

Therefore, L − Luu⊤ =

d X

λi vi vi⊤ ⪰ 0,

i=2

so L ⪰ Luu⊤ . Therefore, we have ALA ⪰ L Auu⊤ A. Taking the quadratic form with u, we obtain u⊤ ALAu ≥ L u⊤ Auu⊤ Au = L (u⊤ Au)2 . Taking expectation, we get u⊤ E[ALA]u = E[u⊤ ALAu] ≥ L E[(u⊤ Au)2 ]. Again by Jensen’s inequality for the convex function t 7→ t2 , 2 E[(u⊤ Au)2 ] ≥ E[u⊤ Au] . Using once more the unbiasedness condition E[A] = Id , we obtain 2 E[u⊤ Au] = (u⊤ E[A]u)2 = (u⊤ Id u)2 = 1. Therefore, u⊤ E[ALA]u ≥ L. 15

Now let ℓ be any admissible constant in (4). Then E[ALA] ⪯ ℓL Id . Taking again the quadratic form with u, we get u⊤ E[ALA]u ≤ u⊤ (ℓLId )u = ℓL∥u∥2 = ℓL. Combining the lower and upper bounds yields L ≤ ℓL. Since L > 0, we conclude that ℓ ≥ 1. (iii) Product bound. We have

d ≥1 since 1 ≤ r ≤ d. r Together with the lower bound ℓ ≥ 1, this implies that for every admissible pair (ω, ℓ), ω≥

ℓω ≥ 1. (iv) Admissibility of ℓ = ω. Since L ⪰ 0 and L = ∥L∥, we have L ⪯ LId . Hence ALA ⪯ LA2 . Taking expectations and using (3), we obtain E[ALA] ⪯ L E[A2 ] ⪯ ωL Id . That is, E[P P ⊤ LP P ⊤ ] ⪯ ωL Id . Therefore, the admissible constants ℓ in (4) contain ω. Hence one may choose ℓ so that ℓ ≤ ω. This completes the proof.

Missing Proofs for Section 3

Proof of Proposition 3.2. We show that the resulting (xk , yk )-sequence coincides with the classical twosequence Nesterov accelerated gradient method 1 N ∇f (ykN ), xN k+1 := yk − L p 1 + 1 + 4t2k tk+1 := , 2 tk − 1 βk := , tk+1 N N N := xN yk+1 k+1 + βk xk+1 − xk , N with initialization xN 0 = x0 , y0 = x0 , and t0 = 1.

(9) (10) (11) (12)

Since r = d and Pk Pk⊤ = Id for all k, we have almost surely (Pk Pk⊤ )2 = Id ,

Pk Pk⊤ LPk Pk⊤ = L.

Taking expectations, E[P P ⊤ ] = Id , E[(P P ⊤ )2 ] = Id , E[P P ⊤ LP P ⊤ ] = L. Hence Assumption 2.2 is satisfied with the admissible choice ω = 1,

ℓ = 1.

By Algorithm 1, the iterates satisfy 1 , 2Lℓ Ak+1 = Ak + ak+1 , Ak ak+1 yk = xk + zk , Ak+1 Ak+1 1 Pk Pk⊤ ∇f (yk ), xk+1 = yk − Lℓ zk+1 = zk − ak+1 Pk Pk⊤ ∇f (yk ). m=

(13) (14) (15) (16) (17)

We now fix ω = ℓ = 1. By (13), 1 . 2L By the scalar update in Algorithm 1, ak+1 > 0 satisfies m=

m(Ak + ak+1 ) =

ω 2 a . 2 k+1

Hence, since ω = 1, 1 1 (Ak + ak+1 ) = a2k+1 . 2L 2 By (14), this is equivalent to 1 1 Ak+1 = a2k+1 2L 2

⇐⇒

Ak+1 = La2k+1 .

By (15), yk =

Ak ak+1 xk + zk . Ak+1 Ak+1

By (16), together with Pk Pk⊤ = Id and ℓ = 1, xk+1 = yk −

1 1 Pk Pk⊤ ∇f (yk ) = yk − ∇f (yk ), Lℓ L

and (17) becomes zk+1 = zk − ak+1 Pk Pk⊤ ∇f (yk ) = zk − ak+1 ∇f (yk ). Thus, in the case r = d with Pk Pk⊤ = Id , the updates become ak+1 Ak xk + zk , Ak+1 Ak+1 1 xk+1 = yk − ∇f (yk ), L zk+1 = zk − ak+1 ∇f (yk ), yk =

Ak+1 = Ak + ak+1 , 17

Ak+1 = La2k+1 .

(18) (19) (20) (21)

We next analyze the scalar sequence {Ak }. Define p tk := LAk+1 ,

k ≥ 0.

Since A0 = 0, at k = 0 we have from (21) A1 = La21 .

A1 = A0 + a1 = a1 , Hence a1 = La21 , so the nonzero solution is a1 = 1/L. Therefore A1 = A0 + a1 =

1 , L

t0 =

LA1 = 1.

For general k ≥ 0, combining Ak+2 = La2k+2 from (21) with Ak+2 = Ak+1 + ak+2 yields La2k+2 = Ak+1 + ak+2 . Multiplying both sides by L and using t2k = LAk+1 ,

t2k+1 = LAk+2 = L2 a2k+2 ,

we obtain t2k+1 = t2k + Lak+2 . On the other hand, since Ak+2 = La2k+2 , we have Lak+2 =

L2 a2k+2 =

LAk+2 = tk+1 ,

where in the first equality we use ak+2 > 0 to select the positive root. Therefore t2k+1 = t2k + tk+1 , that is, t2k+1 − tk+1 − t2k = 0. For each fixed tk , the above quadratic equation has the positive root p 1 + 1 + 4t2k , k ≥ 0, tk+1 = 2 since tk+1 > 0. This is exactly the recursion (10). We now eliminate zk and derive the two-sequence recursion for (xk , yk ). From (18) we have zk =

Ak+1 yk − Ak xk ak+1

so for all k ≥ 0. Using (19), ∇f (yk ) = L yk − xk+1 , and substituting into the update (20) together with (22) gives zk+1 = zk − ak+1 ∇f (yk ) Ak+1 yk − Ak xk − ak+1 L yk − xk+1 . = ak+1

(22)

Since Ak+1 = La2k+1 , we have Lak+1 = Ak+1 /ak+1 , hence zk+1 =

Ak+1 xk+1 − Ak xk Ak+1 yk − Ak xk Ak+1 − yk − xk+1 = . ak+1 ak+1 ak+1

Next, from (18) with k replaced by k + 1 we obtain Ak+1 ak+2 xk+1 + zk+1 . Ak+2 Ak+2

yk+1 =

Substituting the expression for zk+1 from (23) yields Ak+1 ak+2 Ak+1 xk+1 − Ak xk xk+1 + · Ak+2 Ak+2 ak+1 ak+2 ak+2 Ak Ak+1 1+ xk+1 − xk . = Ak+2 ak+1 ak+1 Ak+2

yk+1 =

Define βk :=

ak+2 Ak . ak+1 Ak+2

Then the coefficient of xk is −βk , and the coefficient of xk+1 is Ak+1 ak+2 ak+1 Ak+1 + ak+2 Ak+1 1+ = Ak+2 ak+1 ak+1 Ak+2 ak+1 (Ak+2 − ak+2 ) + ak+2 (Ak + ak+1 ) = = 1 + βk , ak+1 Ak+2 where the second equality is derived from Ak+2 = Ak+1 + ak+2 and Ak+1 = Ak + ak+1 . Thus, we have yk+1 = xk+1 + βk xk+1 − xk . It remains to express βk in terms of {tk } and identify it with (11). By the definition p (k ≥ 0), tk := LAk+1 we have Ak =

t2k−1 , L

Ak+2 =

t2k+1 L

(k ≥ 1).

Moreover, since Ak+1 = La2k+1 , it follows that tk =

LAk+1 =

ak+1 > 0,

L2 a2k+1 = Lak+1 ,

and hence ak+1 =

tk , L

ak+2 =

tk+1 = Lak+2 , tk+1 . L

Therefore, for k ≥ 1, tk+1 t2k−1 · t2 ak+2 Ak βk = = L 2 L = k−1 . ak+1 Ak+2 tk tk+1 tk tk+1 · L L Under the positivity condition, the recursion (10) is equivalent to t2j+1 − tj+1 − t2j = 0, 19

j ≥ 0.

(23)

In particular, for k ≥ 1, applying this with j = k − 1 gives t2k−1 = t2k − tk . Substituting into the expression for βk yields βk =

t2k − tk tk − 1 = , tk tk+1 tk+1

which coincides with (11) for k ≥ 1. For k = 0 we have A0 = 0, so β0 = 0, and also β0 = (t0 − 1)/t1 = 0, so the formula holds for all k ≥ 0. Summarizing, the pair (xk , yk ) generated by (18)–(21) satisfies xk+1 = yk −

1 ∇f (yk ), L

yk+1 = xk+1 +

tk − 1 xk+1 − xk , tk+1

N with the same initialization x0 = y0 = xN 0 = y0 and the same sequence {tk } as in (9)–(12). Thus (xk , yk ) N obeys exactly the same recursion and initial conditions as (xN k , yk ), and hence

xk = xN k,

yk = ykN

for all k ≥ 0.

This completes the proof. Proof of Theorem 3.3. Define 1 2 Ψk := E Ak f (xk ) − f ⋆ + ∥zk − x⋆ ∥ . 2 Fix k ≥ 0, and let Fk := σ(P0 , . . . , Pk−1 ),

gk := ∇f (yk ).

By Algorithm 1, the iterates satisfy 1 , 2Lℓ Ak+1 = Ak + ak+1 , ak+1 Ak xk + zk , yk = Ak+1 Ak+1 1 xk+1 = yk − Pk Pk⊤ ∇f (yk ), Lℓ zk+1 = zk − ak+1 Pk Pk⊤ ∇f (yk ).

(24) (25) (26) (27)

Moreover, by the scalar update in Algorithm 1, ak+1 > 0 satisfies m(Ak + ak+1 ) =

ω 2 a . 2 k+1

Equivalently, by (24),

ω 2 a . 2 k+1 By (26) and (1) applied with x = yk and y = xk+1 , we have mAk+1 =

1 f (xk+1 ) ≤ f (yk ) + ⟨∇f (yk ), xk+1 − yk ⟩ + (xk+1 − yk )⊤ L(xk+1 − yk ) 2 1 1 ⊤ gk , Pk Pk gk + g ⊤ Pk Pk⊤ LPk Pk⊤ gk . = f (yk ) − Lℓ 2L2 ℓ2 k

(28)

Taking conditional expectation and using (2) and (4), 1 1 2 g ⊤ E[Pk Pk⊤ LPk Pk⊤ | Fk ]gk ∥gk ∥ + Lℓ 2L2 ℓ2 k 1 1 2 2 ≤ f (yk ) − ∥gk ∥ + ∥gk ∥ Lℓ 2Lℓ 1 2 = f (yk ) − ∥gk ∥ 2Lℓ 2 = f (yk ) − m ∥gk ∥ ,

E[f (xk+1 ) | Fk ] ≤ f (yk ) −

(29)

By (27), zk+1 − x⋆ = (zk − x⋆ ) − ak+1 Pk Pk⊤ gk . Hence a2 1 1 2 2 2 ∥zk+1 − x⋆ ∥ = ∥zk − x⋆ ∥ − ak+1 zk − x⋆ , Pk Pk⊤ gk + k+1 Pk Pk⊤ gk . 2 2 2 Taking conditional expectation and using (2), (3), ωa2k+1 1 1 2 2 ⋆ 2 ∥gk ∥ . E ∥zk+1 − x ∥ Fk ≤ ∥zk − x⋆ ∥ − ak+1 ⟨zk − x⋆ , gk ⟩ + 2 2 2

(30)

Subtract f ⋆ from both sides of (29), multiply by Ak+1 , and then add (30). This gives 1 ⋆ ⋆ 2 E Ak+1 f (xk+1 ) − f + ∥zk+1 − x ∥ Fk 2 1 2 ⋆ ≤ Ak+1 f (yk ) − f + ∥zk − x⋆ ∥ − ak+1 ⟨zk − x⋆ , gk ⟩ 2 ωa2k+1 2 − Ak+1 m − ∥gk ∥ . 2 By (28), the last coefficient is zero, hence 1 2 E Ak+1 f (xk+1 ) − f ⋆ + ∥zk+1 − x⋆ ∥ Fk 2 1 2 ≤ Ak+1 f (yk ) − f ⋆ + ∥zk − x⋆ ∥ − ak+1 ⟨zk − x⋆ , gk ⟩ . 2

(31)

From (25), Ak+1 yk = Ak xk + ak+1 zk . Therefore Ak (xk − yk ) + ak+1 (x⋆ − yk ) = ak+1 (x⋆ − zk ). Taking inner product with gk , −ak+1 ⟨zk − x⋆ , gk ⟩ = Ak ⟨xk − yk , gk ⟩ + ak+1 ⟨x⋆ − yk , gk ⟩ .

(32)

By convexity of f , f (xk ) ≥ f (yk ) + ⟨gk , xk − yk ⟩ , ⋆

⋆

f = f (x ) ≥ f (yk ) + ⟨gk , x − yk ⟩ . Multiply (33) by Ak , (34) by ak+1 , and add: Ak ⟨xk − yk , gk ⟩ + ak+1 ⟨x⋆ − yk , gk ⟩ ≤ Ak f (xk ) − f (yk ) + ak+1 f ⋆ − f (yk ) .

(33) (34)

Using (32), we obtain −ak+1 ⟨zk − x⋆ , gk ⟩ ≤ Ak f (xk ) − f (yk ) + ak+1 f ⋆ − f (yk ) . Substitute (35) into (31): 1 ⋆ 2 ⋆ E Ak+1 f (xk+1 ) − f + ∥zk+1 − x ∥ Fk 2 1 2 ≤ Ak+1 f (yk ) − f ⋆ + ∥zk − x⋆ ∥ + Ak f (xk ) − f (yk ) + ak+1 f ⋆ − f (yk ) 2 1 2 = Ak f (xk ) − f (yk ) + (Ak+1 − ak+1 ) f (yk ) − f ⋆ + ∥zk − x⋆ ∥ 2 1 2 = Ak f (xk ) − f ⋆ + ∥zk − x⋆ ∥ , 2

(35)

(36)

By nonnegativity, the tower property applies to (36); hence, taking expectation yields Ψk+1 ≤ Ψk . Iterating gives ΨN ≤ Ψ0 . Since A0 = 0 and z0 = x0 , Ψ0 =

1 2 ∥x0 − x⋆ ∥ . 2

Therefore, AN E[f (xN ) − f ⋆ ] ≤ ΨN ≤ Ψ0 = It remains to lower bound AN . Define ρ :=

1 2 ∥x0 − x⋆ ∥ . 2

(37)

2m > 0. ω

From (28), mAk+1 =

ω 2 a , 2 k+1

Ak+1 =

a2k+1 . ρ

or equivalently,

Since Ak+1 = Ak + ak+1 , we have Ak+1 − Ak = Let

s Bk :=

ρAk+1 .

Ak . ρ

Then 2 Bk+1 − Bk2 = Bk+1 .

Hence (Bk+1 − Bk )(Bk+1 + Bk ) = Bk+1 . Since Ak+1 = Ak + ak+1 with ak+1 > 0, we have Ak+1 > Ak ≥ 0, so Bk+1 > Bk ≥ 0 and in particular Bk+1 + Bk > 0. Thus, 1 Bk+1 Bk+1 − Bk = ≥ . Bk+1 + Bk 2 Since B0 = 0, summing gives BN ≥ 22

N . 2

Therefore, ρN 2 m 2 = N . 4 2ω

2 AN = ρBN ≥

(38)

Combining (37) and (38), for all N ≥ 1, 2

E[f (xN ) − f ⋆ ] ≤

1 ω ∥x0 − x⋆ ∥ 2 . ∥x0 − x⋆ ∥ ≤ 2AN m N2

Since m = 1/(2Lℓ), this gives 2

E[f (xN ) − f ⋆ ] ≤ 2Lωℓ

∥x0 − x⋆ ∥ , N2

which is (5). Finally, let R0 := ∥x0 − x⋆ ∥. Solving 2Lωℓ

R02 ≤ϵ N2

gives r N ≥ R0

2Lωℓ . ϵ

This completes the proof.

Missing Proofs for Section 4

Proof of Proposition 4.2. We show that the resulting (xk , yk )-sequence coincides with the classical twosequence accelerated gradient method 1 ∇f (ykN ), L N N N := xN yk+1 k+1 + β xk+1 − xk ,

N xN k+1 := yk −

where

√ β := √

L− L+

√ √

µ µ

(39) (40)

(41)

N with initialization xN 0 = x0 and y0 = x0 . As in the convex case, the assumptions r = d and Pk Pk⊤ = Id imply

(Pk Pk⊤ )2 = Id ,

Pk Pk⊤ LPk Pk⊤ = L

almost surely, so Assumption 2.2 is satisfied with ω = 1,

ℓ = 1.

By Algorithm 2, the iterates satisfy r µ θ= , Lωℓ 1 θ xk + zk , yk = 1+θ 1+θ 1 xk+1 = yk − Pk Pk⊤ ∇f (yk ), Lℓ θ zk+1 = (1 − θ)zk + θyk − Pk Pk⊤ ∇f (yk ). µ

(42) (43)

(44)

We now fix ω = ℓ = 1. By (42), r θ=

µ . L

Using (43)–(44) together with Pk Pk⊤ = Id , we obtain 1 θ xk + zk , 1+θ 1+θ 1 xk+1 = yk − ∇f (yk ), L θ zk+1 = (1 − θ)zk + θyk − ∇f (yk ). µ yk =

(45) (46) (47)

We first rewrite the z-update. From (46), xk+1 = yk −

1 ∇f (yk ) L

=⇒

∇f (yk ) = L yk − xk+1 .

Since µ = θ2 L, it follows that 1 θ θ ∇f (yk ) = 2 L yk − xk+1 = yk − xk+1 . µ θ L θ Substituting this into the update for zk+1 yields 1 zk+1 = (1 − θ)zk + θyk − yk − xk+1 θ 1 1 = (1 − θ)zk + θ − yk + xk+1 . θ θ

(48)

We now derive a convenient representation of zk in terms of {xj }j≤k . From (45), (1 + θ)yk = xk + θzk . Thus θ(1 − θ)zk = (1 − θ) (1 + θ)yk − xk

= (1 − θ2 )yk − (1 − θ)xk . Equivalently, θ(1 − θ)zk + (θ2 − 1)yk = −(1 − θ)xk .

(49)

Dividing (49) by θ > 0 gives 1−θ 1 yk = − xk . (1 − θ)zk + θ − θ θ On the other hand, from (48) we have 1 1 zk+1 = (1 − θ)zk + θ − yk + xk+1 . θ θ Substituting (50) into (51) yields 1−θ 1 xk + xk+1 θ θ 1−θ = xk+1 + xk+1 − xk . θ

zk+1 = −

Using this representation of zk+1 in (45) at step k + 1, we obtain 1 θ 1 θ 1−θ yk+1 = xk+1 + zk+1 = xk+1 + xk+1 + (xk+1 − xk ) , 1+θ 1+θ 1+θ 1+θ θ 24

(50)

(51)

and therefore

1 θ 1−θ yk+1 = + xk+1 + (xk+1 − xk ) 1+θ 1+θ 1+θ 1−θ = xk+1 + (xk+1 − xk ). 1+θ Using θ =

µ/L, we can rewrite √ β=√

L− L+

√ √

µ µ

1−

µ/L µ/L

1−θ . 1+θ

Hence yk+1 = xk+1 + β(xk+1 − xk ), and the x-update (46) is exactly (39) with ykN = yk . Thus the pair (xk , yk ) generated by (45)–(47) satisfies the same recursion (39)–(40) and the same initialN ization x0 = y0 = xN 0 = y0 as the classical accelerated gradient method. Combining the base case k = 0 with the inductive step, we have shown that xN k = xk ,

ykN = yk

for all k ≥ 0.

This completes the proof of the proposition. Proof of Theorem 4.3. First, we show that the parameter θ in Algorithm 2 satisfies θ ∈ (0, 1]. Since L ⪰ 0 and L = ∥L∥, we have L ⪯ LId . Hence Assumption 2.1 implies f (y) ≤ f (x) + ⟨∇f (x), y − x⟩ +

L 2 ∥y − x∥ . 2

Since f is also µ-strongly convex by Assumption 4.1, it follows that µ ≤ L. Moreover, by Proposition 2.3, ωℓ ≥ 1. Therefore, using the definition of θ in Algorithm 2, 0 < θ2 = and hence θ ∈ (0, 1]. Define

µ µ ≤ ≤ 1, Lωℓ L

h i µ 2 Φk := E f (xk ) − f ⋆ + ∥zk − x⋆ ∥ . 2

Fix k ≥ 0, and let Fk := σ(P0 , . . . , Pk−1 ),

gk := ∇f (yk ).

By Algorithm 2, the iterates satisfy r µ θ= , Lωℓ 1 θ yk = xk + zk , 1+θ 1+θ 1 xk+1 = yk − Pk Pk⊤ ∇f (yk ), Lℓ θ zk+1 = (1 − θ)zk + θyk − Pk Pk⊤ ∇f (yk ). µ 25

(52) (53) (54) (55)

As shown above, θ ∈ (0, 1]. In particular, θ > 0, so division by θ below is valid, and 1 − θ ≥ 0. By (54) and (1) applied with x = yk and y = xk+1 , we get 1 f (xk+1 ) ≤ f (yk ) + ⟨∇f (yk ), xk+1 − yk ⟩ + (xk+1 − yk )⊤ L(xk+1 − yk ) 2 1 1 ⊤ g ⊤ Pk Pk⊤ LPk Pk⊤ gk . = f (yk ) − gk , Pk Pk gk + Lℓ 2L2 ℓ2 k Taking conditional expectation and using (2) and (4), 1 1 2 g ⊤ E[Pk Pk⊤ LPk Pk⊤ | Fk ]gk ∥gk ∥ + Lℓ 2L2 ℓ2 k 1 2 ≤ f (yk ) − ∥gk ∥ . 2Lℓ

E[f (xk+1 ) | Fk ] ≤ f (yk ) −

Define m :=

1 . 2Lℓ

(56)

Then 2

E[f (xk+1 ) | Fk ] ≤ f (yk ) − m ∥gk ∥ .

(57)

wk := (1 − θ)(zk − x⋆ ) + θ(yk − x⋆ ).

(58)

Define

From (55), θ Pk Pk⊤ gk − x⋆ µ θ = (1 − θ)(zk − x⋆ ) + θ(yk − x⋆ ) − Pk Pk⊤ gk µ θ = wk − Pk Pk⊤ gk . µ

zk+1 − x⋆ = (1 − θ)zk + θyk −

Hence µ θ2 µ 2 2 2 ∥zk+1 − x⋆ ∥ = ∥wk ∥ − θ wk , Pk Pk⊤ gk + Pk Pk⊤ gk . 2 2 2µ Taking conditional expectation and using (2), (3), hµ i µ ωθ2 2 2 2 ∥gk ∥ . E ∥zk+1 − x⋆ ∥ Fk ≤ ∥wk ∥ − θ ⟨wk , gk ⟩ + 2 2 2µ

(59)

Summing (57) and (59), and subtracting f ⋆ , h i µ 2 E f (xk+1 ) − f ⋆ + ∥zk+1 − x⋆ ∥ Fk 2 µ ωθ2 2 2 ≤ f (yk ) − f ⋆ + ∥wk ∥ − θ ⟨wk , gk ⟩ − m − ∥gk ∥ . 2 2µ By (52) and (56), m=

ωθ2 , 2µ

so the ∥gk ∥ term vanishes. Thus h i µ 2 E f (xk+1 ) − f ⋆ + ∥zk+1 − x⋆ ∥ Fk 2 µ 2 ⋆ ≤ f (yk ) − f + ∥wk ∥ − θ ⟨wk , gk ⟩ . 2 26

(60)

From (53), (1 + θ)yk = xk + θzk

=⇒

θ(zk − yk ) = yk − xk .

Let ak := zk − x⋆ ,

bk := yk − x⋆ .

Then, recalling the definition (58), we have wk = (1 − θ)ak + θbk . Also, µ µ 2 2 ∥wk ∥ = ∥(1 − θ)ak + θbk ∥ 2 2 µ µ µ 2 2 2 = (1 − θ) ∥ak ∥ + θ ∥bk ∥ − θ(1 − θ) ∥ak − bk ∥ . 2 2 2

(61)

Also, −θ ⟨wk , gk ⟩ = −θ ⟨(1 − θ)ak + θbk , gk ⟩ = −θ(1 − θ) ⟨ak , gk ⟩ − θ2 ⟨bk , gk ⟩ = −θ ⟨bk , gk ⟩ − θ(1 − θ) ⟨ak − bk , gk ⟩ = −θ ⟨bk , gk ⟩ + (1 − θ) ⟨xk − yk , gk ⟩ ,

(62)

where in the last step we used ak − bk = zk − yk = (yk − xk )/θ. Since f is µ-strongly convex, for any u, v, f (u) ≥ f (v) + ⟨∇f (v), u − v⟩ +

µ 2 ∥u − v∥ . 2

Applying this with (u, v) = (x⋆ , yk ) gives ⟨gk , bk ⟩ = ⟨gk , yk − x⋆ ⟩ ≥ f (yk ) − f ⋆ +

µ µ 2 2 ∥yk − x⋆ ∥ = f (yk ) − f ⋆ + ∥bk ∥ . 2 2

(63)

Applying it with (u, v) = (xk , yk ) gives ⟨gk , xk − yk ⟩ ≤ f (xk ) − f (yk ) −

µ 2 ∥xk − yk ∥ . 2

(64)

Substitute (61), (62), (63), (64) into (60): h i µ 2 E f (xk+1 ) − f ⋆ + ∥zk+1 − x⋆ ∥ Fk 2 µ µ µ 2 2 2 ⋆ ≤ f (yk ) − f + (1 − θ) ∥ak ∥ + θ ∥bk ∥ − θ(1 − θ) ∥ak − bk ∥ 2 2 2 µ µ 2 2 − θ f (yk ) − f ⋆ + ∥bk ∥ + (1 − θ) f (xk ) − f (yk ) − ∥xk − yk ∥ 2 2 µ µ 2 2 ⋆ ≤ f (yk ) − f + (1 − θ) ∥ak ∥ + θ ∥bk ∥ 2 2 µ 2 ⋆ − θ f (yk ) − f + ∥bk ∥ + (1 − θ)(f (xk ) − f (yk )) . 2 2

The coefficients of f (yk ) and ∥bk ∥ cancel. Using also ak = zk − x⋆ , we obtain h i µ 2 E f (xk+1 ) − f ⋆ + ∥zk+1 − x⋆ ∥ Fk 2 µ 2 ≤ (1 − θ) f (xk ) − f ⋆ + (1 − θ) ∥zk − x⋆ ∥ 2 µ 2 = (1 − θ) f (xk ) − f ⋆ + ∥zk − x⋆ ∥ . 2 27

(65)

By nonnegativity, the tower property applies to (65); hence, taking expectation gives Φk+1 ≤ (1 − θ)Φk . Iterating yields ΦN ≤ (1 − θ)N Φ0 . Moreover, E[f (xN ) − f ⋆ ] ≤ ΦN . By µ-strong convexity, f (x0 ) − f ⋆ ≥

µ 2 ∥x0 − x⋆ ∥ . 2

Since z0 = x0 , this implies Φ0 = f (x0 ) − f ⋆ +

µ 2 ∥z0 − x⋆ ∥ ≤ 2 f (x0 ) − f ⋆ = 2∆0 . 2

Therefore, E[f (xN ) − f ⋆ ] ≤ ΦN ≤ (1 − θ)N Φ0 ≤ 2(1 − θ)N ∆0 which is (6). It remains to derive the complexity bound. Since 1 − θ ≤ e−θ , E[f (xN ) − f ⋆ ] ≤ 2e−N θ ∆0 . Thus, to guarantee E[f (xN ) − f ⋆ ] ≤ ϵ, it suffices that 2∆0 1 log . θ ϵ

N≥ Using the definition of θ in Algorithm 2, this is s N≥

Lωℓ 2∆0 log . µ ϵ

This completes the proof.

Missing Proofs for Section 5

Proof of Proposition 5.1. We verify the three parts of Assumption 2.2 for each sketch: (i) unbiasedness E[P P ⊤ ] = Id , (ii) the second-moment bound E[(P P ⊤ )2 ] ⪯ ωId , and (iii) the matrix-smoothness interaction E[P P ⊤ LP P ⊤ ] ⪯ ℓLId . d×r Haar sketch. Let be the matrix formed by the first r columns pU be Haar-distributed on O(d), let R ∈ R of U , and set P = d/r R. It is standard that E[RR⊤ ] = dr Id , so E[P P ⊤ ] = Id , verifying unbiasedness. For the second moment, d (P P ⊤ )2 = P P ⊤ P P ⊤ = P P ⊤ . r

Taking expectations and using E[P P ⊤ ] = Id gives E[(P P ⊤ )2 ] = dr Id , so we may take ωHaar = dr . For the interaction with L, write L := ∥L∥,

reff :=

tr(L) . ∥L∥

By results in [8], E[P P ⊤ LP P ⊤ ] =

d β (1 − β)L + tr(L) Id , r r 28

d(d−r) where β = (d+2)(d−1) . Using L ⪰ 0 and ∥L∥ = L,

d reff β tr(L) d E[P P LP P ] ⪯ 1−β+β L Id . (1 − β) + L Id = r r L r d ⊤

⊤

Thus (4) in Assumption 2.2 holds with ℓHaar =

d reff 1−β+β . r d

Coordinate sketch. q Let S ∈ Rd×r consist of r distinct columns sampled uniformly from the identity matrix Id , and let P := dr S. Then SS ⊤ is a coordinate projection, and it is straightforward to check that E[SS ⊤ ] = dr Id , so E[P P ⊤ ] = Id . For the second moment, (P P ⊤ )2 = P P ⊤ P P ⊤ =

d P P ⊤. r

Taking expectations, we obtain E[(P P ⊤ )2 ] = dr Id , so we may take ωCoord = dr . For the interaction with L, define ∥ diag(L)∥ δdiag := . ∥L∥ By a result in [8], E[P P ⊤ LP P ⊤ ] =

d r

r−1 d−r L+ diag(L) . d−1 d−1

Using ∥L∥ = L and ∥ diag(L)∥ = δdiag L, d r−1 d−r d r−1 d−r E[P P ⊤ LP P ⊤ ] ⪯ L+ δdiag L Id = + δdiag L Id . r d−1 d−1 r d−1 d−1 Thus we may take ℓCoord =

d r

r−1 d−r + δdiag . d−1 d−1

Gaussian sketch. Finally, assume P ∈ Rd×r has i.i.d. entries Pij ∼ N (0, 1/r). Then E[P P ⊤ ] = Id by √ := construction. Let G r P , so that Gij ∼ N (0, 1) i.i.d. and P P ⊤ = 1r GG⊤ . For the second moment, a result from [8] implies E[GG⊤ Id GG⊤ ] = r(r + 1)Id + r tr(Id ) Id = r(r + 1)Id + rd Id . Hence E[(P P ⊤ )2 ] =

1 d+r+1 1 E[(GG⊤ )2 ] = 2 r(r + 1) + rd Id = Id , 2 r r r

so ωGauss = d+r+1 . r For the interaction with L, a result in [8] implies E[GG⊤ LGG⊤ ] = r(r + 1)L + r tr(L) Id . √ Recalling that P = G/ r and hence P P ⊤ = 1r GG⊤ , we obtain E[P P ⊤ LP P ⊤ ] =

1 r+1 tr(L) E[GG⊤ LGG⊤ ] = L+ Id . 2 r r r

Writing L := ∥L∥,

reff :=

tr(L) , ∥L∥

and using L ⪰ 0 and ∥L∥ = L, we obtain r+1 reff L LId + Id r r r + 1 + reff = L Id . r

E[P P ⊤ LP P ⊤ ] ⪯

Thus, r + 1 + reff . r It remains to prove the stated basic bounds on reff and δdiag . Let λ1 , . . . , λd be the eigenvalues of L. Since L ⪰ 0 and L = ∥L∥, we have ℓGauss =

0 ≤ λi ≤ L

(i = 1, . . . , d),

and since L ̸= 0, at least one eigenvalue equals L. Therefore L ≤ tr(L) =

d X

λi ≤ dL,

i=1

which gives 1 ≤ reff =

tr(L) ≤ d. L

Next, since δdiag = ∥ diag(L)∥/L, we have δdiag = Because L ⪰ 0, for each i,

max1≤i≤d Lii . L

Lii = e⊤ i Lei ≤ ∥L∥ = L,

so δdiag ≤ 1. Also, d

max Lii ≥

1≤i≤d

1X tr(L) L reff Lii = = . d i=1 d d

Dividing by L gives δdiag ≥

reff . d

Since reff ≥ 1, this also implies δdiag ≥

1 . d

This completes the proof. √ Proof of Proposition 5.2. We first prove that ωℓr2 is minimized at r = 1 for each sketch. For the Haar sketch, by Table 2, r p reff d(d − r) 2 ωHaar ℓHaar r = d 1 − β + β , β= . d (d + 2)(d − 1) Since 1−β+β

reff (d − r)(d − reff ) =1− , d (d + 2)(d − 1)

and d − reff ≥ 0 by Proposition 5.1, the quantity inside the square root is nondecreasing in r. Hence p ωHaar ℓHaar r2

is minimized at r = 1. For the Coordinate sketch, p

r ωCoord ℓCoord

r2 = d

r r−1 r−1 d−r + δdiag = d δdiag + (1 − δdiag ). d−1 d−1 d−1

Since δdiag ≤ 1 by Proposition 5.1, the quantity inside the square root is nondecreasing in r. Hence p ωCoord ℓCoord r2 is minimized at r = 1. For the Gaussian sketch, p

ωGauss ℓGauss r2 =

(d + r + 1)(r + 1 + reff ).

Both factors are increasing in r, so this quantity is increasing in r, and therefore it is minimized at r = 1. Substituting r = 1 into the three expressions gives r p p reff + 2 , QC = d δdiag , QG = (d + 2)(reff + 2). QH = d d+2 By Proposition 5.1, 1 ≤ reff ≤ d,

1 ≤ δdiag ≤ 1, d

δdiag ≥

reff . d

Therefore, p √ 3 d ≤ QC ≤ d, 3(d + 2) ≤ QG ≤ d + 2. ≤ QH ≤ d, d+2 Next, the Haar and Gaussian factors satisfy p (d + 2)(reff + 2) QG 2 = p =1+ , QH d d (reff + 2)/(d + 2) r

and hence

QG =

2 d

QH .

For the Haar–Coordinate comparison, using δdiag ≥ reff /d, we obtain Q2C = d2 δdiag ≥ d reff . On the other hand, Q2H = d2 Thus

reff + 2 ≤ d(reff + 2). d+2

Q2H reff + 2 ≤ ≤ 3, Q2C reff

where the last inequality follows from reff ≥ 1. Hence √ QH ≤ 3 QC . Finally, if L = e1 e⊤ 1 , then L = 1,

reff =

tr(L) = 1, L

δdiag =

∥ diag(L)∥ = 1. L

Therefore, r QC = d,

QH = d

This completes the proof. 31

3 = d+2

3 QC . d+2

Experimental Details

G.1

Quadratic instances

We describe the four quadratic instances used in Section 6.1. In all cases, f (x) = Convex diagonal.

1 ⊤ x Lx, 2

d = 1000,

f ⋆ = 0.

Let L = diag 1,

1 1 ,..., ,0 . d−2 d−2 Then L ⪰ 0 but is not positive definite. Its eigenvalues are 1 1 1, ,..., ,0 . d−2 d−2 Hence L = ∥L∥ = 1,

tr(L) = 1 + (d − 2)

1 = 2, d−2

reff =

tr(L) = 2. ∥L∥

Moreover, since the maximum diagonal entry of L is 1, δdiag = Convex dense.

∥ diag(L)∥ = 1. ∥L∥

Let 1 u := √ (1, −1, 1, −1, . . . , 1, −1)⊤ ∈ Rd . d

1 := (1, . . . , 1)⊤ ∈ Rd , Since d is even, u⊤ 1 = 0. Define L=

1 1 1 Id + 1 − 11⊤ . uu⊤ − d−2 d−2 d(d − 2)

Because u and 1 are orthogonal, the eigenspaces split as follows. Along u, the eigenvalue is 1 1 + 1− = 1. d−2 d−2 Along 1, the eigenvalue is 1 1 − d = 0. d − 2 d(d − 2) On the orthogonal complement of span{u, 1}, the eigenvalue is 1 . d−2 Thus the eigenvalues are 1,

1 1 ,..., ,0 , d−2 d−2

and therefore L = ∥L∥ = 1,

tr(L) = 1 + (d − 2)

1 = 2, d−2

reff =

Each diagonal entry equals 1 1 + d−2 d

and hence ∥ diag(L)∥ =

1− 2 , d

1 d−2

−

δdiag = 32

1 2 = , d(d − 2) d ∥ diag(L)∥ 2 = . ∥L∥ d

tr(L) = 2. ∥L∥

Strongly convex diagonal.

Let

1 1 L = diag 1, ,..., d−1 d−1 Then the eigenvalues of L are

1 1 ,..., 1, d−1 d−1

Therefore, µ = λmin (L) =

1 , d−1

L = ∥L∥ = λmax (L) = 1.

Moreover, tr(L) = 1 + (d − 1)

1 = 2, d−1

reff =

tr(L) = 2. ∥L∥

Since the maximum diagonal entry of L is 1, we have δdiag = Strongly convex dense.

∥ diag(L)∥ = 1. ∥L∥

Let 1 d−2 Id + 11⊤ . d−1 d(d − 1)

The vector 1 is an eigenvector with eigenvalue 1 d−2 + · d = 1, d − 1 d(d − 1) while every vector orthogonal to 1 is an eigenvector with eigenvalue 1 . d−1 Therefore, the eigenvalues of L are

Hence µ = λmin (L) =

1 1 ,..., 1, d−1 d−1

1 , d−1

L = ∥L∥ = λmax (L) = 1.

Moreover, tr(L) = 1 + (d − 1)

1 = 2, d−1

reff =

tr(L) = 2. ∥L∥

Finally, each diagonal entry equals 1 d−2 2 + = , d − 1 d(d − 1) d and hence ∥ diag(L)∥ =

2 , d

δdiag =

∥ diag(L)∥ 2 = . ∥L∥ d

G.2

Effect of the sketch dimension

In the main quadratic experiments, we fixed the sketch dimension to r = 1. Here we additionally examine how the sketch dimension affects oracle-axis convergence. More precisely, for each of the four quadratic instances and for each sketch family (Haar, Block-coordinate, and Gaussian), we compare the proposed method with r ∈ {1, 10, 100} against the corresponding full-dimensional accelerated method under the same oracle budget 10,000. All other settings are unchanged: d = 1000, 10 random seeds, and independent Gaussian initialization x0 ∼ N (0, Id ). The purpose of this experiment is to examine whether the theoretical prediction that r = 1 is oracleoptimal is also reflected numerically in the quadratic examples. For convenience, we recall the four quadratic instances and their associated quantities (reff , δdiag ): • Convex diagonal: L = diag 1,

1 1 ,..., ,0 , d−2 d−2

reff = 2,

δdiag = 1.

• Convex dense: 1 1 1 L= Id + 1 − uu⊤ − 11⊤ , d−2 d−2 d(d − 2)

reff = 2,

δdiag =

2 . d

• Strongly convex diagonal:

1 1 L = diag 1, ,..., d−1 d−1

reff = 2,

δdiag = 1.

• Strongly convex dense: L=

1 d−2 Id + 11⊤ , d−1 d(d − 1)

reff = 2,

δdiag =

2 . d

The results are shown√in Figure 3. As summarized in Table 2, the oracle-complexity comparison is governed by the quantity ωℓr2 , and Proposition 5.2 shows that this quantity is minimized at r = 1 for the Haar, Block-coordinate, and Gaussian sketches. The plots are broadly consistent with this prediction. In all Haar and Gaussian panels, and also in the Block-coordinate panels for the two dense instances, the oracle-axis performance deteriorates as r increases from 1 to 10 and 100. The only clear exceptions are the Block-coordinate panels for the convex diagonal and strongly convex diagonal instances; see Figures 3b and 3h. This behavior is also explained by Table 2. For the Block-coordinate sketch, r p r−1 d−r 2 ωCoord ℓCoord r = d + δdiag . d−1 d−1 For the two diagonal instances considered here, we have δdiag = 1, and hence r p r−1 d−r 2 ωCoord ℓCoord r = d + = d, d−1 d−1 independently of r. Therefore, in these two cases, the theory itself does not predict any oracle-complexity improvement from taking a smaller sketch dimension, which explains why the separation among r = 1, 10, 100 is weak in Figures 3b and 3h. Overall, except for these two Block-coordinate diagonal panels, the numerical results are well aligned with the theoretical prediction that r = 1 is oracle-optimal in the present quadratic examples.

(a) Convex diagonal / Haar

(b) Convex coordinate

diagonal

Block-

(d) Convex dense / Haar

(e) Convex dense / Block-coordinate

(f) Convex dense / Gaussian

(g) Strongly convex diagonal / Haar

(h) Strongly convex diagonal / Blockcoordinate

(i) Strongly convex diagonal / Gaussian

(j) Strongly convex dense / Haar

(k) Strongly convex dense / Blockcoordinate

(l) Strongly convex dense / Gaussian

Figure 3: Oracle-axis convergence in the sketch-dimension scan. The horizontal axis shows the cumulative number of oracle calls, and the vertical axis shows the objective gap f (xk ) − f ⋆ on a logarithmic scale. Each row corresponds to one quadratic instance, and each column corresponds to one sketch family. Each panel compares RS-NAG with r ∈ {1, 10, 100} and the corresponding full-dimensional accelerated method.

Table 4: Dataset-dependent quantities for the logistic-regression datasets. Here d is the ambient dimension, n is the number of training samples, and QH , QG , QC denote the r = 1 constants defined in Proposition 5.2 for the Haar, Gaussian, and Coordinate sketches, respectively. Dataset

G.3

phishing 68 a9a 123 w8a 300 mushroom 112 ijcnn1 22 splice 60

15.3670 22.5681 43.7200 21.0717 13.1196 13.5925

reff

15.8190 15.1421 1.5749 22.9350 47.8971 2.2081 44.0115 109.3820 4.4139 21.4480 34.8299 2.0352 14.3123 17.5017 6.5350 14.0456 9.5503 1.1819

δdiag

0.0496 11055 0.1516 32561 0.1329 49749 0.0967 8124 0.6329 49990 0.0253 1000

Matrix smoothness for logistic regression

We derive the matrix smoothness bound used in Section 6.2. Recall the ℓ2 -regularized logistic regression objective n µ 1X log 1 + exp(−yi a⊤ ∥x∥22 , µ > 0, f (x) = i x) + n i=1 2 where ai ∈ Rd and yi ∈ {−1, +1}. Its gradient is n

∇f (x) = −

yi a i 1X + µx. n i=1 1 + exp(yi a⊤ i x)

Let σ(t) = 1/(1 + e−t ). The Hessian is ∇2 f (x) =

n 1X ⊤ ⊤ σ(yi a⊤ i x) 1 − σ(yi ai x) ai ai + µId . n i=1

Since 0 ≤ σ(t)(1 − σ(t)) ≤ we have

1 4

for all t ∈ R,

1 X ai a⊤ i + µId . 4n i=1 Pn ⊤ ⊤ If A ∈ Rn×d denotes the data matrix whose ith row is a⊤ i , then i=1 ai ai = A A. Thus we may take ∇2 f (x) ⪯

1 ⊤ A A + µId . 4n

This gives the matrix smoothness bound used in the logistic-regression experiments.

G.4

Additional logistic-regression results on standard benchmarks

Experimental setup. We evaluate the strongly convex logistic-regression setting on six real-world binaryclassification benchmarks: phishing, a9a, w8a, mushroom, ijcnn1, and splice. The UCI-derived datasets are cited collectively through the UCI Machine Learning Repository [16]; for w8a and ijcnn1, we follow the standard benchmark attributions [30, 31]. For each dataset, we use the ℓ2 -regularized logistic objective in (7), with the regularization parameter set to µ = 1/n, where n is the number of training samples. We compare the full-gradient methods GD and NAG-SC with the randomized-subspace methods RS-GD and RS-NAG-SC. For the randomized-subspace methods, we use three sketch families: Haar, coordinate, and Gaussian. We use 10 random seeds, and initialize each run from a Gaussian random vector. For logistic regression, we use L=

1 ⊤ A A + µId , 4n 36

(a) phishing

(b) a9a

(d) mushroom

(e) ijcnn1

(f) splice

Figure 4: Oracle-axis comparison for ℓ2 -regularized logistic regression on six real-world datasets. The horizontal axis shows oracle calls, and the vertical axis shows f (xk ) − fref on a logarithmic scale, where fref is a reference value computed by L-BFGS-B. We compare GD, NAG-SC, RS-GD, and RS-NAG-SC with Haar, coordinate, and Gaussian sketches. Each plotted curve is the mean over 10 random seeds, and the shaded region shows one standard deviation. For each dataset, µ = 1/n, and the oracle budget is 100,000 and set L = ∥L∥, computed numerically as the largest eigenvalue of L. In the randomized-subspace methods, we use the r = 1 setting throughout, in accordance with the theoretical comparison developed above. Dataset-dependent quantities. For each dataset, we also compute the matrix-smoothness-derived quantities reff and δdiag , together with the corresponding r = 1 constants QH , QG , and QC for the Haar, Gaussian, and coordinate sketches, respectively, as defined in Proposition 5.2. These values are reported in Table 4. Discussion. We observe that RS-NAG-SC consistently achieves strong performance across all six datasets. Moreover, the relative convergence behavior among the Haar, Gaussian, and Coordinate sketches is broadly consistent with the dataset-dependent Q values reported in Table 4. Indeed, on datasets such as a9a, w8a, mushroom, and ijcnn1, where QH and QG are relatively small, the Haar and Gaussian variants tend to converge faster than the Coordinate variant and also faster than NAG-SC. On the other hand, for phishing, where the Coordinate constant is also relatively small, the Coordinate variant performs comparably well to the Haar and Gaussian variants. For splice, where the Coordinate constant is particularly small, the Coordinate variant performs better. Overall, these results suggest that the Q values can serve as a useful practical guide when choosing the sketch distribution before running the method.

G.5

Experimental resources, dataset sources, and terms of use

Computational resources. All reported logistic-regression experiments were run on a single NVIDIA RTX A5000 GPU with 24564 MiB of memory and CUDA support. The machine had an AMD EPYC 7413 CPU and 503 GiB of system memory. The software environment used Python 3.10.12, NumPy 2.2.6, SciPy 1.15.3, scikit-learn 1.7.2, Matplotlib 3.10.8, and PyTorch 2.11.0 with CUDA support. The main and additional logistic-regression experiments each took about 70 hours on this machine. Thus, the reported logistic-regression experiments required about 140 GPU-hours in total. The reported quadratic experiments

were run separately on a local personal computer using CPU only. The local machine had an Apple M1 processor and 8 GB of memory. The main and appendix quadratic experiments each completed within one wall-clock hour in our runs. Dataset sources and terms of use. We use only previously released public benchmark datasets and do not redistribute any dataset. The datasets hiva agnostic and bioresponse were obtained from OpenML [36] with dataset IDs 1039 and 46912, respectively; their OpenML license fields are listed as Public and Public Domain, respectively. The remaining datasets were obtained from the LIBSVM binary-classification dataset page [6], which provides LIBSVM-formatted versions of datasets from existing public benchmark collections and lists source and preprocessing information for each dataset. We cite the corresponding original source papers or dataset records whenever available. We checked the license or terms of use stated on the corresponding source or access pages whenever available. We use the datasets only for research evaluation and do not redistribute any dataset.

Oracle complexity of the basic randomized-subspace gradient method

In Table 1, the row labeled “RS-GD (Kozak et al.)” refers to the basic randomized-subspace gradient iteration xk+1 = xk − η Pk Pk⊤ ∇f (xk ), k ≥ 0, (66) introduced by Kozak et al. [17]. The exact oracle-complexity expressions shown in Table 1 are not stated in this form in Kozak et al. [17]; for convenience, we record below short derivations. Proposition H.1 (Convex rate for the basic randomized-subspace gradient method). Suppose Assumptions 2.1 to 2.2 hold, with ℓ chosen so that ℓ ≤ ω, as assumed throughout the paper. Let {xk }k≥0 be generated by (66) with the constant step-size 1 . 2ωL

η=

Then the expected objective values are nonincreasing: E[f (xk+1 )] ≤ E[f (xk )]

for all k ≥ 0.

Moreover, for every N ≥ 1, 2

2ωL ∥x0 − x⋆ ∥ . N Consequently, defining R0 := ∥x0 − x⋆ ∥, it suffices to take E[f (xN ) − f ⋆ ] ≤

N≥

(67)

2ωLR02 ϵ

to guarantee E[f (xN ) − f ⋆ ] ≤ ϵ. Since each iteration uses r oracle calls, the oracle complexity is L #Oracle = O rωR02 . ϵ

(68)

Proof. Let gk := ∇f (xk ). Since L ⪰ 0 and L = ∥L∥, we have L ⪯ LId . Therefore Assumption 2.1 implies the standard L-smoothness inequality L 2 f (y) ≤ f (x) + ⟨∇f (x), y − x⟩ + ∥y − x∥ for all x, y ∈ Rd . 2 38

In particular, since f is convex and attains its minimum at x⋆ , the standard smooth-convex inequality gives 2 ∥∇f (x)∥ ≤ 2L f (x) − f ⋆ for all x ∈ Rd . (69) We first show monotonicity of the expected objective. Applying Assumption 2.1 with x = xk and y = xk+1 = xk − ηPk Pk⊤ gk , we obtain f (xk+1 ) ≤ f (xk ) − η gk⊤ Pk Pk⊤ gk +

η2 ⊤ g Pk Pk⊤ LPk Pk⊤ gk . 2 k

Taking conditional expectation with respect to Fk , and using E[Pk Pk⊤ | Fk ] = Id ,

E[Pk Pk⊤ LPk Pk⊤ | Fk ] ⪯ ℓLId ,

we obtain 2

E[f (xk+1 ) | Fk ] ≤ f (xk ) − η ∥gk ∥ +

ηℓL η 2 ℓL 2 2 ∥gk ∥ . ∥gk ∥ = f (xk ) − η 1 − 2 2

With η = 1/(2ωL) and ℓ ≤ ω, we have η=

1 1 ≤ , 2ωL 2ℓL

and hence 1−

ηℓL ≥ 0. 2

Therefore, E[f (xk+1 ) | Fk ] ≤ f (xk ). Next, expanding the squared distance to x⋆ , we obtain 2

∥xk+1 − x⋆ ∥ = ∥xk − x⋆ ∥ − 2η xk − x⋆ , Pk Pk⊤ gk + η 2 gk⊤ (Pk Pk⊤ )2 gk . Taking conditional expectation and using E[Pk Pk⊤ | Fk ] = Id ,

E[(Pk Pk⊤ )2 | Fk ] ⪯ ωId ,

we get 2

E[∥xk+1 − x⋆ ∥ | Fk ] ≤ ∥xk − x⋆ ∥ − 2η ⟨xk − x⋆ , gk ⟩ + η 2 ω ∥gk ∥ . By convexity, ⟨xk − x⋆ , gk ⟩ ≥ f (xk ) − f ⋆ , and by (69), 2

∥gk ∥ ≤ 2L(f (xk ) − f ⋆ ). Therefore 2 2 E[∥xk+1 − x⋆ ∥ | Fk ] ≤ ∥xk − x⋆ ∥ − 2η − 2η 2 ωL f (xk ) − f ⋆ . Taking expectation and using the tower property yields 2 2 E[∥xk+1 − x⋆ ∥ ] ≤ E[∥xk − x⋆ ∥ ] − 2η − 2η 2 ωL E[f (xk ) − f ⋆ ]. Summing for k = 0, . . . , N − 1, we obtain −1 NX 2 2η − 2η 2 ωL E[f (xk ) − f ⋆ ] ≤ ∥x0 − x⋆ ∥ . k=0

Since E[f (xk )] is nonincreasing, ⋆

N E[f (xN ) − f ] ≤

N −1 X

E[f (xk ) − f ⋆ ] ≤

k=0

∥x0 − x⋆ ∥ . 2η − 2η 2 ωL

Using η = 1/(2ωL), we have 2η − 2η 2 ωL =

1 . 2ωL

Therefore

2ωL ∥x0 − x⋆ ∥ , N which proves (67). The oracle bound (68) follows immediately because each iteration uses r oracle calls. E[f (xN ) − f ⋆ ] ≤

Proposition H.2 (Strongly convex rate for the basic randomized-subspace gradient method). Suppose Assumptions 2.1 to 2.2 hold, and let {xk }k≥0 be generated by (66) with the constant step-size η= Then, for every k ≥ 0,

1 . ℓL

µ E[f (xk+1 ) − f ⋆ ] ≤ 1 − E[f (xk ) − f ⋆ ]. ℓL

(70)

µN µ N ⋆ f (x0 ) − f ⋆ . f (x0 ) − f ≤ exp − E[f (xN ) − f ] ≤ 1 − ℓL ℓL

(71)

Hence, for every N ≥ 0, ⋆

Consequently, defining ∆0 := f (x0 ) − f ⋆ , it suffices to take N≥

ℓL ∆0 log µ ϵ

to guarantee E[f (xN ) − f ⋆ ] ≤ ϵ. Since each iteration uses r oracle calls, the oracle complexity is L ∆0 #Oracle = O rℓ log . µ ϵ

(72)

Proof. Let gk := ∇f (xk ). Applying Assumption 2.1 with x = xk and y = xk+1 = xk − ηPk Pk⊤ gk , we obtain f (xk+1 ) ≤ f (xk ) − η gk⊤ Pk Pk⊤ gk +

η2 ⊤ g Pk Pk⊤ LPk Pk⊤ gk . 2 k

Taking conditional expectation with respect to Fk , and using E[Pk Pk⊤ | Fk ] = Id ,

E[Pk Pk⊤ LPk Pk⊤ | Fk ] ⪯ ℓLId ,

we obtain 2

E[f (xk+1 ) | Fk ] ≤ f (xk ) − η ∥gk ∥ +

ηℓL η 2 ℓL 2 2 ∥gk ∥ = f (xk ) − η 1 − ∥gk ∥ . 2 2

With η = 1/(ℓL), this yields E[f (xk+1 ) | Fk ] ≤ f (xk ) −

1 2 ∥gk ∥ . 2ℓL

Since f is differentiable and µ-strongly convex, 2

∥∇f (x)∥ ≥ 2µ f (x) − f ⋆

for all x ∈ Rd ,

we obtain E[f (xk+1 ) − f ⋆ | Fk ] ≤ f (xk ) − f ⋆ −

µ µ f (xk ) − f ⋆ = 1 − f (xk ) − f ⋆ . ℓL ℓL

Taking expectation proves (70), and (71) follows by iteration. The oracle bound (72) is immediate because each iteration uses r oracle calls. 40

High-probability and almost-sure guarantees

We now record two consequences of the supermartingale Lyapunov bounds: one in the convex case and one in the strongly convex case. In the convex case, the maximal inequality yields a uniform high-probability bound, and a dyadic Borel– Cantelli argument yields an almost sure eventual rate. In the strongly convex case, after normalization by the linear contraction factor, the Lyapunov process is again a nonnegative supermartingale, and the almost sure argument becomes simpler: no dyadic reduction is needed.

I.1

Convex case

Recall the convex Lyapunov process 1 ⋆ ΦC + ∥zk − x⋆ ∥2 . k := Ak f (xk ) − f 2 Proposition I.1 (Convex case: high-probability and almost sure convergence). Suppose Assumptions 2.1 to 2.2 hold, and let {xk , zk , Ak }k≥0 be generated by Algorithm 1. Then, by (36), the process (ΦC k )k≥0 is a nonnegative supermartingale: C E[ΦC for all k ≥ 0. k+1 | Fk ] ≤ Φk Consequently, the following hold. (i) For every η ∈ (0, 1), with probability at least 1 − η, ΦC Ak f (xk ) − f ⋆ ≤ 0 η

for all k ≥ 0.

(73)

Consequently, since Ak > 0 for all k ≥ 1, f (xk ) − f ⋆ ≤

ΦC 0 η Ak

for all k ≥ 1.

(74)

(ii) Fix any ε > 0. Then, almost surely, there exists a finite random integer Kε such that for all integers k ≥ Kε , l m 1+ε Ak f (xk ) − f ⋆ ≤ ΦC log log2 k . (75) 0 log2 k Consequently, for all k ≥ Kε , 1+ε ΦC 0 log2 k log⌈log2 k⌉ . f (xk ) − f ≤ Ak ⋆

(76)

Moreover, for the sequence (Ak ) generated by Algorithm 1, one has Ak ≥

m 2 k2 k = 2ω 4Lℓω

for all k ≥ 0,

(77)

1+ε log2 k log⌈log2 k⌉ . k2

(78)

and therefore, for all k ≥ Kε , f (xk ) − f

⋆

≤ 4Lℓω ΦC 0

Proof. We first prove part (i). Since (ΦC k ) is a nonnegative supermartingale, Ville’s maximal inequality gives, for every b > 0, ΦC C P sup Φj ≥ b ≤ 0 . b j≥0 41

Taking b = ΦC 0 /η, we obtain

ΦC 0 C P sup Φj ≤ ≥ 1 − η. η j≥0

Because Ak f (xk ) − f ⋆ ≤ ΦC k

a.s. for every k ≥ 0,

this proves (73). Since Ak > 0 for all k ≥ 1, (74) follows by dividing by Ak . We now prove part (ii). Fix ε > 0, and for every integer m ≥ 2, define Nm := 2m ,

1+ε bm := ΦC , 0 m(log m)

and the bad events

Em :=

maxm ΦC > b . m j

0≤j≤2

Applying Ville’s maximal inequality on the finite horizon {0, . . . , 2m } gives P(Em ) ≤ Since

1 ΦC 0 . = bm m(log m)1+ε

∞ X

1 < ∞, m(log m)1+ε m=2 the Borel–Cantelli lemma implies P(Em i.o.) = 0. Therefore, almost surely, there exists a finite random integer mε ≥ 2 such that for all m ≥ mε , C 1+ε max ΦC . j ≤ Φ0 m(log m)

0≤j≤2m

(79)

Now let Kε := 2mε . Fix any integer k ≥ Kε , and set m := ⌈log2 k⌉ . Then m ≥ mε and k ≤ 2m . Hence, by (79), C C 1+ε . Ak f (xk ) − f ⋆ ≤ ΦC k ≤ maxm Φj ≤ Φ0 m(log m) 0≤j≤2

Since m = ⌈log2 k⌉, this proves (75). Dividing by Ak gives (76). The lower bound (77) follows directly from (38). Finally, substituting (77) into (76) yields (78). Comparison with the expectation rate argument is

The expectation-level rate from the convex supermartingale

E[f (xk ) − f ⋆ ] ≤

ΦC 0 . Ak

By contrast, the uniform high-probability bound in Proposition I.1(i) differs only by the multiplicative factor 1/η: ΦC f (xk ) − f ⋆ ≤ 0 for all k ≥ 0 η Ak with probability at least 1 − η. The almost sure eventual bound in Proposition I.1(ii) differs from the expectation rate by the explicit logarithmic factor 1+ε log2 k log⌈log2 k⌉ . No monotonicity of f (xk ) is used anywhere in the proof. 42

I.2

Strongly convex case

Define the strongly convex Lyapunov process by ⋆ ΦSC k := f (xk ) − f +

µ ∥zk − x⋆ ∥2 . 2

Set ρ := 1 − θ. By Theorem 4.3, θ ∈ (0, 1], and hence ρ = 1 − θ ∈ [0, 1).Moreover, by (65), SC E[ΦSC k+1 | Fk ] ≤ ρ Φk

Also, by the definition of ΦSC k ,

f (xk ) − f ⋆ ≤ ΦSC k

for all k ≥ 0. for all k ≥ 0.

(80) (81)

Proposition I.2 (Strongly convex case: high-probability and almost sure convergence). Let ρ := 1 − θ. If ρ = 0, then ΦSC for all k ≥ 1 k =0 almost surely. Consequently, for every η ∈ (0, 1), with probability one, f (x0 ) − f ⋆ ≤

ΦSC 0 , η

f (xk ) − f ⋆ = 0

for all k ≥ 1,

and for every q ∈ (0, 1), (83) holds with Kq = 1. If ρ ∈ (0, 1), then the following hold. (i) For every η ∈ (0, 1), with probability at least 1 − η, f (xk ) − f ⋆ ≤

ΦSC 0 ρk η

for all k ≥ 0.

(82)

(ii) For every number q ∈ (ρ, 1), almost surely there exists a finite random integer Kq such that for all k ≥ Kq , k f (xk ) − f ⋆ ≤ ΦSC (83) 0 q . Proof. First consider the case ρ = 0. Since ΦSC k+1 ≥ 0 and (80) gives E[ΦSC k+1 | Fk ] ≤ 0, SC we have ΦSC k+1 = 0 almost surely for every k ≥ 0. Therefore Φk = 0 almost surely for all k ≥ 1. The claims in the case ρ = 0 follow from (81). It remains to consider the case ρ ∈ (0, 1). Define the normalized process

Mk := ρ−k ΦSC k . Then −k SC E[Mk+1 | Fk ] = ρ−(k+1) E[ΦSC Φk = Mk . k+1 | Fk ] ≤ ρ

Hence (Mk )k≥0 is a nonnegative supermartingale. For part (i), Ville’s maximal inequality gives ΦSC 0 P sup Mj ≥ ≤ η. η j≥0 Therefore, with probability at least 1 − η, Mk ≤

ΦSC 0 η

for all k ≥ 0, 43

that is, ΦSC k ≤

ΦSC 0 ρk η

for all k ≥ 0.

Using (81), we obtain (82). For part (ii), fix any q ∈ (ρ, 1), and define k ρ . ηk := q Since ρ/q < 1, the series gives

P∞

k=0 ηk converges.

Moreover, because E[Mk ] ≤ M0 = ΦSC 0 , Markov’s inequality

ΦSC 0 P Mk > ≤ ηk . ηk Equivalently, SC k P ΦSC ≤ k > Φ0 q Since

∞ k X ρ k=0

k ρ . q

< ∞,

the Borel–Cantelli lemma implies that, almost surely, only finitely many of these events occur. Hence, almost surely, there exists a finite random integer Kq such that for all k ≥ Kq , SC k ΦSC k ≤ Φ0 q .

Using (81), we obtain (83). Comparison with the expectation rate

The expectation-level strongly convex rate is

k E[f (xk ) − f ⋆ ] ≤ ΦSC 0 ρ .

The uniform high-probability bound in Proposition I.2(i) differs only by the multiplicative factor 1/η: f (xk ) − f ⋆ ≤

ΦSC 0 ρk η

for all k ≥ 0

with probability at least 1−η. The almost sure eventual bound in Proposition I.2(ii) replaces the contraction factor ρ by any prescribed factor q ∈ (ρ, 1). Thus the almost sure linear rate can be made arbitrarily close to the expectation linear rate, but it is not exactly identical.

Why the classical two-sequence template is not directly portable

In this appendix, we formalize a limitation of a direct transplantation of the classical two-sequence Nesterov template to randomized-subspace updates. These results should be interpreted as proof-template obstructions rather than impossibility results for all conceivable two-sequence accelerated schemes. Rather, they show that, under the sketch structure satisfied by the Haar and coordinate sketches, the classical proof mechanism is not directly compatible with the naive sketched analogue outside a near-full-dimensional regime. Proposition J.1. Let f : Rd → R be convex and L-smooth, and let x⋆ be a minimizer of f . Let {Pk }k≥0 be an i.i.d. sketch sequence satisfying E[Pk Pk⊤ ] = Id ,

Pk⊤ Pk =

d Ir r

for some 1 ≤ r ≤ d. Consider the direct sketched two-sequence recursion xk+1 = yk − ηPk Pk⊤ ∇f (yk ),

yk+1 = xk+1 + βk (xk+1 − xk ),

where η > 0 is constant and {βk }k≥0 is arbitrary. For each k ≥ 0, let Fk := σ(P0 , . . . , Pk−1 ),

Mk := Pk Pk⊤ .

gk := ∇f (yk ),

Suppose, as in the classical proof template, that there exist sequences {λk }k≥0 ⊂ (0, 1] and {zk }k≥0 ⊂ Rd such that, for every k ≥ 0, yk = (1 − λk )xk + λk zk ,

zk+1 = zk −

η Mk gk . λk

Then, for every k ≥ 0, 1 1 E 2 f (xk+1 ) − f ⋆ + ∥zk+1 − x⋆ ∥2 Fk λk 2η 1 − λk 1 η d ⋆ ⋆ 2 ≤ f (x ) − f + ∥z − x ∥ + (1 + Lη) − 2 ∥gk ∥2 . k k λ2k 2η 2λ2k r

(84)

(85)

In particular, if r ≤ d/2, then the coefficient η d (1 + Lη) − 2 2λ2k r is strictly positive for every η > 0. Therefore, the estimate (85) retains a strictly positive residual coefficient, and hence this affine hidden-sequence template does not recover the classical two-term decrease by the standard argument. Proof. Fix k ≥ 0. Since Pk⊤ Pk = (d/r)Ir , Mk2 = Pk (Pk⊤ Pk )Pk⊤ =

d Mk . r

Moreover, Pk is independent of Fk , and E[Mk ] = Id , so E[Mk2 | Fk ] =

d Id . r

L ∥v − u∥2 2

for all u, v ∈ Rd .

E[Mk | Fk ] = Id , By L-smoothness of f , f (v) ≤ f (u) + ⟨∇f (u), v − u⟩ +

Applying this with u = yk and v = xk+1 = yk − ηMk gk , we obtain f (xk+1 ) ≤ f (yk ) − ηgk⊤ Mk gk +

Lη 2 ⊤ 2 g M gk . 2 k k

Taking conditional expectation and using E[Mk | Fk ] = Id and E[Mk2 | Fk ] = dr Id gives E[f (xk+1 ) | Fk ] ≤ f (yk ) − η∥gk ∥2 +

Lη 2 d ∥gk ∥2 . 2 r

We now estimate the z-term. From (84), zk+1 − x⋆ = (zk − x⋆ ) − 45

η Mk gk , λk

(86)

and hence ∥zk+1 − x⋆ ∥2 = ∥zk − x⋆ ∥2 −

2η η2 ⟨zk − x⋆ , Mk gk ⟩ + 2 ∥Mk gk ∥2 . λk λk

Multiplying by 1/(2η), taking conditional expectation, and using again E[Mk | Fk ] = Id and E[Mk2 | Fk ] = d r Id , we get 1 E ∥zk+1 − x⋆ ∥2 Fk 2η 1 1 η d = ∥zk − x⋆ ∥2 − ⟨zk − x⋆ , gk ⟩ + 2 ∥gk ∥2 . (87) 2η λk 2λk r Subtracting f ⋆ from (86), multiplying by 1/λ2k , and adding (87), we obtain 1 1 ⋆ ⋆ 2 E 2 f (xk+1 ) − f + ∥zk+1 − x ∥ Fk λk 2η 1 1 1 ≤ 2 f (yk ) − f ⋆ + ∥zk − x⋆ ∥2 − ⟨zk − x⋆ , gk ⟩ λk 2η λk η d (1 + Lη) − 2 ∥gk ∥2 . + 2 2λk r

(88)

It remains to control the mixed term. Since yk = (1 − λk )xk + λk zk , we have

1 ⋆ 1 − λk 1 ⋆ (x − yk ) = (x − zk ). (xk − yk ) + λ2k λk λk

Taking inner products with gk gives −

1 1 − λk 1 ⟨zk − x⋆ , gk ⟩ = ⟨x⋆ − yk , gk ⟩ . ⟨xk − yk , gk ⟩ + 2 λk λk λk

(89)

By convexity of f , f ⋆ = f (x⋆ ) ≥ f (yk ) + ⟨gk , x⋆ − yk ⟩ ,

f (xk ) ≥ f (yk ) + ⟨gk , xk − yk ⟩ , so ⟨xk − yk , gk ⟩ ≤ f (xk ) − f (yk ),

⟨x⋆ − yk , gk ⟩ ≤ f ⋆ − f (yk ).

Substituting these bounds into (89) yields −

1 − λk 1 ⋆ 1 ⟨zk − x⋆ , gk ⟩ ≤ f (xk ) − f (yk ) + f − f (yk ) . λk λ2k λk

Substituting (90) into (88), we get 1 1 E 2 f (xk+1 ) − f ⋆ + ∥zk+1 − x⋆ ∥2 Fk λk 2η 1 1 − λk 1 ⋆ 1 ≤ 2 f (yk ) − f ⋆ + f (xk ) − f (yk ) + f − f (yk ) + ∥zk − x⋆ ∥2 2 λk λk λk 2η η d + 2 (1 + Lη) − 2 ∥gk ∥2 . 2λk r The f (yk )-terms cancel, since 1 1 − λk 1 − − = 0. λ2k λ2k λk 46

(90)

Therefore, 1 − λk 1 − λk 1 1 ⋆ f (yk ) − f ⋆ + f (xk ) − f (yk ) + f − f (yk ) = f (xk ) − f ⋆ , λ2k λ2k λk λ2k and (85) follows. Finally, if r ≤ d/2, then d/r ≥ 2, and thus d (1 + Lη) − 2 ≥ 2(1 + Lη) − 2 = 2Lη > 0 r

for every η > 0.

Hence the coefficient of ∥gk ∥2 in (85) is strictly positive. Therefore, (85) retains a strictly positive residual coefficient, and thus the estimate does not collapse to the classical two-term decrease by the standard argument. Proposition J.2. Let f : Rd → R be µ-strongly convex and L-smooth, with µ > 0. Let {Pk }k≥0 be an i.i.d. sketch sequence satisfying d E[Pk Pk⊤ ] = Id , Pk⊤ Pk = Ir r for some 1 ≤ r ≤ d. Consider the direct sketched two-sequence recursion xk+1 = yk − ηPk Pk⊤ ∇f (yk ),

yk+1 = xk+1 + β(xk+1 − xk ),

(91)

where η > 0 and β ∈ [0, 1) are constants. Define θ :=

1−β ∈ (0, 1], 1+β

γ :=

η , θ

zk :=

(1 + θ)yk − xk , θ

and Φk := f (xk ) − f ⋆ + Then yk =

1 θ xk + zk , 1+θ 1+θ

uk := (1 − θ)zk + θyk ,

µ ∥zk − x⋆ ∥2 . 2

zk+1 = uk − γPk Pk⊤ ∇f (yk ).

(92)

Writing gk := ∇f (yk ), Mk := Pk Pk⊤ , and Fk := σ(P0 , . . . , Pk−1 ), one has E[Φk+1 |Fk ] − (1 − θ)Φk 1−θ ⟨gk , xk − yk ⟩ ≤ θ f (yk ) − f ⋆ + (1 − θ) f (yk ) − f (xk ) − µγ⟨gk , yk − x⋆ ⟩ + µγ θ µθ µ(1 − θ) d + ∥yk − x⋆ ∥2 − ∥xk − yk ∥2 + −η + (Lη 2 + µγ 2 ) ∥gk ∥2 . 2 2θ 2r

(93)

If one wants the right-hand side of (93) to collapse into the same two brackets as in the standard stronglyconvex proof, namely E[Φk+1 |Fk ] − (1 − θ)Φk µ ≤ θ f (yk ) − f ⋆ − ⟨gk , yk − x⋆ ⟩ + ∥yk − x⋆ ∥2 2 + (1 − θ) f (yk ) − f (xk ) + ⟨gk , xk − yk ⟩ −

µ(1 − θ) ∥xk − yk ∥2 2θ

(94)

then matching the coefficients of the inner-product terms yields µγ = θ,

that is,

γ=

θ , µ 47

equivalently

η=

θ2 . µ

(95)

Under (95), (93) becomes E[Φk+1 |Fk ] − (1 − θ)Φk µ ≤ θ f (yk ) − f ⋆ − ⟨gk , yk − x⋆ ⟩ + ∥yk − x⋆ ∥2 2 + (1 − θ) f (yk ) − f (xk ) + ⟨gk , xk − yk ⟩ d Lθ2 θ2 µ(1 − θ) 2 ∥xk − yk ∥ + −1 + 1+ ∥gk ∥2 . − 2θ µ 2r µ

(96)

Therefore, in order for the remaining coefficient of ∥gk ∥2 in (96) to be nonpositive, it is necessary that Lθ2 d 1+ ≤ 2. (97) r µ Equivalently, since η = θ2 /µ, d (1 + Lη) ≤ 2. (98) r In particular, if r ≤ d/2, then (97) cannot hold for any β ∈ [0, 1). Therefore, outside a near-full-dimensional regime, the direct sketched analogue of the classical strongly-convex two-sequence Nesterov template is incompatible with the standard hidden-sequence proof mechanism. Proof. Fix k ≥ 0, and let Mk := Pk Pk⊤ ,

Fk := σ(P0 , . . . , Pk−1 ),

gk := ∇f (yk ).

We first rewrite (91) in hidden-sequence form. By the definition of zk , zk = so that

(1 + θ)yk − xk , θ

yk =

1 θ xk + zk . 1+θ 1+θ

zk+1 =

(1 + θ)yk+1 − xk+1 . θ

Also,

Using yk+1 = xk+1 +

1−θ (xk+1 − xk ), 1+θ

we obtain (1 + θ)yk+1 − xk+1 θ 1 1−θ = xk+1 − xk θ θ 1 1−θ xk . = (yk − ηMk gk ) − θ θ

zk+1 =

On the other hand, (1 − θ)zk + θyk = Hence

1−θ 1 1−θ (1 + θ)yk − xk + θyk = yk − xk . θ θ θ

η zk+1 = (1 − θ)zk + θyk − Mk gk = uk − γMk gk , θ

which proves (92). 48

We next derive the one-step Lyapunov estimate for a general γ, and only then identify the value of γ for which the same collapse pattern as in the standard proof occurs. From xk − yk = θ(yk − zk ), we obtain

1 zk − x⋆ = yk − x⋆ − (xk − yk ). θ

Also, since uk = (1 − θ)zk + θyk = yk −

1−θ (xk − yk ), θ

the update zk+1 = uk − γMk gk yields 1−θ (xk − yk ) − γMk gk . θ

zk+1 − x⋆ = yk − x⋆ − Therefore, a direct expansion gives ∥zk+1 − x⋆ ∥2 − (1 − θ)∥zk − x⋆ ∥2 = θ∥yk − x⋆ ∥2 −

2γ(1 − θ) 1−θ ∥xk − yk ∥2 − 2γ⟨Mk gk , yk − x⋆ ⟩ + ⟨Mk gk , xk − yk ⟩ + γ 2 gk⊤ Mk2 gk . (99) θ θ

Now Mk2 = Pk (Pk⊤ Pk )Pk⊤ =

d Mk . r

Since Pk is independent of Fk and E[Mk ] = Id , we have E[Mk2 | Fk ] =

E[Mk | Fk ] = Id ,

d Id . r

Taking conditional expectation in (99), we obtain hµ i (1 − θ)µ E ∥zk+1 − x⋆ ∥2 Fk − ∥zk − x⋆ ∥2 2 2 µθ 1−θ µ(1 − θ) d = ∥yk − x⋆ ∥2 − µγ⟨gk , yk − x⋆ ⟩ + µγ ⟨gk , xk − yk ⟩ − ∥xk − yk ∥2 + µγ 2 ∥gk ∥2 . (100) 2 θ 2θ 2r On the function-value side, by L-smoothness and xk+1 = yk − ηMk gk , we have f (xk+1 ) ≤ f (yk ) − ηgk⊤ Mk gk +

Lη 2 ⊤ 2 g M gk . 2 k k

Taking conditional expectation gives E[f (xk+1 ) | Fk ] ≤ f (yk ) − η∥gk ∥2 +

d Lη 2 ∥gk ∥2 . 2r

(101)

We now expand E[Φk+1 | Fk ] − (1 − θ)Φk : E[Φk+1 |Fk ] − (1 − θ)Φk = E[f (xk+1 ) | Fk ] − f ⋆ − (1 − θ)(f (xk ) − f ⋆ ) hµ i (1 − θ)µ + E ∥zk+1 − x⋆ ∥2 Fk − ∥zk − x⋆ ∥2 2 2 = E[f (xk+1 ) | Fk ] − f (yk ) + θ(f (yk ) − f ⋆ ) + (1 − θ)(f (yk ) − f (xk )) hµ i (1 − θ)µ + E ∥zk+1 − x⋆ ∥2 Fk − ∥zk − x⋆ ∥2 . 2 2 49

(102)

Substituting (100) and (101) into (102), we obtain E[Φk+1 |Fk ] − (1 − θ)Φk 1−θ ⟨gk , xk − yk ⟩ ≤ θ(f (yk ) − f ⋆ ) + (1 − θ)(f (yk ) − f (xk )) − µγ⟨gk , yk − x⋆ ⟩ + µγ θ µθ µ(1 − θ) d + ∥yk − x⋆ ∥2 − ∥xk − yk ∥2 + −η + (Lη 2 + µγ 2 ) ∥gk ∥2 , 2 2θ 2r which is (93). We now derive the value of γ for which the same collapse pattern as in the standard proof occurs. In (93), the inner-product terms are 1−θ ⟨gk , xk − yk ⟩. θ In order to rewrite the inner-product terms in (93) with the same coefficients as in (94), one must impose −µγ⟨gk , yk − x⋆ ⟩ + µγ

µγ = θ,

µγ

1−θ = 1 − θ. θ

γ=

θ . µ

η=

θ2 , µ

These two identities are equivalent, and therefore

Since γ = η/θ, this is equivalent to

which proves (95). Substituting γ = θ/µ and η = θ2 /µ into (93), we obtain E[Φk+1 |Fk ] − (1 − θ)Φk µ ≤ θ f (yk ) − f ⋆ − ⟨gk , yk − x⋆ ⟩ + ∥yk − x⋆ ∥2 2 + (1 − θ) f (yk ) − f (xk ) + ⟨gk , xk − yk ⟩ µ(1 − θ) d 2 2 2 − ∥xk − yk ∥ + −η + (Lη + µγ ) ∥gk ∥2 . 2θ 2r Since −η +

θ2 d d (Lη 2 + µγ 2 ) = − + 2r µ 2r

4 θ θ2 θ2 d Lθ2 L 2+ = −1 + 1+ , µ µ µ 2r µ

this yields (96). Finally, by µ-strong convexity, f (yk ) − f ⋆ − ⟨gk , yk − x⋆ ⟩ +

µ ∥yk − x⋆ ∥2 ≤ 0, 2

and by convexity, f (yk ) − f (xk ) + ⟨gk , xk − yk ⟩ ≤ 0. Thus, in order for the remaining coefficient of ∥gk ∥2 in (96) to be nonpositive, it is necessary that d Lθ2 1+ ≤ 2, r µ which is (97). The equivalent form (98) follows from η = θ2 /µ. Finally, if r ≤ d/2, then d/r ≥ 2, and hence for every θ ∈ (0, 1], d Lθ2 Lθ2 1+ ≥2 1+ > 2. r µ µ Therefore (97) cannot hold for any θ ∈ (0, 1], and hence for any β ∈ [0, 1). This completes the proof.

Record · ID 151726 · SHA-256 9c8e54a4c7b8882c
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.