Stochastic Bias-Reduced Gradient Methods
Year: 2021

Abstract:
We develop a new primitive for stochastic optimization: a low-bias, low-cost estimator of the minimizer x? of any Lipschitz strongly-convex function. In particular, we use a multilevel Monte-Carlo approach due to Blanchet and Glynn [8] to turn any optimal stochastic gradient method into an estimator of x? with bias , variance O(log(1/ )), and an expected sampling cost of O(log(1/ )) stochastic gradient evaluations. As an immediate consequence, we obtain cheap and nearly unbiased gradient estimators for the Moreau-Yoshida envelope of any Lipschitz convex function, allowing us to perform dimension-free randomized smoothing. We demonstrate the potential of our estimator through four applications. First, we develop a method for minimizing the maximum of N functions, improving on recent results and matching a lower bound up to logarithmic factors. Second and third, we recover state-of-the-art rates for projection-efficient and gradient-efficient optimization using simple algorithms with a transparent analysis. Finally, we show that an improved version of our estimator would yield a nearly linear-time, optimal-utility, differentially-private non-smooth stochastic optimization method.

heading: 1 Introduction
text: Consider the fundamental problem of minimizing a µ-strongly convex function F : X ! R given access to a stochastic (sub-)gradient estimator r̂F satisfying E r̂F (x) 2 @F (x) and Ekr̂F (x)k2  G
2 for every x 2 X . Is it possible to transform the unbiased estimator r̂F into a (nearly) unbiased estimator of the minimizer x? := argminx2X F (x)? In particular, can we improve upon the O(G/(µ p T )) bias achieved by T iterations of stochastic gradient descent (SGD)?
In this paper, we answer this question in the affirmative, proposing an optimum estimator x̂?, which (for any fixed > 0) has
bias kEx̂? x?k = O( ) and variance Ekx̂? Ex̂?k2 = O ✓ G 2
µ2 log
✓ G
µ
◆◆ ,
and, in expectation, costs O(log( Gµ )) evaluations of r̂F . 3 Setting = G/(µ
p T ), we obtain the
same bias bound as T iterations of SGD, but with expected cost of only O(log T ) stochastic gradient evaluations (the worst-case cost is T ). Further, the bias can be made arbitrarily small with only logarithmic increase in the variance and the stochastic gradient evaluations of our estimator, and therefore—paralleling the term “nearly linear-time” [27]—we call x̂? nearly unbiased.
Our estimator is an instance of the multilevel Monte Carlo technique for de-biasing estimator sequences [25] and more specifically the method of Blanchet and Glynn [8]. Our key observation is that this method is readily applicable to strongly-convex variants of SGD, or indeed any stochastic optimization method with the same (optimal) rate of convergence.
⇤Stanford University, {asi,jmblpati,yujiajin,sidford}@stanford.edu †Tel Aviv University, ycarmon@tauex.tau.ac.il 3When X = BR(x0) ⇢ Rd, F (x) = 1n P i2[n] F̂ (x; i), and r̂F is the subgradient of a uniformly random F̂ (x; i) we can also get an estimator with bias 0 and expected cost O(log(nd)). See Appendix A.1 for details.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).

heading: 1.1 Estimating proximal points and Moreau-Yoshida envelope gradients
text: Given a convex function f and regularization level , the proximal point of y is Pf, (y) := argminx2Rd f(x) + 2 kx yk 2
. Since computing Pf, amounts to solving a -strongly-convex problem, our technique provides low-bias and cheap proximal point estimators. Proximal points are ubiquitous in optimization [43, 19, 52, 38] and estimating them efficiently with low bias opens up new algorithmic possibilities. One of these possibilities is estimating the gradient of the Moreau-Yoshida envelope f (y) := minx2Rd f(x) + 2 kx yk 2
, which is a -smooth, G2/(2 )-accurate approximation of any G-Lipschitz f (see, e.g., [43, 29] and Appendix B.3). Since rf (y) = (y Pf, (y)), our optimum estimator provides a low-bias estimator for rf (y) with second moment and expected cost greater than those of r̂f by only a logarithmic factor. Thus, for any non-smooth f we can turn r̂f into a gradient estimator for the smooth surrogate f , whose smoothness is independent of the problem dimension, allowing us to perform dimension-free randomized smoothing [20].

heading: 1.2 Applications via accelerated bias-reduced methods
text: Our optimum estimator is a new primitive in stochastic convex optimization and we expect it to find multiple applications. We now describe three such applications: the first improves on previously known complexity bounds while the latter two recover existing bounds straightforwardly. For simplicity of presentation we assume (in the introduction only) Ekr̂fk2  1 and unit domain size. In each application, we wish to minimize an objective function given access to a cheap subgradient estimator r̂f as well as an expensive application-specific operation O (e.g., a projection to a complicated set). Direct use of the standard stochastic gradient method finds an ✏-accurate solution using O(✏ 2) computations of both r̂f and O, and our goal is to improve the O complexity without hurting the r̂f complexity.
To that end, we design stochastic accelerated methods consisting of T iterations, each one involving only a constant number of O and proximal point computations, which we approximate by averaging copies of our optimum estimator.4 Its low bias allows us to bound T ⌧ ✏ 2 as though our proximal points were exact, while maintaining an eO(✏ 2) bound on the total expected number of r̂f calls.5 Thus, we save expensive operations without substantially increasing the gradient estimation cost. Table 1 summarizes each application, and we briefly describe them below.
Minimizing the maximal loss (Section 4). Given N convex, 1-Lipschitz functions f(1), . . . , f(N) we would like to find an ✏-approximate minimizer of their maximum fmax(x) = maxi2[N ] f(i)(x). This problem naturally arises when optimizing worst-case behavior, as in maximum margin classification and robust optimization [53, 15, 45, 6]. We measure complexity by the number of individual function and subgradient evaluations, so that the expensive operation of evaluating f(1), . . . , f(N) at a single point has complexity O(N) and the subgradient method solves this problem with complexity O(N✏ 2). Carmon et al. [13] develop an algorithm for minimizing fmax with complexity eO(N✏ 2/3 + ✏ 8/3), improving on the subgradient method for sufficiently large N . Using our bias-reduced Moreau gradient envelope estimator in a Monteiro-Svaiter-type accelerated proximal point method [12, 11, 38], we obtain improved complexity eO(N✏ 2/3 + ✏ 2). This matches (up to logarithmic factors) a lower bound shown in [13], settling the complexity of minimizing the maximum
4While averaging is parallelizable, our optimum estimator itself is sequential. Consequently, our approach does not yield improve parallelism; see Appendix A.2 for further discussion
5It is easy to turn expected complexity bounds into deterministic ones; see Appendix A.3.
of N non-smooth functions. Our result reveals a surprising fact: for N ⌧ (GR/✏) 4/3, minimizing the maximum of N functions is no harder than minimizing their average.
Projection-efficient optimization via dimension-free randomized smoothing (Section 3). Consider the problem of minimizing a convex function f using an unbiased gradient estimator r̂f over convex set X for which Euclidean projections are expensive to compute (for example, the cone of PSD matrices). When f is L-smooth, a stochastic version of Nesterov’s accelerated gradient descent (AGD) [16] performs only O( p L/✏) projections. For non-smooth f we instead apply AGD to the Moreau envelope smoothing of f (with appropriate = O(✏ 1)) using our nearly-unbiased stochastic estimator for rf . This yields a solution in O(✏ 1) projections and eO(✏ 2) evaluations of r̂f . Our algorithm provides a simple alternative to the recent work of Thekumparampil et al. [51] whose performance guarantees are identical up to a logarithmic factor.
Gradient-efficient composite optimization (Section 5). We would like to minimize (x) = ⇤(x) + f(x), where ⇤ is convex and L-smooth but we can access it only via computing (expensive) exact gradients, while f is a non-smooth convex functions for which we have a (cheap) unbiased subgradient estimator r̂f . Problems of this type include inverse problems with sparsity constraints and regularized loss minimization in machine learning [34]. To save r⇤ computations, it is possible to use composite AGD [41] which solves O( p L/✏) subproblems of the form
minimizex hr⇤(y), xi+ f(x) + 2 kx x
0 k 2 . Lan [34] designed a specialized method, gradient sliding, for which the total subproblem solution cost is O(✏ 2) evaluations of r̂f . We show that a simple alternative—estimating the subproblem solutions via our low-bias optimum estimator—recovers its guarantees up to logarithmic factors.

heading: 1.3 Non-smooth differentially private stochastic convex optimization
text: We now discuss a potential application of our technique that is conditional on the existence of an improved optimum estimator. In it, we minimize the population objective function f(x) = ES⇠P f̂(x;S) under the well-known constraint of differential privacy [22]. Given n i.i.d. samples Si ⇠ P and assuming that each f̂ is 1-Lipschitz, convex and sufficiently smooth, Feldman et al. [23] develop algorithms that obtain the optimal error and compute O(n) subgradients of f̂ . The nonsmooth case is more challenging and the best existing bound is O(n11/8) for the high-dimensional setting d = n [32, 3]. In Section 6 we show that our optimum estimator, combined with recent localization techniques [23], reduces the problem to private mean estimation. Unfortunately, our estimator is heavy-tailed, leading to insufficient utility. Nevertheless, assuming a version of our estimator that has bounded outputs, we give an algorithm that queries eO(n) subgradients for nonsmooth functions, solving a longstanding open problem in private optimization [14, 4]. This motivates the study of improved versions of our estimators that have constant sensitivity.

heading: 1.4 Related work
text: Multilevel Monte-Carlo (MLMC) techniques originate from the literature on parametric integration for solving integral and differential equations [25]. Our approach is based on an MLMC variant put forth by Blanchet and Glynn [8] for estimating functionals of expectations. Among several applications, they propose [8, Section 5.2] an estimator for argminx ES⇠P f̂(x;S) where f̂(·; s) is convex for all s and assuming access to minimizers of empirical objectives of the form P i2[N ] f̂(x; si). The authors provide a preliminary analysis of the estimator’s variance (later elaborated in [9]) using an asymptotic Taylor expansion around the population minimizer. In comparison, we study the more general setting of stochastic gradient estimators and provide a complete algorithm based on SGD, along with a non-asymptotic analysis and concrete settings where our estimator is beneficial.
A number of works have used the Blanchet-Glynn estimator in the context of optimization and machine learning. These applications include estimating the ratio of expectations for semisupervised learning [7], estimating gradients of distributionally robust optimization objectives [35], and estimating gradients in deep latent variable models [47]. Our estimator is similar to that of Levy et al. [35] in that we also have to pick a “critical” doubling probability for the (random) computational budget, which makes the expected cost and variance of our estimators depend logarithmically on the bias.

heading: 1.5 Limitations
text: Our paper demonstrates that our proposed optimum estimator is a useful proof device: it allows us to easily prove upper bounds on the complexity of structured optimization problems, and at least in one case (minimizing the maximum loss) improve over previously known bounds. However, our work does not investigate the practicality of our optimum estimator, as implementation and experiments are outside its scope.
Nevertheless, let us briefly discuss the practical prospects of the algorithms we propose. On the one hand, our optimum estimator itself is fairly easy to implement, adding only a few parameters on top of a basic gradient method. On the other hand, in the settings of Sections 3 and 5, gradient-sliding based methods [34, 51] are roughly as simple to implement and enjoy slightly stronger convergence bounds (better by logarithmic factors) than our optimum estimator. Consequently, in these settings we have no reason to assume that our algorithms are better in practice. In the setting of Section 4 (minimizing the maximum loss) our algorithm does enjoy a stronger guarantee than the previous best bound [13]. However, both our algorithm and [13] are based on an accelerated proximal point method that, in its current form, is not practical [13, Sec. 6.2]. Thus, evaluating the benefit of stochastic bias reduction in the context of minimizing the maximum loss would require us to first develop a practical accelerated proximal point algorithm, which is an open question under active research [see, e.g., 50].
Another limitation of our optimum estimator is that, while it has a bounded second moment, its higher moments are unbounded. While this does not matter for most of our results, the lack of higher moment bounds prevents us from setting the complexity of non-smooth private stochastic convex optimization in Section 6. Finding an optimum estimator that is bounded with high probability—or proving that one does not exist—remains an open question for future work.
Finally, our analyses are limited to convex objective functions. However, while outside the scope of the paper, we believe our results are possibly relevant for non-convex settings as well. In particular, for smooth non-convex functions (and weakly-convex functions [17] more broadly) the problem of computing proximal points with sufficiently high regularization is strongly convex and our estimator applies. Such non-convex proximal points play an important role in non-convex optimization [17] with applications in deep learning [see, e.g., 49]. Applying the optimum-estimator technique in non-convex optimization is therefore a viable direction for future work.

heading: 1.6 Notation
text: We let BR(x) = {y 2 Rd : ky xk  R} denote the ball of radius R around x, where k·k is the Euclidean norm throughout. We write ProjS for the Euclidean projection to S . We write {A} for the indicator of event A, i.e., {A} = 1 when A holds and 0 otherwise. Throughout the paper, r̂f denotes a (stochastic) subgradient estimator for the function f , and X ⇢ Rd denotes the optimization domain, which we always assume is closed and convex. We use Pf, to denote the proximal operator (2) and f to denote the Moreau envelope (3) associated with function f and regularization parameter . Finally, we use Nf and Nr̂f to denote function and subgradient estimator evaluation complexity, respectively.

heading: 2 A multilevel Monte-Carlo optimum estimator
text: In this section, we construct a low-bias estimator for the minimizer of any strongly convex function F : X ! R. This estimator is the key component of our algorithms in the subsequent sections, which use it to approximate proximal points and Moreau envelope gradients. We assume that F is of the form F = f + , where the function is “simple” and that f satisfies the following. Assumption 1. The function f : X ! R is convex (with closed and convex domain X ) and is accessible via an unbiased subgradient estimator r̂f which satisfies Ekr̂f(x)k2  G2 for all x.
Our applications only use of the form (x) = 2 kx x 0 k 2 but our estimator applies more broadly to cases where argminx hv, xi+ (x) + 12⌘kx yk 2 is easy to compute for all v and y.

heading: 2.1 ODC algorithms
text: Our estimator can use, in a black-box fashion, any method for minimizing F with sufficiently fast convergence to x? = argminx2X F (x). We highlight the required convergence property as follows.
Algorithm 1: OPTEST(r̂f, , µ, , 2,X ) . ,
2 are required bias and square error . c is the ODC algorithm constant Tmax = l 4cG2
µ2 min{ 2, 12 2}
m
N = l 32cG2 log(Tmax)
µ2 2
m
for i = 1, . . . , N do x̂ (i) ? = a draw of the estimator (1) return 1N P i2[N ] x̂ (i) ?
Algorithm 2: MORGRADEST(r̂f, y, , , 2,X ) . ,
2 are required bias and square error . is the regularization level . y is the point at which to estimate rf (y) (x) =
2 kx yk
2
x̂? = OPTEST ⇣ r̂f, , , , 2 2 ,X ⌘
return (y x̂?)
Definition 1. An optimal-distance-convergence algorithm ODC takes as input r̂f satisfying Assumption 1, a simple function and a budget T 1. If F = f + is µ-strongly convex with minimizer x?, the algorithm’s output x = ODC(r̂f, , T ) requires at most T evaluations of r̂f to compute and satisfies Ekx x?k2  c G 2
µ2T for some constant c > 0.
Standard lower bound constructions imply that the O( G 2
µ2T ) squared distance convergence rate is indeed optimal; see Appendix A.4 for additional discussion. Conversely, ODC algorithms are readily available in the literature [44, 28] since any point x satisfying EF (x) F (x?) = O(G 2
µT ) (the optimal rate of convergence in strongly convex, Lipschitz optimization) also satisfies Ekx x?k2  O( G 2
µ2T ) by due to the strong convexity of F . We provide a concrete ODC algorithm consisting of a generalization of epoch SGD [28], which allows us to optimize over the composite objective F = f + instead of only f as in the prior study of epoch SGD. Lemma 1. EPOCHSGD (Algorithm 8 in Appendix B.1) is an ODC algorithm with constant c = 32.

heading: 2.2 Constructing an optimum estimator
text: To turn any ODC algorithm into a low-bias, low-cost and near-constant variance optimum estimator, we use the multilevel Monte Carlo (MLMC) technique of Blanchet and Glynn [8]. Given a problem instance r̂f, , an algorithm ODC and a cutoff parameter Tmax 2 N, our estimator x̂? is:
Draw J ⇠ Geom 1 2 2 N and, writing xj := ODC(r̂f, , 2j), compute
x̂? = x0 + ⇢ 2J(xJ xJ 1) 2J  Tmax 0 otherwise.
(1)
We note that for certain ODC algorithms it is possible to extract x0, xJ 1 from the intermediate steps of computing xJ , so that we only need to invoke ODC once. This is particularly simple to do for EPOCHSGD, as we explain in Appendix B.1. The key properties of our estimator are as follows.
Proposition 1. Let f and r̂f satisfy Assumption 1, F = f + be µ-strongly convex with minimizer x? and Tmax 2 N. For any ODC algorithm with constant c, the estimator (1) has bias kEx̂? x?k  p 2c G
µ p Tmax
and variance Ekx̂? Ex̂?k2  16cG 2
µ2 log2(Tmax). Moreover, the expected number of
r̂f evaluations required to compute x̂? is O(log Tmax).
Proof. Let jmax = max{j 2 N | 2j  Tmax} = blog2 Tmax)c. The expectation of x̂? is
Ex̂? = Ex0 + jmaxX
j=1
P(J = j)2j(Exj Exj 1) = Exjmax ,
where the second equality follows from P(J = j) = 2 j and the sum telescoping. Noting that xjmax = ODC(r̂f, , T ) for T = 2jmax Tmax/2, we have that
kExjmax x?k  q Ekxjmax x?k2  p c
G
µ p Tmax/2
by Definition 1. To bound the variance we use ka+ bk2  2kak2 + 2kbk2 and note that
Ekx̂? Ex̂?k2  Ekx̂? x?k2  2Ekx̂? x0k2 + 2Ekx0 x?k2.
The ODC property implies that Ekx0 x?k2  cG2/µ2. For the term Ekx̂? x0k2 we have
Ekx̂? x0k2 = jmaxX
j=1
P(J = j)22jEkxj xj 1k2 = jmaxX
j=1
2jEkxj xj 1k2, and
Ekxj xj 1k2  2Ekxj x?k2 + 2Ekxj 1 x?k2  6c G
2
µ2 2 j .
Substituting, we get Ekx̂? x0k2  6cG 2 µ2 jmax and Ekx̂? Ex̂?k 2  16cG 2
µ2 log2(Tmax). Finally, the expected number of r̂f evaluations is 1 + Pjmax j=1 P(J = j)(2j + 2j 1) = O(jmax).
The function OPTEST in Algorithm 1 computes an estimate of x? with and desired bias and square error 2 by averaging independent draws of the MLMC estimator (1). The following guarantees are immediate from Proposition 1; see Appendix B.2 for a short proof.
Theorem 1. Let f and r̂f satisfy Assumption 1, F = f + be µ-strongly convex with minimizer x? 2 X , and , > 0. The function OPTEST(r̂f, , µ, , 2,X ) outputs x̂? satisfying
kEx̂? x?k  and Ekx̂? x?k2  2
using N r̂f stochastic gradient computations, where
EN r̂f = O
✓ G 2
µ2 2 log2
✓ G
µmin{ , }
◆ + log ✓ G
µmin{ , }
◆◆ .

heading: 2.3 Estimating proximal points and Moreau envelope gradients
text: The proximal point of function f : X ! R with regularization level at point y is
Pf, (y) := argmin x2X f(x) + 2 kx yk
2 . (2)
When f satisfies Assumption 1, we may use OPTEST (with (x) = 2 kx yk 2 and µ = ) to obtain a reduced-bias proximal point estimator. The proximal point Pf, (y) is closely related to the Moreau envelope
f (y) := min x2X f(x) + 2 kx yk
2
(3)
via the relationship rf (y) = (y Pf, (y)) (see Appendix B.3). Therefore, we can use our optimum estimator to turn eO(1) calls to r̂f into a nearly unbiased estimator for rf . We formulate this as: Corollary 2. Let f and r̂f satisfy Assumption 1, let y 2 X and let , , > 0. The function MORGRADEST(r̂f, , y, , 2,X ) outputs r̂f (y) satisfying kEr̂f (y) rf (y)k  and Ekr̂f (y) rf (y)k2  2 and has complexity ENr̂f = O ⇣ G2 2 log 2 ⇣ G min{ , } ⌘ + log ⇣ G min{ , } ⌘⌘ .

heading: 3 Projection-efficient convex optimization
text: In this section, we combine the bias-reduced Moreau envelope gradient estimator with a standard accelerated gradient method to recover the result of Thekumparampil et al. [51]. We consider the problem of minimizing a function f satisfying Assumption 1 over the domain BR(0) subject to the constraint x 2 X , where X ⇢ BR(0) is a complicated convex set that we can only access via (expensive) projections of the form ProjX (x) := argminy2X ky xk. We further assume that an initial point x0 2 X satisfies kx0 x?k  D.
Algorithm 3 applies a variant of Nesterov’s accelerated gradient descent method (related to [2, 1]) on the ( -smooth) Moreau envelope f defined in eq. (3). Since computing the Moreau envelope
Algorithm 3: Stochastic accelerated gradient descent on the Moreau envelope Input: A gradient estimator r̂f satisfying Assumption 1 in BR(0), projection oracle ProjX , and initial point x0 = v0 with kx0 x?k  D. Parameters :Iteration budget T , Moreau regularization , approximation parameters k, 2k
1 for k = 1, · · · , T do 2 yk 1 = k 1 k+1xk + 2 k+1vk 1 3 gk = MORGRADEST(r̂f, yk 1, , k, 2k,BR(0)) 4 xk = ProjX yk 1 1 3 gk 5 vk = ProjBR(0) vk 1 k 6 gk
6 return xT
does not involve projection to X , for sufficiently accurate approximation of rf we require only T = O( p D2/✏) projections to X for finding an O(✏)-suboptimal point of f constrained to X . For that point to be also ✏-suboptimal for f itself, we must choose of the order of G2/✏, so that the number of projections is O(GD/✏).
As noted in [51] computing rf to accuracy O(✏/R) is sufficient for the above guarantee to hold, but doing so using a stochastic gradient method requires O((GD/✏)2) evaluations of r̂f per iteration, and O((GD/✏)3) evaluations in total. To improve this, we employ Algorithm 2 to compute nearlyunbiased estimates for rf and bound the error incurred by their variance. Our result matches the gradient sliding-based technique of Thekumparampil et al. [51] up to polylogarithmic factors while retaining the conceptual simplicity of directly applying AGD on the Moreau envelope. We formally state the guarantees of our method below, and provide a self-contained proof in Appendix C.
Theorem 3. Let f : BR(0) ! R and r̂f satisfy Assumption 1. Let X ✓ BR(0) be a convex set admitting a projection oracle ProjX . Let x0 2 X be an initial point with kx x?k  D for some x? 2 X . With = 2G2 ✏ , k = ✏ 8R , 2 k = 2✏ k+1 , and T = 7GD ✏ Algorithm 3 computes x 2 X with
E [f(x)]  f(x?)+ ✏ with complexity ENr̂f = O ⇣ G2D2 ✏2 log 2 GR ✏ ⌘ and O GD ✏ calls to ProjX .

heading: 4 Accelerated proximal methods and minimizing the maximal loss
text: In this section we apply our estimator in an accelerated proximal point method and use it to obtain an optimal rate for minimizing the maximum of N convex functions (up to logarithmic factors).

heading: 4.1 Accelerated proximal point method via Moreau gradient estimation
text: Algorithm 4 is an Monteiro-Svaiter-type [38, 12] accelerated proximal point method [36, 24] that leverages our reduced-bias Moreau envelope gradient estimator. To explain the method, we contrast it with stochastic AGD on the Moreau envelope (Algorithm 3). First and foremost, Algorithm 3 provides a suboptimality bound on the Moreau envelope f (which for small is far from f ) while Algorithm 4 minimizes f itself.
Second, while Algorithm 3 uses a fixed regularization parameter , Algorithm 4 handles an arbitrary sequence { k} given by a black-box function NEXTLAMBDA. To facilitate our application of the method to minimizing the maximal loss—where gradient estimation is only tractable in small Euclidean balls around a reference point—we include an optional parameter r such that the proximal point movement bound kPf, k+1(yk) ykk  r holds for all k. However, most of our analysis of Algorithm 4 does not require this parameter (i.e., holding for r = 1), making it potentially applicable to other settings that use accelerated proximal point methods [11, 38, 50].
The third and final notable difference between Algorithms 3 and 4 is the method of updating the xk iteration sequence. While a projected stochastic gradient descent step suffices for Algorithm 3, here we require a more direct approximation of function value decrease attained by the exact proximal mapping Pf, (see eq. (2)). For a given accuracy ', we define the '-approximate proximal mapping
eP'f, (y) := any x 2 X such that EF (x)  F (Pf, (y)) + ' for F (z) := f(z) +
2 kz yk
2 . (4)
Algorithm 4: Stochastic accelerated proximal point method Input: Gradient estimator r̂f , function NEXTLAMBDA, initialization x0 = v0 and A0 0. Parameters :Approximation parameters {'k, k, k}, stopping parameters Amax and Kmax,
optional movement bound r > 0. 1 for k = 0, 1, . . . do 2 k+1 = NEXTLAMBDA(xk, vk, Ak) . guaranteeing that kPf, k+1(yk) ykk  r 3 ak+1 =
1 2 k+1
p 1 + 4 k+1Ak and Ak+1 = Ak + ak+1 and Xk = X \ Br(yk)
4 yk = Ak
Ak+1 xk + ak+1 Ak+1 vk and xk+1 = eP'k+1f, k+1(yk) . defined in eq. (4) 5 gk+1 = MORGRADEST(r̂f, yk, k+1, k, 2k,Xk) and vk+1 = ProjX vk 1 2ak+1gk+1 6 if Ak+1 Amax or k + 1 = Kmax then return xk+1
Note that eP0f, = Pf, and that for ' > 0 we can compute eP ' f, with an appropriate SGD variant (such as EPOCHSGD) using O(G2/( ')) evaluations of r̂f .
With the differences between the algorithms explained, we emphasize their key similarity: both algorithms update the vk sequence using our bias reduction method MORGRADEST (Algorithm 2), which holds the key to their efficiency. The following proposition shows that Algorithm 4 has the same bound on Kmax as an exact accelerated proximal point method [12], while requiring at most eO(G2R2✏ 2) stochastic gradient evaluations; see proof in Appendix D.1.
Proposition 2. Let f : X ! R and r̂f satisfy Assumption 1, and let X ✓ BR(x0). For a target accuracy ✏  GR let 'k+1 =
✏ 60 k+1ak+1 , k+1 = ✏ 120R , 2 k+1 = ✏ 60ak+1 , A0 = R G , and
Amax = 9R2 ✏ . If k min 1 Amax = ⌦( ✏R2 ) for all k  Kmax, then lines 4 and 5 of Algorithm 4
have total complexity EN r̂f = O ⇣ Kmax log GR ✏ + G2R2 ✏2 log 2 GR ✏ ⌘ . If in addition kPf, k(yk 1) yk 1k 3r/4 whenever k 2 min then for Kmax = O ⇣ R r 2/3 log GR ✏ + q minR2 ✏ ⌘ , the algorithm’s output xK satisfies f(xK) f(x?)  ✏ with probability at least 2 3 .

heading: 4.2 Minimizing the maximal loss
text: We now consider objectives of the form fmax(x) := maxi2[N ] f(i)(x) where each function f(i) : X ! R is convex and G-Lipschitz. Our approach to minimizing fmax largely follows Carmon et al. [13]; the main difference is that we approximate proximal steps via Algorithm 4 and our reduced-bias bias estimator. The first step of the approach is to replace fmax with the “softmax” function, defined for a given target accuracy ✏ as
fsmax(x) := ✏ 0 log
X
iN
exp f(i)(x)/✏
0 ! , where ✏0 :=
✏
2 logN .
Since fsmax(x) fmax(x) 2 [0, ✏2 ], any ✏ 2 -accurate solution of fsmax is ✏-accurate for fmax.
The second step is to develop an efficient gradient estimator for fsmax; this is non-trivial because fsmax is not a finite sum or expectation. In [13] this is addressed via an “exponentiated softmax” trick; we develop an alternative, rejection sampling-based approach that fits Algorithm 4 more directly (see Algorithm 9). To produce an unbiased estimate for rfsmax(x) for x in a ball of radius r✏ = ✏0/G we require a single rf(i)(x) evaluation (for some i), O(1) evaluations of f(i)(x) in expectation, and evaluation of the N functions f(1)(y), . . . , f(N)(y) for pre-processing. Plugging this estimator into Algorithm 4 with r = r✏, the total pre-processing overhead of lines 4 and 5 is O(KmaxN).
The final step is to find a function NEXTLAMBDA such that kPfsmax, t+1(yk) ykk  r✏ for all k (enabling gradient estimation), and kPfsmax, t+1(yk) ykk 3 4r✏ when k+1 > 2 min (allowing us to bound Kmax with Proposition 2). Here we use the bisection subroutine from [13] as is (see Algorithm 10). By judiciously choosing min—an improvement over the analysis in [13]—we obtain the following complexity guarantee on Nf(i) and N@f(i) , the total numbers of individual function and subgradient evaluations, respectively. (See proof Appendix D.2).
Theorem 4. Let f(1), . . . , f(N) : X ! R be convex and G-Lipschitz and let X ✓ BR(x0). For any ✏ <
1 2GR/ logN , Algorithm 4 (with eP'fsmax, implemented in Algorithm 8, r̂fsmax given by Algorithm 9 and NEXTLAMBDA given by Algorithm 10 with min = eO(✏/(r4/3✏ R2/3)) outputs x 2 X that with probability at least
1 2 is ✏-suboptimal for fmax(x) = maxi2[N ] f(i)(x) and has complexity
ENf(i) = O ✓ N ⇣ GR logN ✏ ⌘2/3 + GR ✏ 2 log2 GR✏ ◆ and EN@f(i) = O ⇣ GR ✏ 2 log2 GR✏ ⌘ .
The rate given by Theorem 4 matches (up to logarithmic factors) the lower bound ⌦(N(GR/✏)2 + (GR/✏)2) shown in [13] and is therefore near-optimal.

heading: 5 Gradient-efficient composite optimization
text: Consider the problem of finding a minimizer of the following convex composite optimization problem
minimize x2X (x) := ⇤(x) + f(x) where ⇤ is L-smooth and f satisfies Assumption 1, (5)
given x0 such that kx0 x?k  R for some x? 2 argminx2X (x). Lan [34] developed a method called “gradient sliding” that finds an ✏-accurate solution to (5) with complexity Nr⇤ = O( p LR2/✏)
evaluations of r⇤(x) and N r̂f = O((GR/✏) 2) evaluations of r̂f(x), which are optimal even for each component separately.6
In this section, we provide an alternative algorithm that matches the complexity of gradient up to logarithmic factors and is conceptually simple. Our approach, Algorithm 5, is essentially composite AGD [41], where at the kth iteration we compute a proximal point (2) with respect to a partial linearization of around yk. In particular, letting ⇤̄k(v) := ⇤(yk)+ hr⇤(yk), v yki and k = 2Lk , we approximate P⇤̄+f, k(vk 1). Similar to Algorithm 4, Algorithm 5 computes two types of approximations: one is an ✏k-approximate proximal point eP✏k⇤̄+f, k(vk 1) as per its definition (4), while the other is our bias-reduced optimum estimator from Algorithm 1. We note, however, that unlike Algorithm 4 which approximates the xk update, here we approximate vk, the “mirror descent” update.
Below we state the formal guarantees for Algorithm 5; we defer its proof to Appendix E.
Algorithm 5: Stochastic composite accelerated gradient descent Input: A problem of the form (5) with ⇤, f , r⇤, r̂f . Parameters :Step size parameters k = 2Lk and k = 2 k+1 , iteration number N , approximation
parameters {✏k, k, 2k} and x0 = v0 satisfying kx0 x?k  R. 1 for k = 1, 2, · · · , N do 2 yk = (1 k)xk 1 + kProjX (vk 1) 3 v̄k = eP✏k⇤̄k+f, k(vk 1) for ⇤̄k(v) := ⇤(yk) + hr⇤(yk), v yki 4 vk = OPTEST(r̂f, k, k, k, 2k,BR(v0) \ X ) for k(z) = k 2 kz vk 1k
2 + ⇤̄k(z) 5 xk = (1 k)xk 1 + kv̄k
6 return xN
Theorem 5. Given problem (5) with solution x?, a point x0 such that kx0 x?k  R and target accuracy ✏ > 0, Algorithm 5 with ✏k = LR/2kN , k = R/16N , 2k = R 2 /4N , and N =
⇥( p LR2/✏) finds an approximate solution x satisfying E (x)  (x?) + ✏ and has complexity
Nr⇤ = O
✓q LR2
✏
◆ and EN
r̂f = O ✓ GR ✏ 2 log2 GR✏ + q LR2 ✏ log GR ✏ ◆ .

heading: 6 Efficient non-smooth private convex optimization
text: We conclude the paper with a potential application of our optimum estimator for differentially private stochastic convex optimization (DP-SCO). In this problem we are given n i.i.d. sample
6The gradient sliding result holds under a relaxed Lipschitz assumption [see 34, eq. (1.2)]. It is straightforward to extend EPOCHSGD, and hence all of our results, to that assumption as well.
si ⇠ P taking values in a set S, and out objective is to privately minimize the population average f(x) = ES⇠P [f̂(x;S)], where f̂ : X ⇥ S ! R, is convex in the first argument and X ⇢ Rd is a convex set, and S is a population of data points. That is, We wish to find x̂ 2 X with small excess loss f(x̂) minx2X f(x) while preserving differential privacy. Definition 2 ([22]). A randomized algorithm A is (↵, )-differentially private ((↵, )-DP) if, for all datasets S,S 0 2 Sn that differ in a single data element and for every event O in the output space of A, we have P [A(S) 2 O]  e↵P [A(S 0) 2 O] + .
DP-SCO has received increased attention over the past few years. Bassily et al. [5] developed (inefficient) algorithms that attain the optimal excess loss 1/ p n + p d log(1/ )/n↵. When each function is O( p n) smooth, Feldman et al. [23] gave algorithms with optimal excess loss and O(n) gradient query complexity. In the non-smooth setting, however, their algorithms require O(n2) subgradients. Subsequently, Asi et al. [3] and Kulkarni et al. [32] developed more efficient algorithms for non-smooth functions which need O(min(n2/ p d, n 5/4 d 1/8 , n 3/2 /d 1/8)) subgradients which is O(n11/8) for the high-dimensional setting d = n. Whether a linear gradient complexity is achievable for DP-SCO in the non-smooth setting is still open.
In this section, we develop an efficient algorithm for non-smooth DP-SCO that queries eO(n) subgradients conditional on the existence of an optimum estimator with the following properties. Definition 3. Let F = f + be µ-strongly convex with minimizer x? and f is G-Lipschitz. For > 0, we say that O is efficient bounded low-bias estimator if it returns x̂? = O (F ) such that kE[x̂? x?]k2  2, kx̂? x?k2  C1G2 log(G/µ )/µ2, and the expected number of gradient queries is C2 log(G/µ ).
Comparing to our MLMC estimator (1) and Proposition 1, we note that the only place our current estimator falls short of satisfying Definition 3 is the probability 1 bound on kx̂? x?k2, which for (1) holds only in expectation. Indeed, for our estimator, kx̂? x?k can be as large as O(G/(µ )), meaning that it is heavy-tailed.
It is not clear whether an EBBOE as defined above exists. Nevertheless, assuming access to such estimator, Algorithm 6 solves the DP-SCO problem with a near-linear amount of gradient computations. The algorithm builds on the recent localization-based optimization methods in [23] which iteratively solve regularized minimization problems.
Algorithm 6: Differentially-private stochastic convex optimization via optimum estimation Input: (s1, . . . , sn) 2 Sn, domain X ⇢ BR(x0), EBBOE O (satisfying Definition 3).
1 Set k = dlog ne, B = 20(log( 1 ) + C2 log 2 n), n̄ = nk , ⌘ = R G min
( 1
p n ,
↵
B log(n) q d log( 1 )
)
2 for i = 1, 2, · · · , k do 3 Let ⌘i = 2 4i⌘ , fi(x) = 1n̄ Pkn̄ j=1+(k 1)n̄ f̂(x; sj), i(x) = kx xi 1k 2 /(⌘in̄) 4 Let x̃i = 1n̄ Pn̄ j=1 O i(Fi) with Fi = fi + i , 2 i = G 2 ⌘ 2 i n̄ 5 Set xi = x̃i + ⇣i where ⇣i ⇠ N(0, 2i Id) with i = 8B( p C1 log n+ 2)⌘i p log(2/ )/↵i
6 return xk
We average multiple draws of the (hypothetical) bounded optimum estimator to solve the regularized problems, and apply private mean estimation procedures to preserve privacy. We defer the proof of the following results Appendix F. Theorem 6 (conditional). Given an efficient bounded low-bias estimator O satisfying Definition 3 for any > 0, then for ↵  log(1/ ), X 2 BR(x0), convex and G-Lipschitz f̂(x; s), Algorithm 6 is (↵, )-DP, queries eO(n) subgradients and has (hiding logarithmic factors in n) E[f(xk)
minx2X f(x)]  GR · eO ✓
1 p n +
p d log3(1/ )
n↵
◆ .
Theorem 6 provides a strong motivation for constructing bounded optimum estimators that satisfy Definition 3 . In Appendix F.3, we discuss the challenges in making our MLMC estimator bounded, as well as some directions to overcome them.

heading: Acknowledgments
text: HA was supported by ONR YIP N00014-19-2288 and the Stanford DAWN Consortium. YC was supported in part by Len Blavatnik and the Blavatnik Family foundation, and the Yandex Machine Learning Initiative for Machine Learning. YJ was supported by a Stanford Graduate Fellowship. AS was supported in part by a Microsoft Research Faculty Fellowship, NSF CAREER Award CCF-1844855, NSF Grant CCF-1955039, a PayPal research award, and a Sloan Research Fellowship.

Review summaries:
The paper develops a new primitive for stochastic optimization – a low-biased and low-cost estimator the minimizer x* for any Lipschitz strongly-convex function. Specially, the author(s) use a multi-level Monte-Carlo method to turn any optimal stochastic gradient method into an estimator x*. The paper also shows the potential of the estimator through four applications. 

This paper proposes a procedure that converts any optimal stochastic gradient method into an estimator of the optimum $x^*$ using the multilevel Monte-Carlo approach. The estimator has a smaller expected cost of stochastic gradient computation compared with the original stochastic gradient method. The authors demonstrate that the proposed estimator can be combined with existing methods to improve the computation complexity and provide alternative algorithms to other methods in different applications.  

The paper proposes a constructive way to obtain an estimator of the minimizer of a Lipschitz strongly convex function $f$. Assuming that we have an access to a stochastic subgradient of $f$, we can run a variant of SGD that guarantees good complexity. Furthermore, the authors demonstrate on several examples how this estimator leads to new algorithms with a better complexity.   I think it is a good paper.