A Unified Theory of Stochastic Proximal Point Methods without Smoothness

## Abstract

Abstract This paper presents a comprehensive analysis of a broad range of variations of the stochastic proximal point method ( SPPM ). Proximal point methods have attracted considerable interest owing to their numerical stability and robustness against imperfect tuning, a trait not shared by the dominant stochastic gradient descent ( SGD ) algorithm. A framework of assumptions that we introduce encompasses methods employing techniques such as variance reduction and arbitrary sampling. A cornerstone of our general theoretical approach is a parametric assumption on the iterates, correction and control vectors. We establish a single theorem that ensures linear convergence under this assumption and the μ 𝜇 \mu italic_μ -strong convexity of the loss function, and without the need to invoke smoothness. This integral theorem reinstates best known complexity and convergence guarantees for several existing methods which demonstrates the robustness of our approach. We expand our study by developing three new variants of SPPM , and through numerical experiments we elucidate various properties inherent to them.

## 1 Introduction

In this paper we consider the stochastic optimization problem

$$ $\min_{x\in\mathbb{R}^{d}}\left\{f(x):={\rm E}_{\xi\sim\mathcal{D}}\left[f_{\xi }(x)\right]\right\},$ (1) $$

where $\xi\sim\mathcal{D}$ is a random variable following distribution $\mathcal{D}$, and ${\rm E}\left[\cdot\right]$ denotes mathematical expectation. Problems of this type are fundamental for statistical supervised learning theory. Here, $x$ is a machine learning model of $d\in\mathbb{N}$ parameters, $\mathcal{D}$ is an unknown distribution of labeled examples, samples $\xi\sim\mathcal{D}$ are available, $f_{\xi}$ is the loss on datapoint $\xi,$ $f$ is the generalization error. In such a setup, an unbiased estimator of the gradient $\nabla f_{\xi}(x)$ is computable, while the gradient $\nabla f(x)$ itself is not. We rely on two assumptions, presented next.

###### Assumption 1 .

Function $f_{\xi}:\mathbb{R}^{d}\to\mathbb{R}$ is differentiable for all samples $\xi\sim\mathcal{D}$.

We implicitly assume that the order of differentiation and expectation can be swapped, which means that $\nabla f(x)\overset{\eqref{eq:main-opt-problem}}{=}\nabla{\rm E}_{\xi\sim
\mathcal{D}}\left[f_{\xi}(x)\right]={\rm E}_{\xi\sim\mathcal{D}}\left[\nabla f
_{\xi}(x)\right].$
This implies that $f$ is differentiable.

###### Assumption 2 .

Function $f_{\xi}:\mathbb{R}^{d}\to\mathbb{R}$ is $\mu$-strongly convex for all samples $\xi\sim\mathcal{D}$, where $\mu>0:$

$$ $f_{\xi}(y)+\left\langle\nabla f_{\xi}(y),x-y\right\rangle+\frac{\mu}{2}{\left \|x-y\right\|^{2}}\leq f_{\xi}(x),\quad\forall x,y\in\mathbb{R}^{d}.$ (2) $$

This implies that $f$ is $\mu$-strongly convex, and hence $f$ has a unique minimizer, which we denote by $x_{\star}$. We know that $\nabla f(x_{\star})=0$. Notably, we do not assume $f$ to be $L$-smooth.

Another type of problems considered in the paper is a minimization of functions $f$ that are averages of a large number of differentiable functions:

$$ $\min_{x\in\mathbb{R}^{d}}\left\{f(x):=\frac{1}{n}\sum_{i=1}^{n}f_{i}(x)\right\}.$ (3) $$

This is a special case of (1), with $\mathcal{D}$ being the uniform distribution over the finite set $[n]$: $\xi=i$ with probability $\frac{1}{n}.$ Problems with this structure commonly emerge in practice during the training of supervised machine learning models via empirical risk minimization. They are known as finite-sum optimization problems. We rely on two assumptions, presented next.

###### Assumption 3 .

Function $f_{i}:\mathbb{R}^{d}\to\mathbb{R}$ is differentiable for all $i\in[n]$.

This implies that $f$ is differentiable.

###### Assumption 4 .

Function $f_{i}:\mathbb{R}^{d}\to\mathbb{R}$ is $\mu_{i}$-strongly convex for all $i\in[n]$, where $\mu_{i}>0$. That is,

$$ $f_{i}(y)+\left\langle\nabla f_{i}(y),x-y\right\rangle+\frac{\mu_{i}}{2}{\left \|x-y\right\|^{2}}\leq f_{i}(x),\quad x,y\in\mathbb{R}^{d}.$ (4) $$

This refines 2, with the the same strong convexity parameter for all functions. This implies that $f$ is $\mu$-strongly convex with $\mu=\min_{i}\mu_{i}$. Hence, $f$ has a unique minimizer $x_{\star}.$

###### Assumption 1 .

###### Assumption 2 .

###### Assumption 3 .

###### Assumption 4 .

## 2 Variations of SPPM

The most common algorithm for finding an $\varepsilon$-accurate solution of problems (1) and (3) is stochastic gradient descent (SGD)Robbins and Monro (); Nemirovski et al. (); Gower et al. (). In general, SGD updates have the form of $x^{k+1}=x^{k}-\gamma g^{k},$ where $g^{k}$ is an unbiased stochastic gradient estimator: ${\rm E}\left[\left.g^{k}\right|x^{k}\right]=\nabla f(x^{k}).$ There are numerous approaches to constructing an estimator with the aim of attaining preferable algorithmic features such as faster convergence rate, lower iteration cost, parallelization, and generalization. One of the most significant issues in SGD variations is the difficulties with a suitable selection of the stepsize. Theoretical results highlight that the stepsize is restricted to small values (Bach and Moulines, ).

Whenever we are able to evaluate the stochastic proximity operator, another option to consider in place of SGD is the stochastic proximal point method (SPPM), whose iterations have the form

$$ ${\rm prox}_{\gamma f_{\xi_{k}}}\left(x_{k}\right):=\arg\min_{x\in\mathbb{R}^{d }}\left\{f_{\xi_{k}}(x)+\frac{1}{2\gamma}{\left\|x-x_{k}\right\|^{2}}\right\}, \quad\xi_{k}\sim\mathcal{D},$ $$

where $\gamma>0$ is a stepsize. Clearly, the proximity operator is well-defined due to the strong convexity 2 on $f_{\xi}(x).$ Given the error tolerance $\varepsilon>0,$ we derive in Section B.2 (see (23); also, see (Asi and Duchi, , Proposition 5.3), (Khaled and Jin, , Theorem 1)) that with $\sigma_{\star}^{2}:={\rm E}_{\xi\sim\mathcal{D}}\left[{\left\|\nabla f_{\xi}(x
_{\star})\right\|^{2}}\right],$ the stepsize $\gamma=\frac{\mu\varepsilon}{\sigma_{\star}^{2}}$, we get ${\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\varepsilon$
provided that

$$ $k\geq\left(\frac{1}{2}+\frac{\sigma_{\star}^{2}}{2\mu^{2}\varepsilon}\right) \log\left(\frac{2{\left\|x_{0}-x_{\star}\right\|^{2}}}{\varepsilon}\right).$ $$

As shown by Gower et al. (), in the same setting, SGD with fixed stepsize reaches an $\varepsilon$-accurate solution after

$$ $k\geq\left(\frac{2L}{\mu}+\frac{2\sigma_{\star}^{2}}{\mu^{2}\varepsilon}\right )\log\left(\frac{4{\left\|x_{0}-x_{\star}\right\|^{2}}}{\varepsilon}\right)$ $$

iterations ($L$ is a bound on the smoothness constant of stochastic functions). Note that although both iteration complexities are dependent on the stochastic noise term, the iteration complexity of SGD additionally hinges on the condition number. $\kappa:=\frac{L}{\mu}\geq 1$ ($\mu$ is a strong convexity parameter of $f$ here). In contrast, the iteration complexity of SPPM remains unaffected by the smoothness constant $L$. Consequently, if we have access to stochastic proximal operator evaluations, we can achieve a faster convergence rate than SGD. Another important aspect is that SPPM still works for large stepsizes. The primary distinction with the result of Khaled and Jin () is that the neighborhood guaranteed in our analysis for SPPM does not blow up to infinity as the stepsize $\gamma$ grows to infinity (see (21) in Theorem 2 and Commentary 2 after it).

Ryu and Boyd (); Asi and Duchi () demonstrated that SPPM exhibits greater resilience in terms of the choice of stepsize compared to SGD. Ryu and Boyd () furnish convergence rates for SPPM and note its stability against learning rate misspecification, a trait not shared by SGD. Asi and Duchi () examine a broader method (AProx) which encompasses SPPM as a particular instance, providing both stability and convergence rates under convexity. Additionally, the convergence rates of SPPM remain consistent with those of SGD across different versions of the algorithms.

We have discussed the versatility in designing the unbiased stochastic estimator $g^{k}$ in SGD, which can be accomplished in various manners. Among these, popular sampling strategies include importance sampling and mini-batching. A comprehensive analysis of these methods within the framework of arbitrary sampling was presented by Gower et al. (). Similarly, an analogous effort is made for SPPM in our paper, where sampling schemes and methods such as SPPM-NS and SPPM-AS are proposed (see Appendix B.3 and B.4). While sampling strategies offer significant utility, both for SGD and SPPM variations, they tend to converge to the vicinity of the solution for fixed stepsizes, reaching the exact solution only in the overparameterized regime.

That being said, the issue of SGD iterates not converging to the optimum prompted the development of variance-reduced methods. These methods enhance the convergence rate of SGD for finite-sum problems by constructing gradient estimators with diminishing variance over time (e.g., SAGA (Defazio et al., ) and SVRG (Johnson and Zhang, )). In other words, variance-reduced methods progressively acquire knowledge of the stochastic gradients at the optimum and mitigate the influence of gradient noise. Consequently, for strongly convex $f$, these variance-reduced methods exhibit linear convergence towards $x_{\star}$ with a fixed step size (see a survey by Gower et al. ()). In the discussion above we concluded that SPPM has advantages over SGD. This motivated Khaled and Jin () to explore variance-reduced variants of SPPM. They proposed the SVRP algorithm with an approximate calculation of the proximal operator and demonstrated that variance reduced variants of SPPM have better convergence guarantees under second-order similarity assumption for the finite-sum setting (3). We present and employ in this paper a more general form of this assumption for the stochastic setting (1) in order to analyze L-SVRP with the proximal operator in the updates. We stipulate the similarity assumption to be valid exclusively at $x_{\star},$ whereas Khaled and Jin () require it to hold for any $y\in\mathbb{R}^{d}$ instead.

###### Assumption 5 (Similarity) .

There exists $\delta\geq 0$ such that

$$ ${\rm E}_{\xi\sim\mathcal{D}}\left[{\left\|\nabla f_{\xi}(x)-\nabla f(x)-\nabla f _{\xi}(x_{\star})\right\|^{2}}\right]\leq\delta^{2}{\left\|x-x_{\star}\right\| ^{2}},\qquad\forall x\in\mathbb{R}^{d}.$ (5) $$

Let us also note here that the standard $\delta$-smoothness assumption implies the second-order similarity assumption by Khaled and Jin (), which in turn implies 5. In Section B.6 we discuss the generality of 5 in detail. In the studies of Szlendak et al. (); Panferov et al. (), federated optimization with compression under the similarity assumption is explored, leading to improved convergence rates achieved through specially designed compression operators.

Point SAGA, also a variant of SPPM, was proposed by Defazio (), where its convergence was analyzed under the individual smoothness assumption. The algorithm requires a large amount of memory for its execution. It inherits this shortcoming from its SGD variance-reduced archetype SAGA by Defazio et al. (). The structure of Point SAGA will not allow us to perform the analysis under the similarity assumption (5). Instead, we will rely on this stronger assumption.

###### Assumption 6 .

We assume that there exists $\nu>0$
such that the inequality

$$ $\frac{1}{n}\sum^{n}_{j=1}{\left\|\nabla f_{j}(x^{j})-\frac{1}{n}\sum^{n}_{i=1} \nabla f_{i}(x^{i})-\nabla f_{j}(x_{\star})\right\|^{2}}\leq\nu^{2}\frac{1}{n} \sum^{n}_{j=1}{\left\|x^{j}-x_{\star}\right\|^{2}}$ (6) $$

holds for all $x^{1},\dots,x^{n}\in\mathbb{R}^{d}$.

Let us note that this condition is weaker than the individual $\nu$-smoothness assumption. In Section B.8 we discuss the generality of 6 in more detail.

Traoré et al. () analyze L-SVRP and Point SAGA algorithms in the setting (3) when each $f_{i}$ is smooth, convex, and either $f$ is convex or satisfies the PŁ-condition. In contrast, in our work we consider a different setting with individually strongly convex functions under 5 or 6 (both are weaker than the individual smoothness assumption). The convergence theory of variance-reduced SPPM methods significantly differs from that of standard SPPM. We suggest the possibility of a unified theory that encompasses both SPPM and its variance-reduced counterparts.

###### Assumption 5 (Similarity) .

###### Assumption 6 .

## 3 Contributions

Numerous efficient adaptations of the SPPM algorithm have emerged, each with its specific applications. Our research stems from the absence of a comprehensive and universally applicable theory. While some connections among existing methods have been established (as demonstrated by Traoré et al. (), who link the classical SPPM method with its variance-reduced counterparts), a cohesive theoretical framework without the smoothness assumption is missing. Understanding the connections among various algorithms rooted in SPPM is becoming increasingly challenging for the community, both in theory and practical applications. While new variants are yet to be discovered, determining concrete principles beyond intuition to guide their discovery remains challenging. Complicating matters further is the use of various assumptions regarding the correction vectors across different literature, each with differing levels of rigor. Key contributions of this study comprise:

$\bullet$ A universal algorithm. We design a universal SPPM-LC algorithm (stochastic proximal point method with learned correction, Algorithm 1) that encompasses $7$ variants of stochastic proximal point method (SPPM, SPPM-NS, SPPM-AS, SPPM*, SPPM-GC, L-SVRP and Point-SAGA; see Table 1) through (random) correction vectors $h_{k}.$ A specific choice of correction vectors allows to retrieve any particular of the mentioned algorithms. Motivated by the work of Gorbunov et al. (), a similar universal algorithm was proposed by Traoré et al. (), but under different assumptions.

$\bullet$ Comprehensive analysis. We introduce a cohesive theoretical framework of assumptions that includes three restrictions imposed on correction vectors $h_{k}$ of the SPPM-LC algorithm (7). These restrictions connect the correction vectors $h_{k}$ to the iterates of the algorithm and have the forms of parametric recursions. The development of this framework constitutes one of our major contributions. Under 1, 2 on the functions and 7 on the correction vectors, we analyze the convergence of SPPM-LC in the case when it is applied to find an $\varepsilon$-accurate solution of the general stochastic optimization problem (1). This is the first comprehensive analysis of the original, sampling-based, and variance-reduced versions of SPPM. It implies the convergence results for all $7$ variations of stochastic proximal point method (SPPM, SPPM-NS, SPPM-AS, SPPM*, SPPM-GC, L-SVRP and Point-SAGA). The sampling-based methods SPPM-NS and SPPM-AS are employed for a less general setting (3), as well as the variance-reduced method Point SAGA, they are analyzed under 3 and 4 on the functions. In addition, the analysis of L-SVRP is based on the similarity 5, and the analysis of Point-SAGA is based on 6. In order to obtain the convergence result for any particular of the mentioned algorithms, one needs to check that (the recursive parametric) 7 holds for it. Traoré et al. () introduce a slightly different theoretical framework of assumptions and analyze SPPM and its variance-reduced versions in a finite-sum setting (3) in another setup when functions are individually convex and $L$-smooth.

$\bullet$ Best known rates retrieved. The rates derived from our comprehensive Theorem 1, under 1, 2 on the functions for SPPM, SPPM*, SPPM-GC, L-SVRP (under additional 5) and Point-SAGA (under additional 6); under 3 and 4 on the functions for SPPM-NS, SPPM-AS, represent the sharpest rates for these methods. Notably, the neighborhood that we guarantee for SPPM in our analysis does not blow up for the stepsize $\gamma\to\infty$ in comparison to the result of Khaled and Jin () (see (21) in Theorem 2 and Commentary 2 after it). Also, we present the analysis for L-SVRP in a simplified way and obtain slightly better bounds on the iteration complexity up to a constant factor than Khaled and Jin ().

$\bullet$ Analysis under general similarity assumptions. The analysis of L-SVRP is based on the similarity 5, and the analysis of Point-SAGA is based on 6. Both of these assumptions are very general and distinguish our approach from the previous ones. Defazio () analyzed Point SAGA under the individual smoothness assumption, which implies 6. Khaled and Jin () analyzed L-SVRP under a more restrictive similarity assumption than 5. Traoré et al. () analyzed L-SVRP and Point SAGA under the individual smoothness assumption, which is much more restrictive than our 5 and 6. Their results were obtained prior to our work, but we conducted our analysis independently: we already had our result when we discovered their paper.

$\bullet$ Analysis in a general setting. We analyze L-SVRP in the stochastic optimization setting (1) while previous works of Khaled and Jin (); Traoré et al. () do this in the less general setting (3).

$\bullet$ New methods. Our comprehensive theory offers complexity bounds for a range of novel (SPPM-LC, SPPM-NS, SPPM-AS, SPPM-GC, SPPM*) and upcoming variations of SPPM. It suffices to confirm that 7 holds, and a complexity estimate is readily provided by Theorem 1. Selected existing and new methods that align with our framework are outlined in Table 1.

$\bullet$ Experiments. Through extensive experimentation, we demonstrate that several of the newly introduced methods, analyzed within our framework, exhibit compelling empirical properties when compared to natural baselines.

## 4 Main result

We are now ready to present our general Algorithm 1, which we call Stochastic Proximal Point Method with Learned Correction (SPPM-LC). Subsequently, we introduce the core assumption on the correction vectors, iterates and control vectors of Algorithm 1 enabling our general analysis, and further state and comment on our unified convergence result (Theorem 1).

Figure: Algorithm 1 Stochastic Proximal Point Method with Learned Correction (SPPM-LC)

### 4.1 Key assumption

We assume that the (random) correction vector $\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}h_{k}$ has zero mean (conditioned on $x_{k}$ and $\phi_{k},$ the $k$-th iterate and control vector, respectively), and that it is connected with the iterates of SPPM-LC via two parametric recursions/inequalities, described next. We introduce versatility by expressing these inequalities parametrically.

###### Assumption 7 (Parametric recursions) .

Let $\{x_{k},\phi_{k}\}_{k\geq 0}$ be the random iterates produced by SPPM-LC. Assume that

$$ ${\rm E}\left[\left.h_{k}\right|x_{k},\phi_{k}\right]=0.$ (7) $$

Further, assume that there exists a function $\sigma^{2}:\mathbb{R}^{m}\to\mathbb{R}_{+}$ and nonnegative constants ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}},{\color[rgb]{
0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.09765625,0.75390625}B_{1}},{\color[rgb]{
0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.75390625,0.09765625}C_{1}},{\color[rgb]{
0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.75390625,0.09765625,0.09765625}A_{2}},{\color[rgb]{
0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.09765625,0.75390625}B_{2}},{\color[rgb]{
0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.75390625,0.09765625}C_{2}}$, with ${\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{2}}<1$, such that

$$ $\displaystyle{\rm E}\left[\left.{\left\|h_{k}-\nabla f_{\xi_{k}}(x_{\star}) \right\|^{2}}\right|x_{k},\phi_{k}\right]$ $\displaystyle\leq$ $\displaystyle{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named] {pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}{\left\|x_{k}-x_{ \star}\right\|^{2}}+{\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor [named]{pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1}}\sigma^{2} _{k}+{\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{1}},$ (8) $\displaystyle{\rm E}\left[\left.\sigma^{2}_{k+1}\right|x_{k+1},\phi_{k}\right]$ $\displaystyle\leq$ $\displaystyle{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named] {pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}}{\left\|x_{k+1}-x _{\star}\right\|^{2}}+{\color[rgb]{0.09765625,0.09765625,0.75390625} \definecolor[named]{pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{2 }}\sigma^{2}_{k}+{\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[ named]{pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{2}},$ (9) $$

where $\sigma_{k}^{2}:=\sigma^{2}(\phi_{k})$.

For brevity, we refer to this assumption as the “$\sigma_{k}^{2}$ assumption”. If 7 holds, then by taking expectation on both sides of (8) and (9) and applying the tower property in each case,
we get

$$ $\displaystyle{\rm E}\left[{\left\|h_{k}-\nabla f_{\xi_{k}}(x_{\star})\right\|^ {2}}\right]$ $\displaystyle\leq$ $\displaystyle{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named] {pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}{\rm E}\left[{ \left\|x_{k}-x_{\star}\right\|^{2}}\right]+{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{1}}{\rm E}\left[\sigma^{2}_{k}\right]+{ \color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{1}},$ (10) $\displaystyle{\rm E}\left[\sigma^{2}_{k+1}\right]$ $\displaystyle\leq$ $\displaystyle{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named] {pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}}{\rm E}\left[{ \left\|x_{k+1}-x_{\star}\right\|^{2}}\right]+{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}{\rm E}\left[\sigma^{2}_{k}\right]+{ \color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{2}},$ (11) $$

The novelty of our approach lies in the introduction of inequalities (8) and (9). We support and validate this assertion by providing numerous examples (in Section 5), demonstrating that these inequalities encapsulate the nature of a broad range of existing SPPM methods as well as some new ones, including standard SPPM alongside its arbitrary sampling and variance-reduced variants. In its essence, we generalize, parameterize and establish as an independent assumption the conditions on correction vectors for SPPM-type methods present in the literature, regardless of the specifics defining the base method from which they stem. Traoré et al. () propose different parameterized recursive inequalities and analyze SPPM and its variance-reduced versions in a finite-sum setting (3) under the condition where functions are individually convex and $L$-smooth. Similar inequalities can be found in the analysis of SGD-type methods (a unified theory developed by Gorbunov et al. ()).

###### Assumption 7 (Parametric recursions) .

### 4.2 Main theorem

We are now prepared to introduce our main convergence result.

###### Theorem 1 .

Let 1 (differentiability) and 2 ($\mu$-strong convexity) hold.
Let $\{x_{k},h_{k}\}$ be the iterates produced by SPPM-LC (Algorithm 1), and assume that they satisfy
7 ($\sigma_{k}^{2}$-assumption).
Choose any $\gamma>0$ and $\alpha>0$ satisfying the inequalities

$$ $\frac{(1+\gamma^{2}{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[ named]{pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}})(1+\alpha{ \color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}})}{(1+\gamma\mu)^{ 2}}<1,\qquad\frac{\gamma^{2}{\color[rgb]{0.09765625,0.09765625,0.75390625} \definecolor[named]{pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1 }}(1+\alpha{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}})}{\alpha(1+\gamma \mu)^{2}}+{\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{2}}<1,$ (12) $$

and define the Lyapunov function

$$ $\Psi_{k}:={\left\|x_{k}-x_{\star}\right\|^{2}}+\alpha\sigma^{2}_{k}.$ (13) $$

Then for all iterates $k\geq 0$ of SPPM-LC we have

$$ ${\rm E}\left[\Psi_{k}\right]\leq\theta^{k}\Psi_{0}+\frac{\zeta}{1-\theta},$ (14) $$

where the parameters $0\leq\theta<1$ and $\zeta\geq 0$ are given by

$$ $\theta=\max\left\{\frac{(1+\gamma^{2}{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{1}})(1+\alpha{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{2}})}{(1+\gamma\mu)^{2}},\frac{\gamma^{2}{ \color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1}}(1+\alpha{\color[ rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb} {0.75390625,0.09765625,0.09765625}A_{2}})}{\alpha(1+\gamma\mu)^{2}}+{\color[ rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb} {0.09765625,0.09765625,0.75390625}B_{2}}\right\},$ (15) $$

$$ $\zeta=\frac{\gamma^{2}{\color[rgb]{0.09765625,0.75390625,0.09765625} \definecolor[named]{pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{1 }}(1+\alpha{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}})}{(1+\gamma\mu)^{ 2}}+\alpha{\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{2}}.$ (16) $$

Theorem 1 proves a linear rate of convergence for a number of stochastic proximal point methods towards a fluctuation neighborhood around the solution, regulated by the additive term in (14). It depends on parameters ${\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{1}}$ and ${\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{2}}.$ The neighborhood vanishes (i.e., $\zeta=0$) iff ${\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{1}}={\color[rgb]{
0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.75390625,0.09765625}C_{2}}=0$. If this happens, then SPPM-LC converges linearly to the solution as $k\to\infty$ for any fixed $\gamma>0,$ satisfying the conditions of Theorem 1. Notice that $\theta\geq{\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{2}}$, and hence the linear rate can not be faster than $({\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{2}})^{k}.$ That is, as shown in Appendix B (also, see Table 2), the main difference between variance-reduced versions of SPPM and its other variants is that the former methods satisfy $\sigma^{2}_{k}$-assumption with ${\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{1}}={\color[rgb]{
0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.75390625,0.09765625}C_{2}}=0$ (and reach the optimum $x_{\star}$), whilst the latter have either ${\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{1}}>0$ or ${\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{2}}>0.$

**Table 1: Compilation of both existing and novel methods that align with our comprehensive analytical framework. Problem (1) encompasses a broader scope compared to problem (3). VR = variance reduced method, AS = arbitrary sampling. Thm $x$ = Thm $1$ + Lemma $x$, Lemma $x$ is in Section A.$x,$ $x\in[8].$ The last column indicates whether the analysis is new or not and whether we recover the previously established rate.**
| Problem | Method | Alg # | VR? | AS? | Section | Result  / Rate |
| --- | --- | --- | --- | --- | --- | --- |
| (1) | SPPM-LC [NEW] | Alg 1 | ✓✗ | ✗ | B.1 | Thm 1 [NEW] |
| (1) | SPPM ^(1)  (Bertsekas, ) | Alg 2 | ✗ | ✗ | B.2 | Thm 2 [NEW] |
| (3) | SPPM-NS [NEW] | Alg 3 | ✗ | ✗ | B.3 | Thm 3 [NEW] |
| (3) | SPPM-AS [NEW] | Alg 4 | ✗ | ✓ | B.4 | Thm 4 [NEW] |
| (1) | SPPM* [NEW] | Alg 5 | ✓ | ✓ | B.5 | Thm 5 [NEW] |
| (1) | SPPM-GC [NEW] | Alg 6 | ✓ | ✗ | B.6 | Thm 6 [NEW] |
| (1) | L-SVRP ^(2) (Khaled and Jin, ) | Alg 7 | ✓ | ✗ | B.7 | Thm 7 [NEW] |
| (3) | Point SAGA ^(3) (Defazio, ) | Alg 8 | ✓ | ✗ | B.8 | Thm 9 [NEW] |

###### Theorem 1 .

## 5 Overview of specific methods and of the framework

In this section we demonstrate the descriptive power of our new framework. We show that most popular existing methods can be expressed in terms of our framework and, in addition to that, we describe several new methods in more detail (see Table 1).

### 5.1 A brief overview

As asserted, the proposed framework is powerful enough to include methods without variance reduction (✗  in the “VR” column) alongside variance-reduced methods (✓  in the “VR” column), methods that fall under the arbitrary sampling paradigm (✓  in the “AS” column). All novel methods introduced in this paper are clearly designated with the label [NEW]. Additionally, to facilitate a thorough understanding of all algorithms discussed, detailed explanations are included in the Appendix. A link is provided for convenient navigation to the supplementary section. The generality of our framework is reflected in Table 1. The “Result / Rate” column of Table 1 refers to a Theorem $x$ which follows from Theorem 1 and Lemma $x,$ $x\in[8].$ The convergence results of the algorithms considered in the paper are outlined in these theorems, offering insights into their performance characteristics. Importantly, in instances where established methods are recovered, the best-known convergence rates are reaffirmed or better results are obtained. This underscores the robustness and reliability of our analytical framework in accurately capturing the behavior of established algorithms.

### 5.2 Parameters of the framework

The algorithms in Table 1 demonstrate specific patterns in relation to the parameters in 7. To elucidate this observation, we provide a summary of these parameters in Table 2. As predicted by Theorem 1, when ${\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{1}}={\color[rgb]{
0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.75390625,0.09765625}C_{2}}=0,$ the corresponding method does not oscillate and converges to the optimum $x_{\star},$ indicating the variance-reduced nature of the algorithm. All parameters referenced in Table 2 are defined in the Appendix, alongside descriptions and analyses of the respective methods.

### 5.3 Five novel algorithms

To showcase the efficiency of our comprehensive framework, we develop three new variants of SPPM which have not previously been addressed in the literature (see Table 1). In this section, we briefly outline the reasoning behind their implementation. Further specifics are available in the Appendix.

SPPM-NS (Algorithm 3). The method is designed for solving the problem (3). Let positive numbers $p_{1},\ldots,p_{n}$ sum up to $1,$ set $i_{k}=i$ with probability $p_{i},$ $i\in[n].$ The step of the method has the form $x_{k+1}={\rm prox}_{\frac{\gamma}{np_{i}}f_{i_{k}}}\left(x_{k}\right).$ It unifies several powerful sampling strategies (as, e.g., importance sampling). Sampling allows to improve the convergence rate and modify the neighborhood.

SPPM-AS (Algorithm 4). The method is also designed for the problem (3). The arbitrary sampling framework was developed for SGD by Gower et al. (). It allows to conduct a sharp unified convergence analysis for various effective sampling and mini-batch strategies. For strongly convex functions, the method with constant stepsize converges linearly to the neighborhood of the solution.

SPPM* (Algorithm 5). This novel algorithm links conventional and variance-reduced SPPM methods. Although not immediately practical, it offers valuable insights into the inner workings of variance reduction. This method addresses the fundamental question: assuming that the gradients $\nabla f_{\xi}(x_{\star})$ are available, can they be leveraged to devise a more potent variant of SPPM? The affirmative answer culminates in the development of SPPM*. The construction of updates in SPPM* involves correction vectors of the form $h_{k}=\nabla f_{\xi_{k}}(x_{\star}).$ In essence, this implies augmenting $x_{k}$ with gradients of the same functions at the optimal point $x_{\star}$, with respect to which the proximal operator is calculated. As evidenced in Table 2, where ${\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{1}}={\color[rgb]{
0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.75390625,0.09765625}C_{2}}=0,$ this method converges directly to $x_{\star}$ without oscillation, rather than converging to a neighborhood of the solution. Practical variance-reduced methods operate by iteratively refining estimates of $\nabla f_{\xi_{k}}(x_{\star}).$ Notably, the term $\sigma_{k}^{2}$ in the Lyapunov function of variance-reduced methods converges to zero.

SPPM-GC (Algorithm 6). This method can be viewed as a practical variant of SPPM* featuring a computable version of the correction vector, $h_{k}=\nabla f_{\xi_{k}}(x_{k})-\nabla f(x_{k}),$ instead of the incomputable correction $\nabla f_{\xi_{k}}(x_{\star}).$ This method follows the same paradigm of iteratively constructing increasingly more refined estimates of $\nabla f_{\xi_{k}}(x_{\star})$ and is also variance reduced. As indicated in Table 2, where ${\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{1}}={\color[rgb]{
0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.75390625,0.09765625}C_{2}}=0,$ this method also converges to $x_{\star},$ and not to some neighborhood of the solution only.

SPPM-LC (Algorithm 1). This new generic algorithm constructs the updates with correction vectors $h_{k}$ of a universal form. We analyze its behavior under a new parametric 7. The convergence result then follows. This is a unified analysis of all stochastic proximal point methods we have encountered so far: SPPM, SPPM-NS, SPPM-AS, SPPM*, SPPM-GC, L-SVRP and Point-SAGA. One only needs to check that this parametric assumption holds in each particular case (see Table 2).

**Table 2: The parameters determining the compliance of the methods listed in Table 1 to Assumption 7. Detailed explanations of the expressions featured in the table are provided in the Appendix.**
| Method | ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{<br>pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}$ | ${\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{<br>pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1}}$ | ${\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{<br>pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{1}}$ | ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{<br>pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}}$ | ${\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{<br>pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{2}}$ | ${\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{<br>pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{2}}$ | Lemma |
| --- | --- | --- | --- | --- | --- | --- | --- |
| SPPM-LC | ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{<br>pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}$ | ${\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{<br>pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1}}$ | ${\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{<br>pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{1}}$ | ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{<br>pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}}$ | ${\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{<br>pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{2}}$ | ${\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{<br>pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{2}}$ | Lemma 1 |
| SPPM | 0 0 | 0 0 | $\sigma_{\star}^{2}$ | 0 0 | 0 0 | 0 0 | Lemma 2 |
| SPPM-NS | 0 0 | 0 0 | $\sigma_{\star,{\rm NS}}^{2}$ | 0 0 | 0 0 | 0 0 | Lemma 3 |
| SPPM-AS | 0 0 | 0 0 | $\sigma_{\star,{\rm AS}}^{2}$ | 0 0 | 0 0 | 0 0 | Lemma 4 |
| SPPM* | 0 0 | 0 0 | 0 0 | 0 0 | 0 0 | 0 0 | Lemma 5 |
| SPPM-GC | $\delta^{2}$ | 0 0 | 0 0 | 0 0 | 0 0 | 0 0 | Lemma 6 |
| L-SVRP | 0 0 | $\delta^{2}$ | 0 0 | $p$ | $1-p$ | 0 0 | Lemma 7 |
| Point SAGA | 0 0 | $\nu^{2}$ | 0 0 | $\frac{1}{n}$ | $\frac{n-1}{n}$ | 0 0 | Lemma 8 |

## 6 Experiments

In this section we describe numerical experiments conducted for the linear regression problem

$$ $\min_{x\in\mathbb{R}^{d}}\left\{\frac{1}{2n}\sum^{n}_{i=1}(a^{\top}_{i}x-b_{i} )^{2}+\lambda_{i}{\left\|x\right\|^{2}}\right\},$ (17) $$

where $a_{i}\in\mathbb{R}^{d},~{}b_{i}\in\mathbb{R}^{d}$ is the $i$-th data pair, each $\lambda_{i}$ is a $\ell_{2}$-regularization parameter. We provide $4$ sets of experiments. The first one is devoted to the comparison of different sampling strategies for SPPM-NS : uniform sampling US, importance sampling IS, variance sampling VS (see Section  B.3), with three selected stepsizes. The second set of experiments demonstrates the change of the radius of neighborhood for SPPM-AS with $\tau$-nice sampling. In the next two sets of experiments we illustrate the main difference between SPPM and SPPM with variance reduction. Also in practice we show the relationship between SPPM-GC, L-SVRP, Poin-SAGA as we did in theory. For the first bunch of experiments we set $n=10$, $d=3$ and the regularization parameters $\lambda_{i}=\nicefrac{{1}}{{2^{i}}}$, where $i\in[d]$. Looking at Figure 1, we can see the different behavior of the methods. In the first two plots with the smallest stepsizes, SPPM-IS has a faster start, but a larger neighborhood than the other considered methods. Unfortunately, we cannot say that SPPM-VS has a much smaller neighborhood radius than SPPM-US or SPPM-IS, but theoretically, the variance is smaller. In the second set of numerical experiments (see Figure 2), we observe a clear correlation between the neighborhood radius and the cardinality of the sampled subset $\tau$. More precisely, the larger $\tau$ is, the smaller the neighborhood radius is. In the third set of numerical experiments (see Figure 3), we set $n=1000$, $d=10$ and each $\lambda_{i}=1$ for the problem (17). On all three plots we observe the superiority of SPPM-star in terms of the stepsize choice. For example, in the third plot with $\gamma=10^{2}$, SPPM-GC diverges. In the fourth set of numerical experiments (see Figure 4) with parameters $n=1000$, $d=10$ and each $\lambda_{i}=1$ for the problem (17), we observe that the performances of SPPM-GC and L-SVRP with $p=1$ are identical, which is supported by our theoretical findings. With decreasing $p$ from $1$ to $\nicefrac{{1}}{{n}}$ we see how the behavior of L-SVRP worsens and matches with the performance of Point-SAGA when $p=\nicefrac{{1}}{{n}}$.

Figure: Figure 1: Comparison of the performance of SPPM-US, SPPM-IS, SPPM-VS and SPPM-AS with $\tau=9$-nice sampling for different selections of stepsize $\gamma\in\{10^{-4},10^{-2},1,10^{2}\}$.
Refer to caption: x1.png

Figure: Figure 2: Comparison of the performance of SPPM-AS with $\tau$-nice sampling with different selections of cardinality $\tau\in\{1,2,5,9,n=10\}$ and stepsize $\gamma\in\{10^{-2},10^{-1},1\}$.
Refer to caption: x5.png

Figure: Figure 3: Comparison of the performance of SPPM-US , SPPM-GC and SPPM-star with different selections of stepsize $\gamma\in\{10^{-2},1,10^{2}\}.$
Refer to caption: x8.png

Figure: Figure 4: Comparison of the performance of SPPM-GC , Point-SAGA , L-SVRP with different selection of probabilities $\gamma\in\{\nicefrac{{1}}{{n}}=10^{-3},5\cdot 10^{-3},10^{2},5\cdot 10^{2},10^ {-1},1\}$ . The stepsizes are taken according to the theory.
Refer to caption: x11.png

## 7 Further discussion

Although our approach is general, we still see some limitations, open problems and several possible directions for future extensions. Generating a similar result in the nonconvex case continues to be an unsolved challenge. Expanding 7 to incorporate iteration-dependent parameters ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}},{\color[rgb]{
0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.09765625,0.75390625}B_{1}},{\color[rgb]{
0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.75390625,0.09765625}C_{1}},{\color[rgb]{
0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.75390625,0.09765625,0.09765625}A_{2}},{\color[rgb]{
0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.09765625,0.75390625}B_{2}},{\color[rgb]{
0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.75390625,0.09765625}C_{2}}$ could facilitate the development of various novel methods as SPPM with decreasing stepsizes. Khaled and Jin () explore federated learning, investigating client sampling to enhance communication efficiency. Another way to do that is to incorporate the utilization of compressed vectors exchanged between a server and the clients. This motivates the problem of providing the analysis for such SPPM-type methods and incorporating them into our framework. It would be interesting to build theory for algorithms with correction vectors $h_{k}$ with a non-zero expected value and unify with our theory. Another potential avenue for future research involves offering a comprehensive analysis of SPPM-type methods incorporating acceleration and momentum.

## Appendix A Extended literature overview

Stochastic gradient descent (SGD) Robbins and Monro (); Nemirovski et al. (); Bottou () is a contemporary common algorithm for solving optimization problems (1) and (3) with updates of the form $x^{k+1}=x^{k}-\gamma g^{k},$ where $g^{k}$ is an unbiased stochastic gradient estimator: ${\rm E}\left[\left.g^{k}\right|x^{k}\right]=\nabla f(x^{k}).$ Theoretical properties of SGD are nowadays studied well in works by Vaswani et al. (); Gower et al. (); Gorbunov et al. (); Khaled and Richtárik (); Demidovich et al. (). The versatility in the design of the estimator $g^{k}$ has led to a development of highly effective variants of SGD based on importance sampling (see works by Needell et al. (); Zhao and Zhang ()), mini-batching (see a work of Konečný et al. ()). These strategies are unified in the arbitrary sampling paradigm proposed by Gower et al. (). The iterates of standard SGD as well as of its sampling variants do not converge to the optimum due to the presence of the stochastic gradient noise. The existence of this issue led to the construction of variance reduction methods in works by Roux et al. (); Defazio et al. (); Johnson and Zhang (); Kovalev et al. (); Nguyen et al. (). Such methods are able to sequentially learn the stochastic gradients at the optimum which allows them to reach a linear convergence to $x_{\star}$ when equipped with constant stepsize in case of strongly convex $f.$ The work by Gorbunov et al. () provided a unified theory for SGD and its variance-reduced, sampling, quantization and coordinate sub-sampling modifications. Problem setting (3) can also be used to describe a Federated Learning setup where $n$ is the number of workers, and $f_{i}$ stands for the loss dependent on the data of worker $i\in[n].$ Workers compute individual stochastic gradients and aggregate them on a master node (see a paper by Konečný and Richtárik ()). A primary bottleneck in such settings is communication. To address it, many techniques are used as quantization (see works by Gupta et al. (); Seide et al. ()), sparsification (see a paper by Alistarh et al. ()), dithering (see a work by Alistarh et al. ()). Many distributed optimization works employ variants of 5 to analyze methods that solve the problem (3) in the strongly convex case (see works by Shamir et al. (); Zhang and Xiao (); Yuan and Li (); Karimireddy et al. ()). Szlendak et al. (); Beznosikov and Gasnikov (); Panferov et al. () achieve better communication complexity guarantees in a Federated Learning setting under the similarity assumption employing specific quantizers and sparsifiers. Selecting the right stepsize is a critical aspect in the implementation of SGD, see the work of Bach and Moulines (). In Section 2 we show that the iteration complexity of SGD depends on condition number $\kappa=\frac{L}{\mu}.$ If we are able to evaluate stochastic proximal operator, another algorithm to use instead of SGD is stochastic proximal point method (SPPM). Its iteration complexity does not depend on the smoothness constant $L.$ Bertsekas () calls it the incremental proximal point method and considers it for solving the problem (3), proving a convergence to the neighborhood when each $f_{i}$ is Lipschitz. In the study by Ryu and Boyd (), convergence rates are presented for SPPM, emphasizing its robustness against learning rate misspecification, a characteristic not exhibited by SGD. Pătraşcu and Necoara () consider SPPM to solve a general stochastic optimization problem (1) in the convex case and provide nonasymptotic convergence guarantees. Asi and Duchi () delve into the exploration of a broader method termed AProx, which includes SPPM as a particular case. Their research entails analysis of convergence rates for convex functions. Variance-reduced versions of SPPM should also have advantages over their SGD counterparts. It motivated Khaled and Jin () to analyze L-SVRP in the finite-sum and federated settings under a stronger version of 5 and strong convexity. Defazio () analyzes Point SAGA under the individual smoothness and strong convexity assumptions. Traoré et al. () consider variance-reduced methods L-SVRP and Point SAGA when each $f_{i}(x)$ is $L$-smooth, convex, and either $f(x)$ is convex or satisfies PŁ-condition. Kim et al. () consider SPPM with momentum.

## Appendix B Special cases

### B.1 Stochastic proximal point method with learned correction ( SPPM-LC )

This section is devoted to a unified analysis of all stochastic proximal point methods we have encountered so far: SPPM, SPPM-NS, SPPM-AS, SPPM*, SPPM-GC, L-SVRP and Point-SAGA. The analysis is based on formulating a new parametric assumption (7) on the behavior of the method, and the convergence result then follows. All that has to be done is to check that this parametric assumption holds in each particular case of interest.

###### Lemma 1 .

7 holds with constants ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}},{\color[rgb]{
0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.09765625,0.75390625}B_{1}},{\color[rgb]{
0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.75390625,0.09765625}C_{1}}\text{ and }{\color[rgb]{
0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.75390625,0.09765625,0.09765625}A_{2}},{\color[rgb]{
0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.09765625,0.75390625}B_{2}},{\color[rgb]{
0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.75390625,0.09765625}C_{2}}.$

###### Proof of Lemma 1 .

Since we assume the iterates produced by SPPM-LC (Algorithm 1) satisfy 7, the statement of Lemma 1 holds automatically.
∎

We can now present the main result of this section.

Theorem 1.
Let 1 (differentiability) and 2 ($\mu$-strong convexity) hold.
Let $\{x_{k},h_{k}\}$ be the iterates produced by SPPM-LC (Algorithm 1), and assume that they satisfy
7 ($\sigma_{k}^{2}$).
Choose any $\gamma>0$ and $\alpha>0$ satisfying the inequalities

$$ $\frac{(1+\gamma^{2}{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[ named]{pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}})(1+\alpha{ \color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}})}{(1+\gamma\mu)^{ 2}}<1,\qquad\frac{\gamma^{2}{\color[rgb]{0.09765625,0.09765625,0.75390625} \definecolor[named]{pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1 }}(1+\alpha{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}})}{\alpha(1+\gamma \mu)^{2}}+{\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{2}}<1,$ (18) $$

and define the Lyapunov function

$$ $\Psi_{k}:={\left\|x_{k}-x_{\star}\right\|^{2}}+\alpha\sigma^{2}_{k}.$ (19) $$

Then for all iterates $k\geq 0$ of SPPM-LC we have

$$ ${\rm E}\left[\Psi_{k}\right]\leq\theta^{k}\Psi_{0}+\frac{\zeta}{1-\theta},$ $$

where the parameters $0\leq\theta<1$ and $\zeta\geq 0$ are given by

$$ $\theta=\max\left\{\frac{(1+\gamma^{2}{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{1}})(1+\alpha{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{2}})}{(1+\gamma\mu)^{2}},\frac{\gamma^{2}{ \color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1}}(1+\alpha{\color[ rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb} {0.75390625,0.09765625,0.09765625}A_{2}})}{\alpha(1+\gamma\mu)^{2}}+{\color[ rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb} {0.09765625,0.09765625,0.75390625}B_{2}}\right\},\qquad\zeta=\frac{\gamma^{2}{ \color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{1}}(1+\alpha{\color[ rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb} {0.75390625,0.09765625,0.09765625}A_{2}})}{(1+\gamma\mu)^{2}}+\alpha{\color[ rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb} {0.09765625,0.75390625,0.09765625}C_{2}}.$ $$

###### Proof of Theorem 1 ..

We use Lemma 9 (see Appendix C) and $\sigma^{2}_{k}$-assumption. In particular, by combining inequality 7 (see (9)) and Lemma 9 (see (70)), for all $\gamma>0$ we get

$$ $\displaystyle{\rm E}\left[\Psi_{k+1}\right]$ $\displaystyle\overset{\eqref{eq:Lyapunov_func_sigma_k}}{=}$ $\displaystyle{\rm E}\left[{\left\|x_{k+1}-x_{\star}\right\|^{2}}+\alpha\sigma^ {2}_{k+1}\right]$ $\displaystyle=$ $\displaystyle{\rm E}\left[{\left\|x_{k+1}-x_{\star}\right\|^{2}}\right]+\alpha {\rm E}\left[\sigma^{2}_{k+1}\right]$ $\displaystyle\overset{\eqref{eq:sigma_k_1}}{\leq}$ $\displaystyle\frac{(1+\gamma^{2}{\color[rgb]{0.75390625,0.09765625,0.09765625} \definecolor[named]{pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1 }})}{(1+\gamma\mu)^{2}}{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right ]+\frac{\gamma^{2}{\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[ named]{pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1}}}{(1+\gamma \mu)^{2}}{\rm E}\left[\sigma^{2}_{k}\right]+\frac{\gamma^{2}{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}}{(1+\gamma\mu)^{2}}+\alpha{\rm E}\left [\sigma^{2}_{k+1}\right]$ $\displaystyle\overset{\eqref{eq:sigma_k_assumption_1}}{\leq}$ $\displaystyle\frac{(1+\gamma^{2}{\color[rgb]{0.75390625,0.09765625,0.09765625} \definecolor[named]{pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1 }})}{(1+\gamma\mu)^{2}}{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right ]+\frac{\gamma^{2}{\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[ named]{pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1}}}{(1+\gamma \mu)^{2}}{\rm E}\left[\sigma^{2}_{k}\right]+\frac{\gamma^{2}{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}}{(1+\gamma\mu)^{2}}$ $\displaystyle\qquad+\alpha{\color[rgb]{0.75390625,0.09765625,0.09765625} \definecolor[named]{pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2 }}{\rm E}\left[{\left\|x_{k+1}-x_{\star}\right\|^{2}}\right]+\alpha{\color[rgb ]{0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}{\rm E}\left[\sigma^{2}_{k}\right]+ \alpha{\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{2}}$ $\displaystyle\overset{\eqref{eq:sigma_k_1}}{\leq}$ $\displaystyle\frac{(1+\gamma^{2}{\color[rgb]{0.75390625,0.09765625,0.09765625} \definecolor[named]{pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1 }})}{(1+\gamma\mu)^{2}}{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right ]+\frac{\gamma^{2}{\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[ named]{pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1}}}{(1+\gamma \mu)^{2}}{\rm E}\left[\sigma^{2}_{k}\right]+\frac{\gamma^{2}{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}}{(1+\gamma\mu)^{2}}$ $\displaystyle\qquad+\frac{(1+\gamma^{2}{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{1}})\alpha{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{2}}}{(1+\gamma\mu)^{2}}{\rm E}\left[{\left \|x_{k}-x_{\star}\right\|^{2}}\right]+\frac{\gamma^{2}{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{1}}\alpha{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{2}}}{(1+\gamma\mu)^{2}}{\rm E}\left[\sigma ^{2}_{k}\right]+\frac{\gamma^{2}{\color[rgb]{0.09765625,0.75390625,0.09765625} \definecolor[named]{pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{1 }}\alpha{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}}}{(1+\gamma\mu)^{2}}$ $\displaystyle\qquad+\alpha{\color[rgb]{0.09765625,0.09765625,0.75390625} \definecolor[named]{pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{2 }}{\rm E}\left[\sigma^{2}_{k}\right]+\alpha{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}$ $\displaystyle=$ $\displaystyle\frac{(1+\gamma^{2}{\color[rgb]{0.75390625,0.09765625,0.09765625} \definecolor[named]{pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1 }})(1+\alpha{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}})}{(1+\gamma\mu)^{ 2}}{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]+\left(\frac{\gamma ^{2}{\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1}}(1+\alpha{\color[ rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb} {0.75390625,0.09765625,0.09765625}A_{2}})}{\alpha(1+\gamma\mu)^{2}}+{\color[ rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb} {0.09765625,0.09765625,0.75390625}B_{2}}\right)\alpha{\rm E}\left[\sigma^{2}_{ k}\right]$ $\displaystyle\qquad+\frac{\gamma^{2}{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}(1+\alpha{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{2}})}{(1+\gamma\mu)^{2}}+\alpha{\color[rgb ]{0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}$ $\displaystyle\leq$ $\displaystyle\underbrace{\max\left\{\frac{(1+\gamma^{2}{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{1}})(1+\alpha{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{2}})}{(1+\gamma\mu)^{2}},\frac{\gamma^{2}{ \color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1}}(1+\alpha{\color[ rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb} {0.75390625,0.09765625,0.09765625}A_{2}})}{\alpha(1+\gamma\mu)^{2}}+{\color[ rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb} {0.09765625,0.09765625,0.75390625}B_{2}}\right\}}_{:=\theta}{\rm E}\left[V_{k}\right]$ $\displaystyle+$ $\displaystyle\underbrace{\frac{\gamma^{2}{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}(1+\alpha{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{2}})}{(1+\gamma\mu)^{2}}+\alpha{\color[rgb ]{0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}}_{:=\zeta}$ $\displaystyle=$ $\displaystyle\theta{\rm E}\left[\Psi_{k}\right]+\zeta.$ $$

By unrolling the recurrence (we can also just apply 3), we get

$$ $\displaystyle{\rm E}\left[\Psi_{k}\right]\leq\theta^{k}\Psi_{0}+\sum^{k-1}_{t= 0}\theta^{k-1-t}\zeta=\theta^{k}\Psi_{0}+\sum^{k-1}_{t=0}\theta^{t}\zeta\leq \theta^{k}\Psi_{0}+\sum^{\infty}_{t=0}\theta^{t}\zeta=\theta^{k}\Psi_{0}+\frac {1}{1-\theta}\zeta.$ $$

∎

###### Lemma 1 .

###### Proof of Lemma 1 .

###### Proof of Theorem 1 ..

### B.2 Stochastic proximal point method ( SPPM )

We aim to solve problem (1) via the stochastic proximal point method (SPPM), formalized as Algorithm 2. The main step has the form

$$ ${\rm prox}_{\gamma f_{\xi_{k}}}\left(x_{k}\right):=\arg\min_{x\in\mathbb{R}^{d }}\left\{f_{\xi_{k}}(x)+\frac{1}{2\gamma}{\left\|x-x_{k}\right\|^{2}}\right\},$ $$

where $\gamma>0$ is a stepsize.

Figure: Algorithm 2 Stochastic Proximal Point Method (SPPM)

Notice that if $\gamma$ is kept “too large”, then $x_{k+1}\approx\arg\min_{x}f_{\xi_{k}}(x)$, and hence the method just “samples” the minimizers of the stochastic functions, and does not progress towards finding $x_{\star}$ (unless, of course, $x_{\star}$ happens to be shared by all functions $f_{\xi_{k}}$). As we shall see, the situation is different if $\gamma$ is kept sufficiently small.

###### Lemma 2 ( SPPM ) .

7 holds for the iterates of SPPM (Algorithm 2) with

$$ ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{1}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}=\sigma^{2}_{\star},\quad\text{and} \quad{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}=0.$ $$

###### Proof of Lemma 2 ..

Recall that the iterates of SPPM have the form

$$ $x_{k+1}={\rm prox}_{\gamma f_{\xi_{k}}}\left(x_{k}\right).$ $$

Thus, $h_{k}=0$. Let $\sigma^{2}\equiv 0$. Clearly, (7) holds. Furthermore,

$$ ${\rm E}\left[\left.{\left\|h_{k}-\nabla f_{\xi_{k}}(x_{\star})\right\|^{2}} \right|x_{k},\sigma_{k}^{2}\right]={\rm E}\left[{\left\|\nabla f_{\xi_{k}}(x_{ \star})\right\|^{2}}\right]:=\sigma^{2}_{\star}.$ $$

Hence, 7 holds with

$$ ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{1}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}=\sigma^{2}_{\star},\quad\text{and} \quad{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}=0.$ $$

∎

The convergence of SPPM is captured by the following theorem.

###### Theorem 2 .

Let 1 (differentiability) and 2 ($\mu$-strong convexity) hold and define

$$ $\sigma_{\star}^{2}:={\rm E}_{\xi\sim\mathcal{D}}\left[{\left\|\nabla f_{\xi}(x _{\star})\right\|^{2}}\right].$ (20) $$

Let $x_{0}\in\mathbb{R}^{d}$ be an arbitrary starting point. Then for any $k\geq 0$ and any $\gamma>0$, the iterates of SPPM (Algorithm 2) satisfy

$$ $\displaystyle{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\left (\frac{1}{1+\gamma\mu}\right)^{2k}{\left\|x_{0}-x_{\star}\right\|^{2}}+\frac{ \gamma\sigma_{\star}^{2}}{\gamma\mu^{2}+2\mu}.$ (21) $$

Commentary:

- 1.
Interpolation regime. Consider the interpolation regime, characterized by $\sigma_{\star}^{2}=0$. Since we can use arbitrarily large $\gamma>0$, we obtain an arbitrarily fast convergence rate:
${\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\left(\frac{1}{1+
\gamma\mu}\right)^{2k}{\left\|x_{0}-x_{\star}\right\|^{2}}.$
(22)
Indeed, $\left(\frac{1}{1+\gamma\mu}\right)^{2k}$ can be made arbitrarily small for any fixed $k\geq 1$, even $k=1$, by choosing $\gamma$ large enough. However, this is not surprising, since now $f$ and all functions $f_{\xi}$ share a single minimizer, $x_{\star}$, and hence it is possible to find it by sampling a single function $f_{\xi}$ and minimizing it, which is what the prox does, as long as $\gamma$ is large enough.
- 2.
A single step travels far. Observe that for $\gamma=\frac{1}{\mu}$, we have $\frac{\gamma\sigma_{\star}^{2}}{\gamma\mu^{2}+2\mu}=\frac{\sigma_{\star}^{2}}{
3\mu^{2}}$. In fact, the convergence neighborhood $\frac{\gamma\sigma_{\star}^{2}}{\gamma\mu^{2}+2\mu}$ is bounded above by three times this quantity irrespective of the choice of the stepsize. Indeed,
$\frac{\gamma\sigma_{\star}^{2}}{\gamma\mu^{2}+2\mu}\leq\min\left\{\frac{\sigma
_{\star}^{2}}{\mu^{2}},\frac{\gamma\sigma_{\star}^{2}}{\mu}\right\}\leq\frac{
\sigma_{\star}^{2}}{\mu^{2}}.$
This means that no matter how far the starting point $x_{0}$ is from the optimal solution $x_{\star}$, if we choose the stepsize $\gamma$ to be large enough, then we can get a decent-quality solution after a single iteration of SPPM already! Indeed, if we choose $\gamma$ large enough so that
$\left(\frac{1}{1+\gamma\mu}\right)^{2}{\left\|x_{0}-x_{\star}\right\|^{2}}\leq\delta,$
where $\delta>0$ is chosen arbitrarily, then for $k=1$ we get
$\displaystyle{\rm E}\left[{\left\|x_{1}-x_{\star}\right\|^{2}}\right]$
$\displaystyle\leq$
$\displaystyle\delta+\frac{\sigma_{\star}^{2}}{\mu^{2}}.$
- 3.
Iteration complexity. We have seen above that accuracy arbitrarily close to (but not reaching) $\frac{\sigma_{\star}^{2}}{\mu^{2}}$ can be achieved via a single step of the method, provided the stepsize $\gamma$ is large enough. Assume now that we aim for $\varepsilon$ accuracy where $\varepsilon\leq\frac{\sigma_{\star}^{2}}{\mu^{2}}$. Using the inequality $1-t\leq\exp(-t)$ which holds for all $t>0$, we get
$\left(\frac{1}{1+\gamma\mu}\right)^{2k}=\left(1-\frac{\gamma\mu}{1+\gamma\mu}
\right)^{2k}\leq\exp\left(-\frac{2\gamma\mu k}{1+\gamma\mu}\right).$
Therefore, provided that
$k\geq\frac{1+\gamma\mu}{2\gamma\mu}\log\left(\frac{2{\left\|x_{0}-x_{\star}
\right\|^{2}}}{\varepsilon}\right),$
we get
$\left(\frac{1}{1+\gamma\mu}\right)^{2k}{\left\|x_{0}-x_{\star}\right\|^{2}}
\leq\frac{\varepsilon}{2}.$ Furthermore, as long as $\gamma\leq\frac{2\varepsilon\mu}{2\sigma_{\star}^{2}-\varepsilon\mu^{2}}$ (this is true provided that the more restrictive but also more elegant-looking condition $\gamma\leq\varepsilon\frac{\mu}{\sigma_{\star}^{2}}$ holds), we get
$\frac{\gamma\sigma_{\star}^{2}}{\gamma\mu^{2}+2\mu}\leq\frac{\varepsilon}{2}.$
Putting these observations together, we conclude that with the stepsize $\gamma=\varepsilon\frac{\mu}{\sigma_{\star}^{2}}$, we get
$\displaystyle{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\varepsilon$
provided that
$k\geq\frac{1+\gamma\mu}{2\gamma\mu}\log\frac{2{\left\|x_{0}-x_{\star}\right\|^
{2}}}{\varepsilon}=\left(\frac{\sigma_{\star}^{2}}{2\varepsilon\mu^{2}}+\frac{
1}{2}\right)\log\left(\frac{2{\left\|x_{0}-x_{\star}\right\|^{2}}}{\varepsilon
}\right).$
(23)

###### Proof of Theorem 2 .

Recall that in the case of SPPM, we have ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=0,{\color[rgb]{
0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.09765625,0.75390625}B_{1}}=0,{\color[rgb]{
0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.75390625,0.09765625}C_{1}}=\sigma^{2}_{\star},{\color[rgb]{
0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.75390625,0.09765625,0.09765625}A_{2}}=0,{\color[rgb]{
0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.09765625,0.75390625}B_{2}}=0,{\color[rgb]{
0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{
0.09765625,0.75390625,0.09765625}C_{2}}=0$ (it is the result of Lemma 2). From Theorem 1, choosing any $\alpha>0,$ $\theta=\frac{1}{(1+\gamma\mu)^{2}}$, $\zeta=\frac{\gamma^{2}\sigma^{2}_{\star}}{(1+\gamma\mu)^{2}}$ (see (15) and (16)), we obtain

$$ $\frac{\zeta}{1-\theta}=\frac{\gamma^{2}\sigma^{2}_{\star}}{(1+\gamma\mu)^{2}(1 -\theta)}=\frac{\gamma^{2}\sigma^{2}_{\star}}{(1+\gamma\mu)^{2}-1}=\frac{ \gamma\sigma_{\star}^{2}}{\gamma\mu^{2}+2\mu}.$ $$

For any $k\geq 0$ and any $\gamma>0$, the iterates of SPPM (Algorithm 2) satisfy

$$ $\displaystyle{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\left (\frac{1}{1+\gamma\mu}\right)^{2k}{\left\|x_{0}-x_{\star}\right\|^{2}}+\frac{ \gamma\sigma_{\star}^{2}}{\gamma\mu^{2}+2\mu}.$ $$

∎

###### Lemma 2 ( SPPM ) .

###### Proof of Lemma 2 ..

###### Theorem 2 .

###### Proof of Theorem 2 .

### B.3 Stochastic proximal point method with nonuniform sampling ( SPPM-NS )

Applying SPPM directly to the optimization formulation (3), Theorem 2 applies with $\mu=\min_{i}\mu_{i}$ since $f$ is $\mu$-strongly convex for $\mu=\min_{i}\mu_{i}$. If $\min_{i}\mu_{i}$ is small, the convergence rate becomes weak. In this section we shall describe a trick which enables us to obtain a dependence on the average of the strong convexity constants instead of their minimum.

Choose positive numbers $p_{1},\dots,p_{n}$ summing up to 1, and let

$$ $\tilde{f}_{i}(x):=\frac{1}{np_{i}}f_{i}(x),\qquad i\in[n].$ $$

Note that (3) can be reformulated in the form

$$ $\min_{x\in\mathbb{R}^{d}}\left\{f(x):=\sum_{i=1}^{n}p_{i}\tilde{f}_{i}(x) \right\}.$ (24) $$

We rely on Assumptions 3 and 4.

This is a more refined version of 2, where we assumed the strong convexity parameter was the same for all functions. This implies that $f$ is $\mu$-strongly convex with $\mu=\min_{i}\mu_{i}$. Hence, $f$ has a unique minimizer, which we shall denote $x_{\star}$.

We can now apply SPPM to the reformulated problem (24) instead, with $\mathcal{D}$ being the nonuniform distribution over the finite set $[n]$ given by the parameters $p_{1},\dots,p_{n}$ as follows: $\xi=i$ with probability $p_{i}>0.$ This method is called stochastic proximal point method with nonuniform sampling (SPPM-NS):

$$ $x_{k+1}={\rm prox}_{\gamma\tilde{f}_{i_{k}}}\left(x_{k}\right)={\rm prox}_{ \frac{\gamma}{np_{i}}f_{i_{k}}}\left(x_{k}\right),$ $$

where $i_{k}=i$ with probability $p_{i}>0$.

Figure: Algorithm 3 Stochastic Proximal Point Method with Nonuniform Sampling (SPPM-NS)

Define

$$ $\mu_{\rm NS}:=\min_{i}\frac{\mu_{i}}{np_{i}},\qquad\sigma_{\star,{\rm NS}}^{2} :=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{np_{i}}{\left\|\nabla f_{i}(x_{\star}) \right\|^{2}}.$ (25) $$

###### Lemma 3 ( SPPM-NS ) .

7 holds for the iterates of SPPM-NS (Algorithm 3) with

$$ ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{1}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}=\sigma^{2}_{\star,{\rm NS}},\quad\text {and}\quad{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}=0.$ $$

###### Proof.

Recall that the iterates of (SPPM-NS) have the form

$$ $x_{k+1}={\rm prox}_{\gamma\tilde{f}_{i_{k}}}\left(x_{k}\right)={\rm prox}_{ \frac{\gamma}{np_{i}}f_{i_{k}}}\left(x_{k}\right),$ $$

where $i_{k}=i$ with probability $p_{i}>0.$ Therefore, that $h_{k}=0$. Let $\sigma^{2}\equiv 0$. Clearly, (7) holds. Furthermore,

$$ ${\rm E}\left[\left.{\left\|h_{k}-\nabla\tilde{f}_{i_{k}}(x_{\star})\right\|^{2 }}\right|x_{k},\sigma_{k}^{2}\right]={\rm E}\left[{\left\|\nabla\tilde{f}_{i_{ k}}(x_{\star})\right\|^{2}}\right]:=\sigma^{2}_{\star,{\rm NS}}.$ $$

Hence, 7 holds with

$$ ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{1}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}=\sigma^{2}_{\star,{\rm NS}},\quad\text {and}\quad{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}=0.$ $$

∎

###### Theorem 3 .

Let 3 and 4 hold. Let $x_{0}\in\mathbb{R}^{d}$ be an arbitrary starting point. Let $\mu_{\rm NS}$ and $\sigma^{2}_{\star,{\rm NS}}$ be as in (25). Then for any $k\geq 0$ and any $\gamma>0$, the iterates of SPPM-NS (Algorithm 3) satisfy

$$ $\displaystyle{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\left (\frac{1}{1+\gamma\mu_{\rm NS}}\right)^{2k}{\left\|x_{0}-x_{\star}\right\|^{2} }+\frac{\gamma\sigma_{\star,{\rm NS}}^{2}}{\gamma\mu_{\rm NS}^{2}+2\mu_{\rm NS }}.$ (26) $$

Commentary:

- (a)
Uniform sampling. If we choose $p_{i}=\frac{1}{n}$ for all $i\in[n],$ we shall refer to Algorithm 3 as stochastic proximal point method with uniform sampling (SPPM-US). In this case,
$\mu_{\rm NS}=\mu_{\rm US}:=\min_{i}\mu_{i},\qquad\sigma_{\star,{\rm NS}}^{2}=
\sigma_{\star,{\rm US}}^{2}:=\frac{1}{n}\sum_{i=1}^{n}{\left\|\nabla f_{i}(x_{
\star})\right\|^{2}}.$
- (b)
Importance sampling: optimizing the linear rate. If we choose $p_{i}=\frac{\mu_{i}}{\sum_{j=1}^{n}\mu_{j}}$ for all $i\in[n],$ we shall refer to Algorithm 3 as stochastic proximal point method with importance sampling (SPPM-IS). In this case,
$\mu_{\rm NS}=\mu_{\rm IS}:=\frac{1}{n}\sum_{i=1}^{n}\mu_{i},\qquad\sigma_{
\star,{\rm NS}}^{2}=\sigma_{\star,{\rm IS}}^{2}:=\frac{\sum_{i=1}^{n}\mu_{i}}{
n}\sum_{i=1}^{n}\frac{{\left\|\nabla f_{i}(x_{\star})\right\|^{2}}}{n\mu_{i}}.$
This choice maximizes the value of $\mu_{\rm NS}$ (and hence the first part of the convergence rate) over the choice of the probabilities.
- (c)
Variance sampling: optimizing the variance.
If we choose $p_{i}=\frac{{\left\|\nabla f_{i}(x_{\star})\right\|}}{\sum_{j=1}^{n}{\left\|
\nabla f_{j}(x_{\star})\right\|}}$ for all $i\in[n]$, we shall refer to Algorithm 3 as stochastic proximal point method with variance sampling (SPPM-VS). In this case,
$\mu_{\rm NS}=\mu_{\rm VS}:=\frac{1}{n}\sum_{i=1}^{n}{\left\|\nabla f_{i}(x_{
\star})\right\|}\left(\min_{i}\frac{\mu_{i}}{{\left\|\nabla f_{i}(x_{\star})
\right\|}}\right),$
$\sigma_{\star,{\rm NS}}^{2}=\sigma_{\star,{\rm VS}}^{2}:=\left(\frac{1}{n}\sum
_{i=1}^{n}{\left\|\nabla f_{i}(x_{\star})\right\|}\right)^{2}.$
This choice minimizes the value of $\sigma_{\star,{\rm NS}}$ (and hence the second part of the convergence rate) over the choice of the probabilities.

###### Proof of Theorem 3 .

From Lemma 3, we have that 7 holds for the iterates of SPPM-NS (Algorithm 3) with

$$ ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{1}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}=\sigma^{2}_{\star,{\rm NS}},\quad\text {and}\quad{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}=0.$ $$

From Theorem 1, choosing any $\alpha>0,$ $\theta=\frac{1}{\left(1+\gamma\mu_{\rm NS}\right)^{2}},$ $\zeta=\frac{\gamma^{2}\sigma^{2}_{\star,{\rm NS}}}{\left(1+\gamma\mu_{{\rm NS}
}\right)^{2}}$ (see (15) and (16)), we obtain that, for any $k\geq 0$ and any $\gamma>0$, the iterates of SPPM-NS (Algorithm 3) satisfy

$$ $\displaystyle{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\left (\frac{1}{1+\gamma\mu_{\rm NS}}\right)^{2k}{\left\|x_{0}-x_{\star}\right\|^{2} }+\frac{\gamma\sigma_{\star,{\rm NS}}^{2}}{\gamma\mu_{\rm NS}^{2}+2\mu_{\rm NS }}.$ $$

Alternatively, we can derive Theorem 3 from Theorem 2. It is easy to show that each function $\tilde{f}_{i}$ is $\tilde{\mu}_{i}$-strongly convex with $\tilde{\mu}_{i}:=\frac{\mu_{i}}{np_{i}}$. Hence, $\tilde{f}_{i}$ is also $\mu_{\rm NS}$-strongly convex for all $i$. It now only remains to apply Theorem 2.
∎

###### Lemma 3 ( SPPM-NS ) .

###### Proof.

###### Theorem 3 .

###### Proof of Theorem 3 .

### B.4 Stochastic proximal point method with arbitrary sampling ( SPPM-AS )

We consider the optimization problem (3)
and rely on 3 (differentiability) and 4 (strong convexity).

Let $\mathcal{S}$ be a probability distribution over the $2^{n}$ subsets of $[n]$. Given a random set $S\sim\mathcal{S}$, we define

$$ $p_{i}:=\mathbf{Prob}\left(i\in S\right),\qquad i\in[n].$ (27) $$

We will restrict our attention to proper and nonvacuous random sets.

###### Assumption 8 .

$S$ is proper (i.e., $p_{i}>0$ for all $i\in[n]$) and nonvacuous (i.e., $\mathbf{Prob}\left(S=\varnothing\right)=0$).

Given $\varnothing\neq C\subseteq[n]$ and $i\in[n]$, we define

$$ $v_{i}(C):=\begin{cases}\frac{1}{p_{i}}&i\in C\\ 0&i\notin C\end{cases},$ (28) $$

and

$$ $f_{C}(x):=\frac{1}{n}\sum_{i=1}^{n}v_{i}(C)f_{i}(x)\overset{\eqref{eq:v_i(S)-9 0u}}{=}\sum_{i\in C}\frac{1}{np_{i}}f_{i}(x).$ (29) $$

Note that $v_{i}(S)$ is a random variable and $f_{S}$ is a random function. By construction, ${\rm E}_{S\sim\mathcal{S}}\left[v_{i}(S)\right]=1$ for all $i\in[n]$, and hence

$$ ${\rm E}_{S\sim\mathcal{S}}\left[f_{S}(x)\right]={\rm E}_{S\sim\mathcal{S}} \left[\frac{1}{n}\sum_{i=1}^{n}v_{i}(S)f_{i}(x)\right]=\frac{1}{n}\sum_{i=1}^{ n}{\rm E}_{S\sim\mathcal{S}}\left[v_{i}(S)\right]f_{i}(x)=\frac{1}{n}\sum_{i=1 }^{n}f_{i}(x)=f(x).$ $$

Therefore, the optimization problem (3) is equivalent to the stochastic optimization problem

$$ $\min_{x\in\mathbb{R}^{d}}\left\{f(x):={\rm E}_{S\sim\mathcal{S}}\left[f_{S}(x) \right]\right\}.$ (30) $$

Further, if for each $C\subset[n]$ we let $p_{C}:=\mathbf{Prob}\left(S=C\right)$, $f$ can be written in the equivalent form

$$ $f(x)={\rm E}_{S\sim\mathcal{S}}\left[f_{S}(x)\right]=\sum_{C\subseteq[n]}p_{C} f_{C}(x)=\sum_{C\subseteq[n],p_{C}>0}p_{C}f_{C}(x).$ (31) $$

Applying SPPM to (30), we arrive at
stochastic proximal point method with arbitrary sampling (SPPM-AS) (Algorithm 4):

$$ $x_{k+1}={\rm prox}_{\gamma f_{S_{k}}}\left(x_{k}\right),$ $$

where $S_{k}\sim\mathcal{S}$.

Figure: Algorithm 4 Stochastic Proximal Point Method with Arbitrary Sampling (SPPM-AS)

Define

$$ $\mu_{\rm AS}:=\min_{C\subseteq[n],p_{C}>0}\sum_{i\in C}\frac{\mu_{i}}{np_{i}}, \qquad\sigma_{\star,{\rm AS}}^{2}:=\sum_{C\subseteq[n],p_{C}>0}p_{C}{\left\| \sum_{i\in C}\frac{1}{np_{i}}\nabla f_{i}(x_{\star})\right\|^{2}}.$ (32) $$

###### Lemma 4 .

7 holds for the iterates of SPPM-AS (Algorithm 4) with

$$ ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{1}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}=\sigma^{2}_{\star,{\rm AS}},\quad\text {and}\quad{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}=0.$ $$

###### Proof.

Recall that the iterates of SPPM-AS (Algorithm 4) have the form

$$ $x_{k+1}={\rm prox}_{\gamma f_{S_{k}}}\left(x_{k}\right).$ $$

Thus, $h_{k}=0$. Let $\sigma^{2}\equiv 0$. Clearly, (7) holds. Furthermore,

$$ ${\rm E}\left[\left.{\left\|h_{k}-\nabla f_{C}(x_{\star})\right\|^{2}}\right|x_ {k},\sigma_{k}^{2}\right]={\rm E}\left[{\left\|\nabla f_{C}(x_{\star})\right\| ^{2}}\right]:=\sigma^{2}_{\star,{\rm AS}},$ $$

since

$$ $\displaystyle\sigma_{\star}^{2}$ $\displaystyle:=$ $\displaystyle{\rm E}_{\xi\sim\mathcal{D}}\left[{\left\|\nabla f_{\xi}(x_{\star })\right\|^{2}}\right]$ $\displaystyle\overset{\eqref{eq:SPPM-AS-reform2}}{=}$ $\displaystyle\sum_{C\subseteq[n],p_{C}>0}p_{C}{\left\|\nabla f_{C}(x_{\star}) \right\|^{2}}$ $\displaystyle\overset{\eqref{eq:f_C}}{=}$ $\displaystyle\sum_{C\subseteq[n],p_{C}>0}p_{C}{\left\|\sum_{i\in C}\frac{1}{np _{i}}\nabla f_{i}(x_{\star})\right\|^{2}}$ $\displaystyle:=$ $\displaystyle\sigma_{\star,{\rm AS}}^{2}.$ $$

Hence, 7 holds with

$$ ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{1}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}=\sigma^{2}_{\star,{\rm AS}},\quad\text {and}\quad{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}=0.$ $$

∎

The convergence of SPPM-AS is captured by the following theorem.

###### Theorem 4 .

Let 3 (differentiability) and 4 (strong convexity) hold. Let $S$ be a random set satisfying 8. Let $x_{0}\in\mathbb{R}^{d}$ be an arbitrary starting point, $\mu_{\rm AS}$ and $\sigma_{\star,{\rm AS}}^{2}$ be as in (32). Then for any $k\geq 0$ and any $\gamma>0$, the iterates of SPPM-AS (Algorithm 4) satisfy

$$ $\displaystyle{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\left (\frac{1}{1+\gamma\mu_{\rm AS}}\right)^{2k}{\left\|x_{0}-x_{\star}\right\|^{2} }+\frac{\gamma\sigma_{\star,{\rm AS}}^{2}}{\gamma\mu_{\rm AS}^{2}+2\mu_{\rm AS }}.$ (33) $$

Commentary:

- (a)
Full sampling. Let $S=[n]$ with probability 1 (“full sampling”; abbreviated as “FS”). Then SPPM-AS applied to (30) becomes PPM for minimizing $f$. So, besides Theorem 2, Theorem 4 is the third theorem capturing the convergence of PPM in the differentiable and strongly convex regime.
Moreover, we have $p_{i}=1$ for all $i\in[n]$ and the expressions (32) take on the form
$\mu_{\rm AS}=\mu_{\rm FS}:=\frac{1}{n}\sum_{i=1}^{n}\mu_{i},\qquad\sigma_{
\star,{\rm AS}}^{2}=\sigma_{\star,{\rm FS}}^{2}:=0.$
Note that $\mu_{\rm FS}$ a lower bound on the the strong convexity constant of $f$ – one that can be computed from the strong convexity constants $\mu_{i}$ of the constituent functions $f_{i}$.
- (b)
Nonuniform sampling. Let $S=\{i\}$ with probability $q_{i}>0$, where $\sum_{i}q_{i}=1$. This leads to SPPM-NS. Then $p_{i}:=\mathbf{Prob}\left(i\in S\right)=q_{i}$ for all $i\in[n]$, and the expressions (32) take on the form
$\mu_{\rm AS}=\mu_{\rm NS}:=\min_{i}\frac{\mu_{i}}{np_{i}},\qquad\sigma_{\star,
{\rm AS}}^{2}=\sigma_{\star,{\rm NS}}^{2}:=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{
np_{i}}{\left\|\nabla f_{i}(x_{\star})\right\|^{2}},$
which recovers the quantities from Section B.3; see (25).
- (c)
Nice sampling. Choose $\tau\in[n]$ and let $S$ be a random subset of $[n]$ of size $\tau$ chosen uniformly at random. We call the resulting method SPPM-NICE. Then $p_{i}:=\mathbf{Prob}\left(i\in S\right)=\frac{\tau}{n}$ for all $i\in[n]$. Moreover, $p_{C}=\frac{1}{\binom{n}{\tau}}$ whenever $|C|=\tau$ and $p_{C}=0$ otherwise. So,
$\mu_{\rm AS}=\mu_{\rm NICE}(\tau):=\min_{C\subseteq[n],|C|=\tau}\frac{1}{\tau}
\sum_{i\in C}\mu_{i}$
(34)
and
$\sigma_{\star,{\rm AS}}^{2}=\sigma_{\star,{\rm NICE}}^{2}(\tau):=\sum_{C
\subseteq[n],|C|=\tau}\frac{1}{\binom{n}{\tau}}{\left\|\frac{1}{\tau}\sum_{i
\in C}\nabla f_{i}(x_{\star})\right\|^{2}}.$
(35)
It can be shown that $\mu_{\rm NICE}(\tau)$ is a nondecreasing function of $\tau$. So, as the minibatch size $\tau$ increases, the strong convexity constant $\mu_{\rm NICE}(\tau)$ can only improve. Since $\mu_{\rm NICE}(1)=\min_{i}\mu_{i}$ and $\mu_{\rm NICE}(n)=\frac{1}{n}\sum_{i=1}^{n}\mu_{i}$, the value of $\mu_{\rm NICE}(\tau)$ interpolates these two extreme cases as $\tau$ varies between $1$ and $n$.
- (d)
Block sampling. Let $C_{1},\dots,C_{b}$ be a partition of $[n]$ into $b$ nonempty blocks. For each $i\in[n]$, let $B(i)$ indicate which block does $i$ belong to. In other words, $i\in C_{j}$ iff $B(i)=j$. Let $S=C_{j}$ with probability $q_{j}>0$, where $\sum_{j}q_{j}=1$. We call the resulting method SPPM-BS. Then $p_{i}:=\mathbf{Prob}\left(i\in S\right)=q_{B(i)},$ and hence the expressions (32) take on the form
$\mu_{\rm AS}=\mu_{\rm BS}:=\min_{j\in[b]}\frac{1}{nq_{j}}\sum_{i\in C_{j}}\mu_
{i},\qquad\sigma_{\star,{\rm AS}}^{2}=\sigma_{\star,{\rm BS}}^{2}:=\sum_{j\in[
b]}q_{j}{\left\|\sum_{i\in C_{j}}\frac{1}{np_{j}}\nabla f_{i}(x_{\star})\right
\|^{2}}.$
We now consider two extreme cases:
–
If $b=1$, then SPPM-BS = SPPM-FS = PPM. Let’s see, as a sanity check, whether we recover the right rate as well. We have $q_{1}=1$, $C_{1}=[n]$, $p_{i}:=\mathbf{Prob}\left(i\in S\right)=1$ for all $i\in[n]$, and the expressions for $\mu_{\rm AS}$ and $\sigma_{\star,{\rm BS}}^{2}$ simplify to
$\mu_{\rm BS}=\mu_{\rm FS}:=\frac{1}{n}\sum_{i=1}^{n}\mu_{i},\qquad\sigma_{
\star,{\rm BS}}^{2}=\sigma_{\star,{\rm FS}}^{2}:=0.$
So, indeed, we recover the same rate as SPPM-FS.
–
If $b=n$, then SPPM-BS = SPPM-NS. Let’s see, as a sanity check, whether we recover the right rate as well. We have $C_{i}=\{i\}$ and $p_{i}:=\mathbf{Prob}\left(i\in S\right)=q_{i}$ for all $i\in[n]$, and the expressions for $\mu_{\rm AS}$ and $\sigma_{\star,{\rm BS}}^{2}$ simplify to
$\mu_{\rm BS}=\mu_{\rm NS}:=\min_{i\in[n]}\frac{\mu_{i}}{np_{i}},\qquad\sigma_{
\star,{\rm BS}}^{2}=\sigma_{\star,{\rm NS}}^{2}:=\frac{1}{n}\sum_{i=1}^{n}
\frac{1}{np_{i}}{\left\|\nabla f_{i}(x_{\star})\right\|^{2}}.$
So, indeed, we recover the same rate as SPPM-NS.
- (e)
Stratified sampling. Let $C_{1},\dots,C_{b}$ be a partition of $[n]$ into $b$ nonempty blocks, as before. For each $i\in[n]$, let $B(i)$ indicate which block $i$ belongs to. In other words, $i\in C_{j}$ iff $B(i)=j$. Now, for each $j\in[b]$ pick $\xi_{j}\in C_{j}$ uniformly at random, and define $S=\cup_{j\in[b]}\{\xi_{j}\}$. We call the resulting method SPPM-SS. Clearly, $p_{i}:=\mathbf{Prob}\left(i\in S\right)=\frac{1}{|C_{B(i)}|}$. The expressions (32) take on the form
$\mu_{\rm AS}=\mu_{\rm SS}:=\min_{(i_{1},\dots,i_{b})\in C_{1}\times\cdots
\times C_{b}}\sum_{j=1}^{b}\frac{\mu_{i_{j}}|C_{j}|}{n}$
and
$\sigma_{\star,{\rm AS}}^{2}=\sigma_{\star,{\rm SS}}^{2}:=\left(\frac{1}{\prod_
{j=1}^{b}|C_{j}|}\right)\sum_{(i_{1},\dots,i_{b})\in C_{1}\times\cdots\times C
_{b}}{\left\|\sum_{j=1}^{b}\frac{|C_{j}|}{n}\nabla f_{i_{j}}(x_{\star})\right
\|^{2}}.$
We now consider two extreme cases:
–
If $b=1$, then SPPM-SS = SPPM-US. Let’s see, as a sanity check, whether we recover the right rate as well. We have $C_{1}=[n]$, $|C_{1}|=n$, $\left(\prod_{j=1}^{b}\frac{1}{|C_{j}|}\right)=\frac{1}{n}$ and hence
$\mu_{\rm SS}=\mu_{\rm US}:=\min_{i}\mu_{i},\qquad\sigma_{\star,{\rm SS}}^{2}=
\sigma_{\star,{\rm US}}^{2}:=\frac{1}{n}\sum_{i=1}^{n}{\left\|\nabla f_{i}(x_{
\star})\right\|^{2}}.$
So, indeed, we recover the same rate as SPPM-US.
–
If $b=n$, then SPPM-SS = SPPM-FS. Let’s see, as a sanity check, whether we recover the right rate as well. We have $C_{i}=\{i\}$ for all $i\in[n]$, $\left(\prod_{j=1}^{b}\frac{1}{|C_{j}|}\right)=1$, and hence
$\mu_{\rm SS}=\mu_{\rm FS}:=\frac{1}{n}\sum_{i=1}^{n}\mu_{i},\qquad\sigma_{
\star,{\rm SS}}^{2}=\sigma_{\star,{\rm FS}}^{2}:=0.$
So, indeed, we recover the same rate as SPPM-FS.

###### Proof of Theorem 4 .

From Lemma 4 we have that 7 holds for the iterates of SPPM-AS (Algorithm 4) with

$$ ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{1}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}=\sigma^{2}_{\star,{\rm AS}},\quad\text {and}\quad{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}=0.$ $$

From Theorem 1, choosing $\alpha>0,$ $\theta=\frac{1}{(1+\gamma\mu)^{2}},$ $\zeta=\frac{\gamma^{2}\sigma_{*,{\rm AS}}^{2}}{(1+\gamma\mu)^{2}}$ (see (15) and (16)), we obtain that, for any $k\geq 0$ and any $\gamma>0$, the iterates of SPPM-AS (Algorithm 4) satisfy

$$ $\displaystyle{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\left (\frac{1}{1+\gamma\mu_{\rm AS}}\right)^{2k}{\left\|x_{0}-x_{\star}\right\|^{2} }+\frac{\gamma\sigma_{\star,{\rm AS}}^{2}}{\gamma\mu_{\rm AS}^{2}+2\mu_{\rm AS }}.$ $$

Alternatively, we can derive Theorem 4 from Theorem 2. Let $C$ be any (necessarily nonempty) subset of $[n]$ such that $p_{C}>0$. Recall that in view of (29) we have

$$ $f_{C}(x)=\sum_{i\in C}\frac{1}{np_{i}}f_{i}(x);$ $$

i.e., $f_{C}$ is a conic combination of the functions $\{f_{i}\;:\;i\in C\}$ with weights $w_{i}=\frac{1}{np_{i}}$. Since each $f_{i}$ is $\mu_{i}$-strongly convex, Lemma 10 says that $f_{C}$
$\mu_{C}$-strongly convex with

$$ $\mu_{C}:=\sum_{i\in C}\frac{\mu_{i}}{np_{i}}.$ $$

So, every such $f_{C}$ is $\mu$-strongly convex with

$$ $\mu=\mu_{\rm AS}:=\min_{C\subseteq[n],p_{C}>0}\sum_{i\in C}\frac{\mu_{i}}{np_{ i}}.$ $$

Further, the quantity $\sigma_{\star}^{2}$ from (20) is equal to

$$ $\displaystyle\sigma_{\star}^{2}$ $\displaystyle:=$ $\displaystyle{\rm E}_{\xi\sim\mathcal{D}}\left[{\left\|\nabla f_{\xi}(x_{\star })\right\|^{2}}\right]$ $\displaystyle\overset{\eqref{eq:SPPM-AS-reform2}}{=}$ $\displaystyle\sum_{C\subseteq[n],p_{C}>0}p_{C}{\left\|\nabla f_{C}(x_{\star}) \right\|^{2}}$ $\displaystyle\overset{\eqref{eq:f_C}}{=}$ $\displaystyle\sum_{C\subseteq[n],p_{C}>0}p_{C}{\left\|\sum_{i\in C}\frac{1}{np _{i}}\nabla f_{i}(x_{\star})\right\|^{2}}$ $\displaystyle:=$ $\displaystyle\sigma_{\star,{\rm AS}}^{2}.$ $$

It now only remains to apply Theorem 2.
∎

###### Assumption 8 .

###### Lemma 4 .

###### Proof.

###### Theorem 4 .

###### Proof of Theorem 4 .

### B.5 Stochastic proximal point method with optimal gradient correction ( SPPM* )

We showed that SPPM converges up to a neighborhood of size $\frac{\gamma\sigma_{\star}^{2}}{\gamma\mu^{2}+2\mu}$, where

$$ $\sigma_{\star}^{2}:={\rm E}_{\xi\sim\mathcal{D}}\left[{\left\|\nabla f_{\xi}(x _{\star})\right\|^{2}}\right].$ $$

We now describe a simple trick which gets rid of the neighborhood when $\sigma_{\star}^{2}>0$. The trick is of a conceptual nature: as is, it is practically useless. However, it will serve as an inspiration for a trick that can be implemented.

Figure: Algorithm 5 Stochastic Proximal Point Method with Optimal Gradient Correction (SPPM*)

Let us reformulate the problem by adding a smart zero.

For each $\xi\sim\mathcal{D}$, define

$$ $\tilde{f}_{\xi}(x):=f_{\xi}(x)-\left\langle\nabla f_{\xi}(x_{\star}),x\right\rangle,$ (36) $$

and instead of solving (1), consider solving the problem

$$ $\min_{x\in\mathbb{R}^{d}}\left\{\tilde{f}(x):={\rm E}_{\xi\sim\mathcal{D}} \left[\tilde{f}_{\xi}(x)\right]\right\}.$ (37) $$

Observations:

- •
We do not know $\nabla f_{\xi}(x_{\star})$, and hence formulation (37) is not of practical interest.
- •
It is easy to see that $f=\tilde{f}$, and hence problems (1) and (37) are equivalent. Indeed,
$\displaystyle\tilde{f}(x)\quad\overset{\eqref{eq:opt-problem-SPPM-shift}}{=}
\quad{\rm E}_{\xi\sim\mathcal{D}}\left[\tilde{f}_{\xi}(x)\right]$
$\displaystyle\overset{\eqref{eq:smart_zero_STAR}}{=}$
$\displaystyle{\rm E}_{\xi\sim\mathcal{D}}\left[f_{\xi}(x)-\left\langle\nabla f
_{\xi}(x_{\star}),x\right\rangle\right]$
$\displaystyle=$
$\displaystyle{\rm E}_{\xi\sim\mathcal{D}}\left[f_{\xi}(x)\right]-{\rm E}_{\xi
\sim\mathcal{D}}\left[\left\langle\nabla f_{\xi}(x_{\star}),x\right\rangle\right]$
$\displaystyle=$
$\displaystyle{\rm E}_{\xi\sim\mathcal{D}}\left[f_{\xi}(x)\right]-\left\langle
\underbrace{{\rm E}_{\xi\sim\mathcal{D}}\left[\nabla f_{\xi}(x_{\star})\right]
}_{=\nabla f(x_{\star})=0},x\right\rangle\quad\overset{\eqref{eq:main-opt-
problem}}{=}\quad f(x).$
- •
All stochastic gradients of $\tilde{f}$ at $x_{\star}$ are zero. Indeed,
$\nabla\tilde{f}(x)\overset{\eqref{eq:smart_zero_STAR}}{=}\nabla f_{\xi}(x)-
\nabla f_{\xi}(x_{\star}),$
and hence $\nabla\tilde{f}_{\xi}(x_{\star})=0$.
- •
It is easy to see that since $f_{\xi}$ is differentiable and $\mu$-strongly convex, then so is $\tilde{f}_{\xi}$.

Hiding the prox. Further, recall that if for some $x\in\mathbb{R}^{d}$ and differentiable and convex $\phi$ we let $x_{+}:={\rm prox}_{\phi}\left(x\right)$, then $x_{+}=x-\nabla\phi(x_{+})$. Therefore, steps 4 and 5 of the method can be written in the equivalent form

$$ $x_{k+1}=x_{k}+\gamma{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor [named]{pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}h_{k}}-\gamma \nabla f_{\xi_{k}}(x_{k+1})=x_{k}-\gamma\left(\nabla f_{\xi_{k}}(x_{k+1})-{ \color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}\nabla f_{\xi_{k}}(x_{ \star})}\right).$ $$

###### Lemma 5 ( SPPM* ) .

7 holds for the iterates of SPPM* (Algorithm 5) with

$$ ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{1}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}=0,\quad\text{and}\quad{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{2}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}=0.$ $$

###### Proof.

Recall that the iterates of SPPM* have the form

$$ $x_{k+1}={\rm prox}_{\gamma f_{\xi_{k}}}\left(x_{k}+\gamma\nabla f_{\xi_{k}}(x_ {\star})\right).$ $$

Thus, $h_{k}=\nabla f_{\xi_{k}}(x_{\star})$. Let $\sigma^{2}\equiv 0$. Clearly, (7) holds. Furthermore,

$$ ${\rm E}\left[\left.{\left\|h_{k}-\nabla f_{\xi_{k}}(x_{\star})\right\|^{2}} \right|x_{k},\sigma_{k}^{2}\right]=0.$ $$

Hence, 7 holds with

$$ ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{1}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}=0,\quad\text{and}\quad{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{2}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}=0.$ $$

∎

The convergence of SPPM* is captured by the following theorem.

###### Theorem 5 .

Let 1 (differentiability), 2 ($\mu$-strong convexity) hold. Let $x_{0}\in\mathbb{R}^{d}$ be an arbitrary starting point. Then for any $k\geq 0$ and any $\gamma>0$, the iterates of SPPM* (Algorithm 5) satisfy

$$ $\displaystyle{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\left (\frac{1}{1+\gamma\mu}\right)^{2k}{\left\|x_{0}-x_{\star}\right\|^{2}}.$ (38) $$

Commentary:

- •
The convergence neighborhood is fully removed; the method converges to the exact solution! The result is identical to (22); i.e., to the rate of SPPM in the interpolation regime.
–
The method converges to $x_{\star}$ for any fixed $\gamma>0$ as long as $k\to\infty$.
–
The method converges to $x_{\star}$ for any fixed $k\geq 1$ as long as $\gamma\to\infty$.
- •
The method is practically useless since it relies on the knowledge of the optimal stochastic gradients $\nabla f_{\xi}(x_{\star})$ for all $\xi$ as hyper-parameters of the method. Needless to say, these vectors are rarely known.
- •
If we follow (36), add “smart zero” to (30), and apply SPPM to the resulting formulation, we automatically get a “star” variant of SPPM-AS, and Theorem 5 captures the complexity of this method.

###### Proof of Theorem 5 .

From Lemma 5 we know that 7 holds for the iterates of SPPM* (Algorithm 5) with

$$ ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{1}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}=0,\quad\text{and}\quad{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{2}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}=0.$ $$

Therefore, from Theorem 1, by choosing any $\alpha>0,$ $\theta=\left(\frac{1}{1+\gamma\mu}\right)^{2},$ $\zeta=0$ (see (15) and (16)), we have

$$ $\displaystyle{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\left (\frac{1}{1+\gamma\mu}\right)^{2k}{\left\|x_{0}-x_{\star}\right\|^{2}}.$ $$

∎

###### Lemma 5 ( SPPM* ) .

###### Proof.

###### Theorem 5 .

###### Proof of Theorem 5 .

### B.6 Stochastic proximal point method with gradient correction ( SPPM-GC )

We consider the stochastic optimization problem (1), i.e.,

$$ $\min_{x\in\mathbb{R}^{d}}\left\{f(x):={\rm E}_{\xi\sim\mathcal{D}}\left[f_{\xi }(x)\right]\right\},$ $$

and rely on 1 (differentiability of $f_{\xi}$) and 2 ($\mu$-strong convexity of $f_{\xi}$). Recall that this implies strong convexity of $f.$ Hence $f$ has a unique minimizer, which we denote $x_{\star}$.

We have already described the SPPM* method – this a variant of SPPM without the neighborhood term in the convergence bound. In order to run it, we need to know $\nabla f_{\xi}(x_{\star})$ for all $\xi$, which is of course something we almost never know; one exception to this is the interpolation regime, defined by the assumption that $\nabla f_{\xi}(x_{\star})=0$ for all $\xi$.

In this section we describe a practical method inspired by SPPM*. The method is based on the following ideas:

- •
While we do not know $\nabla f_{\xi}(x_{\star})$, what if this quantity could be in some appropriate/useful sense approximated by some vector $g_{\xi}$ we do know?
- •
One option is to require a quantity such as
${\left\|g_{\xi}-\nabla f_{\xi}(x_{\star})\right\|^{2}}\qquad\text{or}\qquad{
\rm E}_{\xi\sim\mathcal{D}}\left[{\left\|g_{\xi}-\nabla f_{\xi}(x_{\star})
\right\|^{2}}\right]$
to be “small” in some sense. However, where can such vectors come from? And what should “small” mean?
- •
One idea is to make an extra assumption on the functions $f_{\xi}$ that would somehow automatically guarantee the existence and availability of such vectors. So, we trade off easy availability of vectors $g_{\xi}$ for a limitation on the class of problems we solve this way.

Recall the Similarity 5.
There exists $\delta\geq 0$ such that

$$ ${\rm E}_{\xi\sim\mathcal{D}}\left[{\left\|\nabla f_{\xi}(x)-\nabla f(x)-\nabla f _{\xi}(x_{\star})\right\|^{2}}\right]\leq\delta^{2}{\left\|x-x_{\star}\right\| ^{2}},\qquad\forall x\in\mathbb{R}^{d}.$ $$

The above assumption says that we can consider the vectors $g_{\xi}=\nabla f_{\xi}(x)-\nabla f(x)$ for any $x$, and that by “small enough” we require a bound by $\delta^{2}{\left\|x-x_{\star}\right\|^{2}}$. Why these particular choices make sense will become clear from the convergence proof. Let us now make some observations about the class of functions satisfying 5:

- 1.
The approximation of $\nabla f_{\xi}(x_{\star})$ by $\nabla f_{\xi}(x)-\nabla f(x)$ gets better as $x$ gets closer to $x_{\star}$.
- 2.
The smaller $\delta$ is, the better the approximation.
- 3.
If $\delta=0$, we get perfect approximation.
- 4.
Since
${\rm E}_{\xi\sim\mathcal{D}}\left[\nabla f_{\xi}(x)-\nabla f_{\xi}(x_{\star})
\right]=\nabla f(x),$
the left-hand side of (5) is the variance of the random vector $\nabla f_{\xi}(x)-\nabla f_{\xi}(x_{\star})$ as an estimator of $\nabla f(x)=\nabla f(x)-\nabla f(x_{\star})$.
- 5.
It follows that (5) can be equivalently written in the form
${\rm E}_{\xi\sim\mathcal{D}}\left[{\left\|\nabla f_{\xi}(x)-\nabla f_{\xi}(x_{
\star})\right\|^{2}}\right]-{\left\|\nabla f(x)-\nabla f(x_{\star})\right\|^{2
}}\leq\delta^{2}{\left\|x-x_{\star}\right\|^{2}},\qquad\forall x\in\mathbb{R}^
{d}.$
(39)
- 6.
Note that (39) holds if the following stronger condition holds:
${\rm E}_{\xi\sim\mathcal{D}}\left[{\left\|\nabla f_{\xi}(x)-\nabla f_{\xi}(x_{
\star})\right\|^{2}}\right]\leq\delta^{2}{\left\|x-x_{\star}\right\|^{2}},
\qquad\forall x\in\mathbb{R}^{d}.$
(40)
- 7.
Furthermore, (40) holds if the following even stronger condition holds: there exists $\delta\geq 0$ such that
${\left\|\nabla f_{\xi}(x)-\nabla f_{\xi}(x_{\star})\right\|}\leq\delta{\left\|
x-x_{\star}\right\|},\qquad\forall x\in\mathbb{R}^{d}$
(41)
for all $\xi$.
- 8.
Finally, (41) holds if there exists $\delta\geq 0$ such that $\nabla f_{\xi}$ is $\delta$-Lipschitz for all $\xi$:
${\left\|\nabla f_{\xi}(x)-\nabla f_{\xi}(y)\right\|}\leq\delta{\left\|x-y
\right\|},\qquad\forall x,y\in\mathbb{R}^{d}.$
(42)

For each $\xi\in[n]$ and any “parameter” $v\in\mathbb{R}^{d}$, define

$$ $\tilde{f}_{\xi}^{v}(x):=f_{\xi}(x)-\left\langle\nabla f_{\xi}(v)-\nabla f(v),x \right\rangle,$ (43) $$

and instead of solving (1), consider solving the problem

$$ $\min_{x\in\mathbb{R}^{d}}\left\{\tilde{f}^{v}(x):={\rm E}_{\xi\sim\mathcal{D}} \left[\tilde{f}_{\xi}^{v}(x)\right]\right\}.$ (44) $$

Figure: Algorithm 6 Stochastic Proximal Point Method with Gradient Correction (SPPM-GC)

The algorithm can be interpreted as:

- •
a practical variant of SPPM* in which we use the computable correction $\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}h_{k}=\nabla f_{\xi_{k}}
(x_{k})-\nabla f(x_{k})$ instead of incomputable correction $\nabla f_{\xi_{k}}(x_{\star})$,
- •
SPPM applied to the reformulated problem (44), with the “control” vector $v=x_{k}$ at iteration $k$:
$v=x_{k},\qquad x_{k+1}={\rm prox}_{\gamma\tilde{f}_{\xi_{k}}^{v}}\left(x_{k}\right)$

Hiding the prox. Further, recall that if for some $x\in\mathbb{R}^{d}$ and differentiable and convex $\phi$ we let $x_{+}:={\rm prox}_{\phi}\left(x\right)$, then $x_{+}=x-\nabla\phi(x_{+})$. Therefore, steps 4 and 5 of the method can be written in the equivalent form

$$ $x_{k+1}=x_{k}+\gamma{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor [named]{pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}h_{k}}-\gamma \nabla f_{\xi_{k}}(x_{k+1})=x_{k}-\gamma\left(\nabla f_{\xi_{k}}(x_{k+1})-{ \color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}\nabla f_{\xi_{k}}(x_{k} )+\nabla f(x_{k})}\right).$ $$

###### Lemma 6 ( SPPM-GC ) .

Suppose Assumption 5 holds with $\delta>0.$ 7 holds for the iterates of SPPM-GC (Algorithm 6) with

$$ ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=\delta^{2},{ \color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}=0,\quad\text{and}\quad{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{2}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}=0.$ $$

###### Proof.

Recall that the iterates of SPPM-GC have the form

$$ $x_{k+1}={\rm prox}_{\gamma f_{\xi_{k}}}\left(x_{k}+\gamma h_{k}\right),$ $$

where $h_{k}=\nabla f_{\xi_{k}}(x_{k})-\nabla f(x_{k})$. Let $\sigma^{2}\equiv 0$. Since Assumption 5 holds, we get

$$ ${\rm E}\left[\left.{\left\|h_{k}-\nabla f_{\xi_{k}}(x_{\star})\right\|^{2}} \right|x_{k},\sigma_{k}^{2}\right]={\rm E}\left[\left.{\left\|\nabla f_{\xi_{k }}(x_{k})-\nabla f(x_{k})-\nabla f_{\xi_{k}}(x_{\star})\right\|^{2}}\right|x_{ k}\right]\leq\delta^{2}{\left\|x_{k}-x_{\star}\right\|^{2}}.$ $$

Hence, 7 holds with

$$ ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=\delta^{2},{ \color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}=0,\quad\text{and}\quad{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{2}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}=0.$ $$

∎

###### Theorem 6 .

Let 1 (differentiability), 2 ($\mu$-strong convexity), and 5 ($\delta$-similarity) hold. Choose any $x_{0}\in\mathbb{R}^{d}$. Then for any $\gamma>0$, and all $k\geq 0$, we have

$$ $\displaystyle{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\left (\frac{1+\gamma^{2}\delta^{2}}{\left(1+\gamma\mu\right)^{2}}\right)^{k}{\left \|x_{0}-x_{\star}\right\|^{2}}.$ (45) $$

Commentary:

- 1.
Perfect similarity. If $\delta=0$, for all $\gamma>0$ we get
$\displaystyle{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\left
(\frac{1}{1+\gamma\mu}\right)^{2k}{\left\|x_{0}-x_{\star}\right\|^{2}}.$
(46)
We get convergence even with $k=1$ provided that $\gamma$ is chosen large enough. This rate is identical to what Theorem 2 predicts in the interpolation regime. However, the methods are different, and we do not need to assume interpolation regime here. Instead, we assume perfect similarity ($\delta=0$), and the availability of the gradient of $f$. So, while the rates are exactly the same, both the methods and the assumptions are different.
- 2.
General case. It can be shown that the expression $\frac{1+\gamma^{2}\delta^{2}}{\left(1+\gamma\mu\right)^{2}}$ is minimized for $\gamma=\frac{\mu}{\delta^{2}}$. With this choice of the stepsize we get
$\displaystyle{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\left
(\frac{\delta^{2}}{\delta^{2}+\mu^{2}}\right)^{k}{\left\|x_{0}-x_{\star}\right
\|^{2}}=\left(1-\frac{\mu^{2}}{\delta^{2}+\mu^{2}}\right)^{k}{\left\|x_{0}-x_{
\star}\right\|^{2}}.$
(47)
This means that
$\displaystyle k\geq\left(1+\frac{\delta^{2}}{\mu^{2}}\right)\log\left(\frac{{
\left\|x_{0}-x_{\star}\right\|^{2}}}{\varepsilon}\right)$
(48)
iterations suffice to guarantee ${\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\varepsilon$.

###### Proof of Theorem 6 .

From Lemma 6 we have that 7 holds for the iterates of SPPM-GC (Algorithm 6) with

$$ ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=\delta^{2},{ \color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}=0,\quad\text{and}\quad{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{2}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}=0,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}=0.$ $$

From Theorem 1, choosing any $\alpha>0,$ $\theta=\frac{1+\gamma^{2}\delta^{2}}{\left(1+\gamma\mu\right)^{2}},$ $\zeta=0$ (see (15) and (16)), we have

$$ $\displaystyle{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\left (\frac{1+\gamma^{2}\delta^{2}}{\left(1+\gamma\mu\right)^{2}}\right)^{k}{\left \|x_{0}-x_{\star}\right\|^{2}}.$ $$

∎

###### Lemma 6 ( SPPM-GC ) .

###### Proof.

###### Theorem 6 .

###### Proof of Theorem 6 .

### B.7 Loopless stochastic variance reduced proximal point method ( L-SVRP / SPPM-LGC )

We consider the stochastic optimization problem (1), i.e.,

$$ $\min_{x\in\mathbb{R}^{d}}\left\{f(x):={\rm E}_{\xi\sim\mathcal{D}}\left[f_{\xi }(x)\right]\right\},$ $$

and rely on 1 (differentiability of $f_{\xi}$) and 2 ($\mu$-strong convexity of $f_{\xi}$). Recall that this implies strong convexity of $f.$ Hence $f$ has a unique minimizer, which we denote $x_{\star}$.

Note that SPPM-GC needs to compute $\nabla f(x_{k})$ in iteration $k$. This can be very costly or even impossible to do in practice. To make this more clear, consider the problem

$$ $\min_{x\in\mathbb{R}^{d}}\left\{f(x)=\frac{1}{n}\sum_{i=1}^{n}f_{i}(x)\right\}$ $$

as a special case of (1).

- •
One worker. Assume we have a single machine solving this problem. Moreover, assume it takes one unit of time to this machine to compute $\nabla f_{i}$ for any $i$, and $n$ units of time to compute $\nabla f$. If the computation of $\nabla f$ is the bottleneck (i.e., if it is more expensive than the evaluation of the proximity operator of $f_{i}$), then an attempt to design a method addressing this bottleneck would be justified.
- •
Parallel workers. Assume we have $n$ workers able to work in parallel. Then $\nabla f$ can be computed in $1$ unit of time if communication among the workers is instantaneous. However, it may still be desirable to avoid having to compute the gradient:
–
The server aggregating the $n$ gradients computed by the workers may have limited capacity, and it make take more time for it to be able to compute the average of a very large number of vectors.
–
Some workers may be not available at all times.

These considerations justify the desire to reduce the reliance of SPPM-GC on the computation of $\nabla f$. The key idea is to compute the gradient only periodically, i.e., to be “lazy” about computing the gradient. In particular, we flip a biased coin, and compute a new gradient if the coin lands the right way. Otherwise, we use the previously computed gradient instead. We shall formalize this in the next section.

We are now ready to present the stochastic proximal point method with lazy gradient correction (SPPM-LGC). In the literature, the method is known under the name loopless stochastic variance reduced proximal (L-SVRP) point method.

Figure: Algorithm 7 Loopless Stochastic Variance Reduced Proximal Point Method (L-SVRP / SPPM-LGC)

Note that as intended, L-SVRP indeed reduces to SPPM-GC when $w_{0}=x_{0}$ and $p=1$.

#### Hiding the prox.

Further, recall that if for some $x\in\mathbb{R}^{d}$ and differentiable and convex $\phi$ we let $x_{+}:={\rm prox}_{\phi}\left(x\right)$, then $x_{+}=x-\nabla\phi(x_{+})$. Therefore, steps 4 and 5 of the method can be written in the equivalent form

$$ $x_{k+1}=x_{k}+\gamma h_{k}-\gamma\nabla f_{\xi_{k}}(x_{k+1})=x_{k}-\gamma\left (\nabla f_{\xi_{k}}(x_{k+1})-\nabla f_{\xi_{k}}({\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}w_{k}})+\nabla f({\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}w_{k}})\right).$ $$

#### L-SVRP vs L-SVRG.

The name L-SVRP was intentionally coined to resemble the name L-SVRG, which is a method proposed in and studied by Kovalev et al. (). This method has the form

$$ $x_{k+1}=x_{k}-\gamma\left(\nabla f_{\xi_{k}}(x_{k})-\nabla f_{\xi_{k}}({\color [rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb }{0.75390625,0.09765625,0.09765625}w_{k}})+\nabla f({\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}w_{k}})\right),$ $$

with Step 6 being identical. That is, the only difference here is that L-SVRP involves $\nabla f_{\xi_{k}}(x_{k+1})$ while L-SVRG uses $\nabla f_{\xi_{k}}(x_{k})$ in the same place. So, “P” in L-SVRP refers to the proximal nature of the term $\nabla f_{\xi_{k}}(x_{k+1})$, while “G” in L-SVRG refers to the gradient nature of the corresponding term $\nabla f_{\xi_{k}}(x_{k})$.

#### Loopless vs loopy structure.

The word “loopless” refers to the way the control vector ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}w_{k+1}}$ is updated in Step 6. The alternative to this, used in the famous SVRG method of Johnson and Zhang (), is to update ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}w_{k+1}}$ once every $m$ iterations, where $m$ is an appropriately chosen parameter. This change introduces an outer loop into the method, and makes it look a bit more cumbersome. More importantly, the loopless nature of L-SVRG is useful in three ways:

- (i)
leads to a somewhat sharper analysis,
- (ii)
makes the method easier to analyze, and
- (iii)
allows for easier to extensions / modifications.

The last two points are more important than the first one.

###### Lemma 7 ( L-SVRP ) .

Suppose Assumption 5 holds with $\delta>0.$ 7 holds for the iterates of L-SVRP (Algorithm 7) with

$$ ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{1}}=\delta^{2},{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}=0,\quad\text{and}\quad{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{2}}=p,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}=1-p,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}=0.$ $$

###### Proof of Lemma 7 .

Recall that the iterates of L-SVRP have the form

$$ $x_{k+1}={\rm prox}_{\gamma f_{\xi_{k}}}\left(x_{k}+\gamma h_{k}\right),$ $$

where $h_{k}$ is defined as
$h_{k}=\nabla f_{\xi_{k}}(w_{k})-\nabla f(w_{k}),$ and $w_{k}$ is updated in a loopless fashion.
Let $\phi_{k}=w_{k}$. Then

$$ ${\rm E}\left[\left.h_{k}\right|x_{k},\phi_{k}\right]={\rm E}\left[\left.\nabla f _{\xi_{k}}(w_{k})-\nabla f(w_{k})\right|x_{k},w_{k}\right]=0,$ $$

and hence (7) holds. If, moreover, Assumption 5 holds, then

$$ $\displaystyle{\rm E}\left[\left.{\left\|h_{k}-\nabla f_{\xi_{k}}(x_{\star}) \right\|^{2}}\right|x_{k},w_{k}\right]$ $\displaystyle\overset{\eqref{eq:GDY78d8g7d-xx}}{\leq}$ $\displaystyle\delta^{2}{\left\|w_{k}-x_{\star}\right\|^{2}},$ $$

which means that (8) holds with ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=0$, ${\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1}}=\delta^{2}$ and ${\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{1}}=0$ if we let $\sigma^{2}(z):={\left\|z-x_{\star}\right\|^{2}}$ (since then
$\sigma_{k}^{2}=\sigma^{2}(w_{k})={\left\|w_{k}-x_{\star}\right\|^{2}}$). On the other hand, from the proof of auxiliary Lemma 12 we know that

$$ $\displaystyle{\rm E}\left[\left.{\left\|w_{k+1}-x_{\star}\right\|^{2}}\right|x _{k+1},w_{k}\right]$ $\displaystyle\overset{\eqref{eq:nuD*7g7f}}{=}$ $\displaystyle p{\left\|x_{k+1}-x_{\star}\right\|^{2}}+(1-p){\left\|w_{k}-x_{ \star}\right\|^{2}},$ $$

which means that (9) holds with ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}}=p$, ${\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{2}}=1-p$, and ${\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{2}}=0$. In summary, 7 holds with

$$ ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{1}}=\delta^{2},{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}=0,\quad\text{and}\quad{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{2}}=p,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}=1-p,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}=0.$ $$

∎

###### Theorem 7 .

Let 1 (differentiability), 2 ($\mu$-strong convexity), and 5 ($\delta$-similarity) hold. Choose any $x_{0},w_{0}\in\mathbb{R}^{d}$. Then for any $p\in(0,1]$, $\gamma>0$, $\alpha>0$ and all $k\geq 0$, we have

$$ $\displaystyle{\rm E}\left[\Psi_{k}\right]\leq\max\left\{\frac{1+\alpha p}{(1+ \gamma\mu)^{2}},\frac{1+\alpha p}{(1+\gamma\mu)^{2}}\frac{\gamma^{2}\delta^{2} }{\alpha}+1-p\right\}^{k}\Psi_{0},$ (49) $$

where

$$ $\displaystyle\Psi_{k}:={\left\|x_{k}-x_{\star}\right\|^{2}}+\alpha{\left\|w_{k }-x_{\star}\right\|^{2}}.$ (50) $$

If $\delta=0$, we have the more precise result

$$ $\displaystyle{\rm E}\left[{\left\|x_{k+1}-x_{\star}\right\|^{2}}\right]\leq \frac{1}{(1+\gamma\mu)^{2}}{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}} \right].$ (51) $$

Commentary:

- 1.
Convergence vs divergence. Clearly, it is possible for the maximum in (49) to not be smaller than 1. In this case, the theorem gives a meaningless result. Whether or not the value is smaller than 1 depends on the choice of the parameters $\alpha$, $p$ and $\gamma$ in relation to the strong convexity constant $\mu$. For example, it’s clear that if $\gamma$ and $p$ are fixed and $\alpha$ is too large, the expression $\frac{1+\alpha p}{(1+\gamma\mu)^{2}}$ might exceed 1, rendering the rate vacuous.
- 2.
Optimal choice of $\alpha$. Note that $\alpha\mapsto\frac{1+\alpha p}{(1+\gamma\mu)^{2}}$ is linear and increasing, and $\alpha\mapsto\frac{1+\alpha p}{(1+\gamma\mu)^{2}}\frac{\gamma^{2}\delta^{2}}{
\alpha}+1-p$ is convex and decreasing (make sure you understand why!). Moreover, while the first function has a finite value at $\alpha=0$, the second function blows up as $\alpha$ approaches 0 0 from the right. This means that the maximum of these two functions will be minimized at the point where the graphs of the two functions intersect, i.e., at $\alpha$ satisfying
$\frac{1+\alpha p}{(1+\gamma\mu)^{2}}=\frac{1+\alpha p}{(1+\gamma\mu)^{2}}\frac
{\gamma^{2}\delta^{2}}{\alpha}+1-p.$
(52)
In the $p=1$ case (L-SVRP reduces to SPPM-GC in this regime), the equation simplifies to $1=\frac{\gamma^{2}\delta^{2}}{\alpha}$, i.e., the optimal solution is $\alpha=\gamma^{2}\delta^{2}$, and (49) reduces to
$\displaystyle{\rm E}\left[\Psi_{k}\right]\leq\left(\frac{1+\gamma^{2}\delta^{2
}}{(1+\gamma\mu)^{2}}\right)^{k}\Psi_{0}.$
(53)
This is the same result we obtained in (47) for the SPPM-GC method, up to the choice of the Lyapunov function. However, if we initialize with $w_{0}=x_{0}$, then
$\Psi_{k}=(1+\gamma\mu){\left\|x_{k}-x_{\star}\right\|^{2}},$
(54)
and plugging this into (53) gives
$\displaystyle{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\left
(\frac{1+\gamma^{2}\delta^{2}}{(1+\gamma\mu)^{2}}\right)^{k}{\left\|x_{0}-x_{
\star}\right\|^{2}},$
(55)
which is exactly (47).
The $0<p<1$ case turns out to be a more cumbersome. After a bit of algebra, one obtains that equation (52) is equivalent to the quadratic equation
$p\alpha^{2}+b\alpha+c=0,$
where $b=1-(1+\gamma\mu)^{2}(1-p)-p\gamma^{2}\delta^{2}$ and $c=\gamma^{2}\delta^{2}$. The roots of the quadratic are
$\alpha=\frac{-b\pm\sqrt{b^{2}-4pc}}{2p}.$
It seems it is awkward to work with this expression since the resulting rate will become hard to parse and interpret. So, we’ll give up on working with perfectly optimal $\alpha$. Nevertheless, we will show how to choose some (slightly suboptimal) $\alpha$ in Corollary 8 which also gives the right complexity result.

Admittedly, it’s not easy to understand how good is the rate provided by (49). The following corollary sheds light on what is achievable.

###### Corollary 8 .

If we choose $\alpha=\frac{\gamma\mu}{p}$ and $\gamma=\frac{p}{p\frac{\delta^{2}}{\mu}+(1-p)\mu}$, then for any $\varepsilon>0$ we have

$$ $k\geq\left(\frac{1}{p}+\frac{\delta^{2}}{\mu^{2}}\right)\log\left(\frac{\Psi_{ 0}}{\varepsilon}\right)\qquad\Rightarrow\qquad{\rm E}\left[\Psi_{k}\right]\leq\varepsilon.$ (56) $$

###### Proof.

Since $\alpha=\frac{\gamma\mu}{p}$, we have

$$ $A(\gamma):=\frac{1+\alpha p}{(1+\gamma\mu)^{2}}=\frac{1}{1+\gamma\mu}$ $$

and

$$ $B(\gamma):=\frac{1+\alpha p}{(1+\gamma\mu)^{2}}\frac{\gamma^{2}\delta^{2}}{ \alpha}+1-p=\frac{1}{1+\gamma\mu}\frac{p\gamma\delta^{2}}{\mu}+1-p.$ $$

Plugging this into (50) leads to

$$ $\displaystyle{\rm E}\left[\Psi_{k}\right]\leq\max\left\{A(\gamma),B(\gamma) \right\}^{k}\Psi_{0},$ (57) $$

We will now select stepsize $\gamma$ which minimizes

$$ $\gamma\mapsto\max\{A(\gamma),B(\gamma)\}.$ $$

Notice that $\gamma\mapsto A(\gamma)$ is decreasing to zero on $(0,\infty)$, with $A(\gamma)$ blowing up to $\infty$ as $\gamma$ approaches zero from the right. Further, $\gamma\mapsto B(\gamma)$ is increasing in $(0,\infty)$. This means that $\max\{A(\gamma),B(\gamma)\}$ is minimized at the point where the graphs of the two functions intersect, i.e., at $\gamma$ satisfying $A(\gamma)=B(\gamma)$. Direct calculation shows that the solution of this is

$$ $\gamma=\gamma_{\star}:=\frac{p}{p\frac{\delta^{2}}{\mu}+(1-p)\mu},$ (58) $$

and hence

$$ ${\rm E}\left[\Psi_{k}\right]\leq\max\left\{A(\gamma_{\star}),B(\gamma_{\star}) \right\}^{k}\Psi_{0}=A(\gamma_{\star})^{k}\Psi_{0}=\left(\frac{1}{1+\gamma_{ \star}\mu}\right)^{k}\Psi_{0}=\left(1-\frac{\gamma_{\star}\mu}{1+\gamma_{\star }\mu}\right)^{k}\Psi_{0}.$ $$

This implies that

$$ $k\geq\left(1+\frac{1}{\gamma_{\star}\mu}\right)\log\left(\frac{\Psi_{0}}{ \varepsilon}\right)\qquad\Rightarrow\qquad{\rm E}\left[\Psi_{k}\right]\leq\varepsilon.$ (59) $$

Plugging $\gamma_{\star}$ into this iteration complexity result gives

$$ $1+\frac{1}{\gamma_{\star}\mu}\overset{\eqref{eq:08y0hf9hdfd}}{=}1+\frac{\frac{ p\delta^{2}}{\mu}+(1-p)\mu}{p\mu}=\frac{1}{p}+\frac{\delta^{2}}{\mu^{2}}.$ $$

Plugging this back into (59) gives the final result

$$ $k\geq\left(\frac{1}{p}+\frac{\delta^{2}}{\mu^{2}}\right)\log\left(\frac{\Psi_{ 0}}{\varepsilon}\right)\qquad\Rightarrow\qquad{\rm E}\left[\Psi_{k}\right]\leq\varepsilon.$ (60) $$

∎

Commentary:

- 1.
Comparing to SPPM with Gradient Correction (i.e., $p=1$). Recall that L-SVRP reduces to SPPM-GC when $w_{0}=x_{0}$ and $p=1$. Therefore, one would expect the rates to be the same. First, notice that $w_{k}=x_{k}$ for all $k$, and as a result, the Lyapunov function (50) reduces to
$\Psi_{k}=(1+\gamma\mu){\left\|x_{k}-x_{\star}\right\|^{2}},$
(61)
and hence (57) reduces to
$\displaystyle(1+\gamma\mu){\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}
\right]\overset{\eqref{eq:Psi-simpified}}{=}{\rm E}\left[\Psi_{k}\right]\leq
\max\left\{\left(\frac{1}{1+\gamma\mu}\right)^{k},\left(\frac{1}{1+\gamma\mu}
\frac{\gamma\delta^{2}}{\mu}\right)^{k}\right\}\Psi_{0}.$
If we choose $\gamma\leq\frac{\mu}{\delta^{2}}$, then the first term in the max dominates, and we get
$\displaystyle\max\left\{\left(\frac{1}{1+\gamma\mu}\right)^{k},\left(\frac{1}{
1+\gamma\mu}\frac{\gamma\delta^{2}}{\mu}\right)^{k}\right\}\Psi_{0}$
$\displaystyle=$
$\displaystyle\left(\frac{1}{1+\gamma\mu}\right)^{k}\Psi_{0}$
$\displaystyle\overset{\eqref{eq:Psi-simpified}}{=}$
$\displaystyle\left(\frac{1}{1+\gamma\mu}\right)^{k}(1+\gamma\mu){\left\|x_{0}-
x_{\star}\right\|^{2}}.$
Combining the above observations, we get
$\displaystyle{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\left
(\frac{1}{1+\gamma\mu}\right)^{k}{\left\|x_{0}-x_{\star}\right\|^{2}}=\left(1-
\frac{1}{1+\frac{1}{\gamma\mu}}\right)^{k}{\left\|x_{0}-x_{\star}\right\|^{2}}.$
It is easy to see that this means that if we choose $k\geq\left(1+\frac{1}{\gamma\mu}\right)\log\left(\frac{{\left\|x_{0}-x_{\star}
\right\|^{2}}}{\varepsilon}\right),$ then ${\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\varepsilon$. The best rate is obtained for the largest allowed stepsize, i.e., for $\gamma=\frac{\mu}{\delta^{2}}$ (recall that this was the optimal stepsize choice for SPPM-GC established in Section B.6), in which case we conclude that
$k\geq\left(1+\frac{\delta^{2}}{\mu^{2}}\right)\log\left(\frac{{\left\|x_{0}-x_
{\star}\right\|^{2}}}{\varepsilon}\right)\qquad\Rightarrow\qquad{\rm E}\left[{
\left\|x_{k}-x_{\star}\right\|^{2}}\right]\leq\varepsilon.$
If more similarity (i.e., smaller $\delta$) or more strong convexity (i.e., larger $\mu$) is present, fewer iterations are needed to solve the problem. Note that we get the same result as in (48); so, we do not lose anything by doing the analysis in the $\delta>0$ case using the Lyapunov approach.
- 2.
Comparison to the result of Khaled and Jin (). Choosing $\alpha=\frac{\gamma\mu}{p},$ $\gamma=\frac{\mu}{2\delta^{2}},$ $p=\frac{1}{n},$ we retrieve the convergence guarantees of Khaled and Jin (). Indeed, from Corollary 8 we have that
$A(\gamma)=\frac{1}{1+\gamma\mu},\quad B(\gamma)=\frac{1}{1+\gamma\mu}\frac{
\gamma\delta^{2}p}{\mu}+1-p.$
The condition on the stepsize states that $\frac{\gamma\delta^{2}}{\mu}\leq\frac{1}{2}.$ It implies that $B(\gamma)\leq 1-\frac{p}{2}.$ Let $\rho=\min\left\{\frac{\gamma\mu}{1+\gamma\mu},\frac{p}{2}\right\}.$ Clearly, ${\rm E}\left[{\left\|x_{k}-x_{\star}\right\|}^{2}\right]\leq{\rm E}\left[\Psi_
{k}\right].$ Further, as $w_{0}=x_{0},$
${\rm E}\left[\Psi_{0}\right]={\left\|x_{0}-x_{\star}\right\|}^{2}+\frac{\gamma
\mu}{p}{\left\|w_{0}-x_{\star}\right\|}^{2}=\left(1+\frac{\gamma\mu}{p}\right)
{\left\|x_{0}-x_{\star}\right\|}^{2}.$
For any $k\geq 0,$ we obtain
${\rm E}\left[{\left\|x_{k}-x_{\star}\right\|}^{2}\right]\leq\left(1+\frac{
\gamma\mu}{p}\right)\left(1-\rho\right)^{k}{\left\|x_{0}-x_{\star}\right\|}^{2},$
therefore,
${\rm E}\left[{\left\|x_{K}-x_{\star}\right\|}^{2}\right]\leq\left(1+\frac{
\gamma\mu}{p}\right)\exp\left\{-\rho K\right\}{\left\|x_{0}-x_{\star}\right\|}
^{2}.$
If we run L-SVRP for
$K\geq\frac{1}{\rho}\log\left(\frac{{\left\|x_{0}-x_{\star}\right\|}^{2}\left(1
+\frac{\gamma\mu}{p}\right)}{\varepsilon}\right).$
Making the substitutions $\rho=\min\left\{\frac{\gamma\mu}{1+\gamma\mu},\frac{p}{2}\right\},$ $\gamma=\frac{\mu}{2\delta^{2}},$ $p=\frac{1}{n},$ we arrive at
$K\geq 2\max\left\{\frac{1}{2}+\frac{\delta^{2}}{\mu^{2}},n\right\}\log\left(
\frac{{\left\|x_{0}-x_{\star}\right\|}^{2}\left(1+\frac{\mu^{2}n}{2\delta^{2}}
\right)}{\varepsilon}\right),$
which is even slightly better than the result by Khaled and Jin ().

###### Proof of Theorem 7 ..

From Lemma 7 we know that 7 holds for the iterates of L-SVRP (Algorithm 7) with

$$ ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{1}}=\delta^{2},{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}=0,\quad\text{and}\quad{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{2}}=p,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}=1-p,{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}=0.$ $$

From Theorem 1 choosing any $\alpha>0,$ $\theta=\max\left\{\frac{1+\alpha p}{(1+\gamma\mu)^{2}},\frac{1+\alpha p}{(1+
\gamma\mu)^{2}}\frac{\gamma^{2}\delta^{2}}{\alpha}+1-p\right\},$ $\zeta=0$ (see (15) and (16)), we get

$$ $\displaystyle{\rm E}\left[\Psi_{k}\right]\leq\max\left\{\frac{1+\alpha p}{(1+ \gamma\mu)^{2}},\frac{1+\alpha p}{(1+\gamma\mu)^{2}}\frac{\gamma^{2}\delta^{2} }{\alpha}+1-p\right\}^{k}\Psi_{0},$ $$

where

$$ $\displaystyle\Psi_{k}:={\left\|x_{k}-x_{\star}\right\|^{2}}+\alpha{\left\|w_{k }-x_{\star}\right\|^{2}}.$ $$

∎

###### Lemma 7 ( L-SVRP ) .

###### Proof of Lemma 7 .

###### Theorem 7 .

###### Corollary 8 .

###### Proof.

###### Proof of Theorem 7 ..

### B.8 Point SAGA ( Point SAGA )

The main motivation is to give one more example of a SPPM method based on the idea of gradient correction which does not need to compute the full/exact gradient of $f$ in each iteration. We consider another well-known variance-reduced SPPM method called Point SAGA. However, we will revert back to the finite-sum optimiziation problem

$$ $\min_{x\in\mathbb{R}^{d}}\left\{f(x)=\frac{1}{n}\sum_{i=1}^{n}f_{i}(x)\right\}.$ $$

Figure: Algorithm 8 Point SAGA (Point SAGA)

Commentary:

- 1.
Compared to Algorithm 7, Algorithm 8 uses additional memory to store the table of control vectors $w^{i}_{k}$ or computed gradients $\nabla f_{i}(w^{i}_{k})$.
- 2.
In each iteration, the method update only one “column” in a memory table replacing the old control vector/gradient with the corresponding new one.

The structure of Point SAGA will not allow us perform the analysis under the similarity assumption (5). Instead, we will rely on the stronger 6.

We assume that there exists $\nu>0$
such that the inequality

$$ $\frac{1}{n}\sum^{n}_{j=1}{\left\|\nabla f_{j}(x^{j})-\frac{1}{n}\sum^{n}_{i=1} \nabla f_{i}(x^{i})-\nabla f_{j}(x_{\star})\right\|^{2}}\leq\nu^{2}\frac{1}{n} \sum^{n}_{j=1}{\left\|x^{j}-x_{\star}\right\|^{2}}$ $$

holds for all $x^{1},\dots,x^{n}\in\mathbb{R}^{d}$.

This inequality can be written in the form

$$ $\frac{1}{n}\sum^{n}_{j=1}{\left\|\nabla f_{j}(x^{j})-\nabla f_{j}(x_{\star}) \right\|^{2}}-{\left\|\frac{1}{n}\sum^{n}_{i=1}\nabla f_{i}(x^{i})\right\|^{2} }\leq\nu^{2}\frac{1}{n}\sum^{n}_{j=1}{\left\|x^{j}-x_{\star}\right\|^{2}}, \quad\forall x^{j}\in\mathbb{R}^{d},~{}~{}j\in[n].$ (62) $$

Thus, (62) holds, if the following condition is assumed

$$ $\frac{1}{n}\sum^{n}_{j=1}{\left\|\nabla f_{j}(x^{j})-\nabla f_{j}(x_{\star}) \right\|^{2}}\leq\nu^{2}\frac{1}{n}\sum^{n}_{j=1}{\left\|x^{j}-x_{\star}\right \|^{2}},\quad\forall x^{j}\in\mathbb{R}^{d},~{}~{}j\in[n].$ (63) $$

Moreover, (63) is equivalent to the following condition: for all $j\in[n]$, we have

$$ ${\left\|\nabla f_{j}(x)-\nabla f_{j}(x_{\star})\right\|^{2}}\leq\nu^{2}{\left \|x-x_{\star}\right\|^{2}},\quad\forall x\in\mathbb{R}^{d}.$ (64) $$

Finally, $\eqref{eq:similarity_point_saga_4}$ holds, if each $f_{j}$ is $\nu$-smooth, i.e.

$$ ${\left\|\nabla f_{j}(x)-\nabla f_{j}(y)\right\|^{2}}\leq\nu^{2}{\left\|x-y \right\|^{2}},\quad\forall x,y\in\mathbb{R}^{d}.$ (65) $$

In summary, we have the following relations between the above conditions:

$$ $\eqref{eq:similarity_point_saga_5}\quad\Rightarrow\quad\eqref{eq:similarity_ point_saga_4}\quad\equiv\quad\eqref{eq:similarity_point_saga_3}\quad \Rightarrow\quad\eqref{eq:similarity_point_saga_2}\quad\equiv\quad\eqref{eq: similarity_point_saga}.$ $$

###### Lemma 8 .

Suppose Assumption 6 holds with $\nu>0.$ Then 7 holds for the iterates of Point SAGA (Algorithm 8) with

$$ ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{1}}=\nu^{2},{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}=0,\quad\text{and}\quad{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{2}}=\frac{1}{n},{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}=\frac{n-1}{n},{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}=0.$ $$

###### Proof of Lemma 8 .

Let $\sigma_{k}=\frac{1}{n}\sum\limits_{i=1}^{n}{\left\|w_{k}^{i}-x_{\star}\right\|
^{2}},$ $\phi_{k}=\left(w_{k}^{1},\ldots,w_{k}^{n}\right).$ Recalling that

$$ $h_{k}:=\nabla f_{i_{k}}(w^{i_{k}}_{k})-\frac{1}{n}\sum\limits_{j=1}^{n}\nabla f _{j}(w^{j}_{k}),$ (66) $$

we have

$$ $\displaystyle{\rm E}\left[\left.h_{k}\right|x_{k},\phi_{k}\right]$ $\displaystyle=$ $\displaystyle{\rm E}\left[\left.\frac{1}{n}\sum\limits_{j=1}^{n}\nabla f_{j}(w ^{j}_{k})-\frac{1}{n}\sum\limits_{j=1}^{n}\nabla f_{j}(w^{j}_{k})\right|x_{k}, \phi_{k}\right]=0.$ (67) $$

Further,

$$ $\displaystyle{\rm E}\left[\left.{\left\|h_{k}-\nabla f_{i_{k}}(x_{\star}) \right\|^{2}}\right|x_{k},\phi_{k}\right]$ $\displaystyle\overset{\eqref{eq:gdhdhdjdj}}{=}$ $\displaystyle{\rm E}\left[\left.{\left\|\nabla f_{i_{k}}(w^{i_{k}}_{k})-\frac{ 1}{n}\sum\limits_{j=1}^{n}\nabla f_{j}(w^{j}_{k})-\nabla f_{i_{k}}(x_{\star}) \right\|^{2}}\right|x_{k},\phi_{k}\right]$ $\displaystyle=$ $\displaystyle\frac{1}{n}\sum\limits_{i=1}^{n}{\left\|\nabla f_{i}(w^{i}_{k})- \frac{1}{n}\sum\limits_{j=1}^{n}\nabla f_{j}(w^{j}_{k})-\nabla f_{i}(x_{\star} )\right\|^{2}}$ $\displaystyle\overset{\eqref{eq:similarity_point_saga}}{\leq}$ $\displaystyle\frac{\nu^{2}}{n}\sum\limits_{i=1}^{n}{\left\|w_{k}^{i}-x_{\star} \right\|^{2}}.$ $$

Therefore, we have that ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=0,$ ${\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1}}=\nu^{2},$ ${\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{1}}=0.$

$$ $\displaystyle{\rm E}\left[\left.\sigma_{k+1}^{2}\right|x_{k+1},\phi_{k}\right]$ $\displaystyle=$ $\displaystyle{\rm E}\left[\left.\frac{1}{n}\sum_{i=1}^{n}{\left\|w_{k+1}^{i}-x _{\star}\right\|}^{2}\right|x_{k+1},\phi_{k}\right]$ $\displaystyle=$ $\displaystyle\frac{1}{n}\sum_{i_{k}=1}^{n}\left[\frac{1}{n}{\left\|x_{k+1}-x_{ \star}\right\|}^{2}+\frac{1}{n}\sum_{j\neq i_{k}}{\left\|w_{k}^{j}-x_{\star} \right\|}^{2}\right]$ $\displaystyle=$ $\displaystyle\frac{1}{n}{\left\|x_{k+1}-x_{\star}\right\|}^{2}+\frac{n-1}{n} \sigma_{k}^{2}.$ $$

Therefore, we have that ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{2}}=\frac{1}{n},$ ${\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{2}}=\frac{n-1}{n},$ ${\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{
pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{2}}=0.$
∎

The convergence of Point SAGA is captured by the following theorem.

###### Theorem 9 .

Let 1, 2 and 6 hold. Chose any $x_{0},w^{1}_{0},\ldots,w^{n}_{0}\in\mathbb{R}^{d}$. Then for any $\gamma>0$, and all $k\geq 0$, we have

$$ ${\rm E}\left[\Psi_{k}\right]\leq\max\left\{\left(\frac{1}{1+\gamma\mu}\right)^ {k},\left(\frac{1}{1+\gamma\mu}\frac{\gamma\nu^{2}}{\mu n}+1-\frac{1}{n}\right )^{k}\right\}\Psi_{0},$ (68) $$

where

$$ $\Psi_{k}:={\left\|x_{k}-x_{\star}\right\|^{2}}+\gamma\mu\sum_{i=1}^{n}{\left\| w^{i}_{k}-x_{\star}\right\|^{2}}.$ (69) $$

Clearly, it is possible for the maximum in (68) to not be smaller than $1.$ In this case, the theorem produces the meaningless result. Whether or not the value is smaller than $1$ dependes on the choice of $\gamma$ with respect to the strong convexity constant $\mu,$ the number of individual functions $n,$ the similarity constant $\nu.$

###### Proof of Theorem 9 .

From Lemma 8 we have that 7 holds for the iterates of Point SAGA (Algorithm 8) with

$$ ${\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}=0,{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{1}}=\nu^{2},{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}=0,\quad\text{and}\quad{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{2}}=\frac{1}{n},{\color[rgb]{ 0.09765625,0.09765625,0.75390625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.09765625,0.75390625}B_{2}}=\frac{n-1}{n},{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{2}}=0.$ $$

From Theorem 1, choosing $\alpha=\gamma\mu n,$ $\theta=\max\left\{\frac{1}{1+\gamma\mu},\frac{1}{1+\gamma\mu}\frac{\gamma\nu^{
2}}{\mu n}+1-\frac{1}{n}\right\},$ $\zeta=0$ (see (15) and (16)), we obtain

$$ ${\rm E}\left[\Psi_{k}\right]\leq\max\left\{\left(\frac{1}{1+\gamma\mu}\right)^ {k},\left(\frac{1}{1+\gamma\mu}\frac{\gamma\nu^{2}}{\mu n}+1-\frac{1}{n}\right )^{k}\right\}\Psi_{0},$ $$

where

$$ $\Psi_{k}:={\left\|x_{k}-x_{\star}\right\|^{2}}+\gamma\mu\sum_{i=1}^{n}{\left\| w^{i}_{k}-x_{\star}\right\|^{2}}.$ $$

∎

###### Corollary 10 .

If we choose $\gamma=\frac{1}{\frac{\gamma^{2}}{\mu}+(n-1)\mu},$ then, for any $\varepsilon>0,$ we have

$$ $k\geq\left(n+\frac{\nu^{2}}{\mu^{2}}\right)\log\left(\frac{\Psi_{0}}{ \varepsilon}\right)\qquad\Rightarrow\qquad{\rm E}\left[\Psi_{k}\right]\leq\varepsilon.$ $$

###### Proof of Corollary 10 .

Notice that $A(\gamma):=\frac{1}{1+\gamma\mu}$ is decreasint for $\gamma>0,$ and $B(\gamma):=\frac{\gamma\nu^{2}}{1+\gamma\mu}+1-\frac{1}{n}$ is increasing for $\gamma>0.$ This means that $\max\left\{\right\}$ is minimized at $\gamma:=\gamma_{\star}$ where $A(\gamma)=B(\gamma).$ Direct calculation shows that the solution of this is

$$ $\gamma=\gamma_{\star}:=\frac{1}{\frac{\nu^{2}}{\mu}+\left(n-1\right)\mu},$ $$

and hence

$$ ${\rm E}\left[\Psi_{k}\right]\leq\max\left\{A(\gamma_{\star}),B(\gamma_{\star}) \right\}^{k}\Psi_{0}=A(\gamma_{\star})^{k}\Psi_{0}=\left(\frac{1}{1+\gamma_{ \star}\mu}\right)^{k}\Psi_{0}=\left(1-\frac{\gamma_{\star}\mu}{1+\gamma_{\star }\mu}\right)^{k}\Psi_{0}.$ $$

This implies that

$$ $k\geq\left(1+\frac{1}{\gamma_{\star}\mu}\right)\log\left(\frac{\Psi_{0}}{ \varepsilon}\right)\qquad\Rightarrow\qquad{\rm E}\left[\Psi_{k}\right]\leq\varepsilon.$ $$

Plugging $\gamma_{\star}$ into this iteration complexity result gives

$$ $1+\frac{1}{\gamma_{\star}\mu}=1+\frac{\frac{\nu^{2}}{\mu n}+\left(1-\frac{1}{n }\right)\mu}{\frac{\mu}{n}}=n+\frac{\nu^{2}}{\mu^{2}}.$ $$

We obtain the final result

$$ $k\geq\left(n+\frac{\delta^{2}}{\mu^{2}}\right)\log\left(\frac{\Psi_{0}}{ \varepsilon}\right)\qquad\Rightarrow\qquad{\rm E}\left[\Psi_{k}\right]\leq\varepsilon.$ $$

∎

###### Lemma 8 .

###### Proof of Lemma 8 .

###### Theorem 9 .

###### Proof of Theorem 9 .

###### Corollary 10 .

###### Proof of Corollary 10 .

## Appendix C Auxiliary lemma for SPPM-LC

Our main result relies on a single lemma only.

###### Lemma 9 .

Let 1, 2 and Assumption 7 hold. Then

$$ ${\rm E}\left[{\left\|x_{k+1}-x_{\star}\right\|^{2}}\right]\leq\frac{(1+\gamma^ {2}{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}})}{(1+\gamma\mu)^{ 2}}{\rm E}\left[{\left\|x_{k}-x_{\star}\right\|^{2}}\right]+\frac{\gamma^{2}{ \color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1}}}{(1+\gamma\mu)^{2 }}{\rm E}\left[\sigma^{2}_{k}\right]+\frac{\gamma^{2}{\color[rgb]{ 0.09765625,0.75390625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.09765625,0.75390625,0.09765625}C_{1}}}{(1+\gamma\mu)^{2}}.$ (70) $$

###### Proof.

By combining 1 and 2, we get

$$ $\displaystyle{\left\|x_{k+1}-x_{\star}\right\|^{2}}$ $\displaystyle\overset{\lx@cref{creftypecap~refnum}{fact:ascent_fixed_point}}{=}$ $\displaystyle{\left\|{\rm prox}_{\gamma f_{\xi_{k}}}(x_{k}+\gamma h_{k})-{\rm prox }_{\gamma f_{\xi_{k}}}(x_{\star}+\gamma\nabla f_{\xi_{k}}(x_{\star}))\right\|^ {2}}$ $\displaystyle\overset{\lx@cref{creftypecap~refnum}{fact:prox-contraction}}{\leq}$ $\displaystyle\frac{1}{(1+\gamma\mu)^{2}}{\left\|x_{k}+\gamma h_{k}-(x_{\star}+ \gamma\nabla f_{\xi_{k}}(x_{\star}))\right\|^{2}}$ $\displaystyle=$ $\displaystyle\frac{1}{(1+\gamma\mu)^{2}}{\left\|x_{k}-x_{\star}+\gamma\left(h_ {k}-\nabla f_{\xi_{k}}(x_{\star})\right)\right\|^{2}}.$ $$

Thus, we have that

$$ ${\left\|x_{k+1}-x_{\star}\right\|^{2}}\leq\frac{1}{(1+\gamma\mu)^{2}}\left({ \left\|x_{k}-x_{\star}\right\|^{2}}+2\gamma\left\langle h_{k}-\nabla f_{\xi_{k }}(x_{\star}),x_{k}-x_{\star}\right\rangle+\gamma^{2}{\left\|h_{k}-\nabla f_{ \xi_{k}}(x_{\star})\right\|^{2}}\right).$ $$

We can use it since it holds irrespective of the choice of $h_{k}$. Taking conditional expectation on both sides, we get

$$ $\displaystyle{\rm E}\left[\left.{\left\|x_{k+1}-x_{\star}\right\|^{2}}\right|x _{k},\phi_{k}\right]$ $\displaystyle\leq$ $\displaystyle\frac{1}{(1+\gamma\mu)^{2}}{\rm E}\left[\left.{\left\|x_{k}-x_{ \star}\right\|^{2}}\right|x_{k},\phi_{k}\right]$ (71) $\displaystyle\qquad+\frac{2\gamma}{(1+\gamma\mu)^{2}}\left\langle{\rm E}\left[ \left.h_{k}-\nabla f_{\xi_{k}}(x_{\star})\right|x_{k},\phi_{k}\right],x_{k}-x_ {\star}\right\rangle$ $\displaystyle\qquad+\frac{\gamma^{2}}{(1+\gamma\mu)^{2}}{\rm E}\left[\left.{ \left\|h_{k}-\nabla f_{\xi_{k}}(x_{\star})\right\|^{2}}\right|x_{k},\phi_{k} \right].$ $$

Note that

$$ $\displaystyle{\rm E}\left[\left.{\left\|x_{k}-x_{\star}\right\|^{2}}\right|x_{ k},\phi_{k}\right]$ $\displaystyle=$ $\displaystyle{\left\|x_{k}-x_{\star}\right\|^{2}}.$ (72) $$

Further,

$$ $\displaystyle{\rm E}\left[\left.h_{k}-\nabla f_{\xi_{k}}(x_{\star})\right|x_{k },\phi_{k}\right]$ $\displaystyle=$ $\displaystyle{\rm E}\left[\left.h_{k}\right|x_{k},\phi_{k}\right]-{\rm E}\left [\left.\nabla f_{\xi_{k}}(x_{\star})\right|x_{k},\phi_{k}\right]$ (73) $\displaystyle=$ $\displaystyle{\rm E}\left[\left.h_{k}\right|x_{k},\phi_{k}\right]-\nabla f(x_{ \star})$ $\displaystyle=$ $\displaystyle{\rm E}\left[\left.h_{k}\right|x_{k},\phi_{k}\right]$ $\displaystyle=$ $\displaystyle 0,$ $$

where the last equality follows from Assumption 7, relation (7). Relation (8) of Assumption 7 says that

$$ $\displaystyle{\rm E}\left[\left.{\left\|h_{k}-\nabla f_{\xi_{k}}(x_{\star}) \right\|^{2}}\right|x_{k},\phi_{k}\right]$ $\displaystyle\leq$ $\displaystyle{\color[rgb]{0.75390625,0.09765625,0.09765625}\definecolor[named] {pgfstrokecolor}{rgb}{0.75390625,0.09765625,0.09765625}A_{1}}{\left\|x_{k}-x_{ \star}\right\|^{2}}+{\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor [named]{pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1}}\sigma^{2} _{k}+{\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{1}}.$ (74) $$

Plugging (74), (73) and (72) into (71), we obtain

$$ ${\rm E}\left[\left.{\left\|x_{k+1}-x_{\star}\right\|^{2}}\right|x_{k},\phi_{k} \right]\leq\frac{1}{(1+\gamma\mu)^{2}}\left(\left(1+\gamma^{2}{\color[rgb]{ 0.75390625,0.09765625,0.09765625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.75390625,0.09765625,0.09765625}A_{1}}\right){\left\|x_{k}-x_{\star}\right\|^ {2}}+\gamma^{2}{\color[rgb]{0.09765625,0.09765625,0.75390625}\definecolor[ named]{pgfstrokecolor}{rgb}{0.09765625,0.09765625,0.75390625}B_{1}}\sigma^{2}_ {k}+\gamma^{2}{\color[rgb]{0.09765625,0.75390625,0.09765625}\definecolor[named ]{pgfstrokecolor}{rgb}{0.09765625,0.75390625,0.09765625}C_{1}}\right).$ (75) $$

It only remains to take expectation on both sides and apply the tower property. ∎

###### Lemma 9 .

###### Proof.

## Appendix D Auxiliary lemmas

### D.1 SPPM-AS

###### Lemma 10 .

Let $\phi_{1},\dots,\phi_{m}:\mathbb{R}^{d}\to\mathbb{R}$ be differentiable functions, with $\phi_{i}$ being $\mu_{i}$-strongly convex for all $i\in[n]$. Further, let $w_{1},\dots,w_{m}$ be positive scalars. Then the function $\phi:=\sum_{i=1}^{m}w_{i}\phi_{i}$ is $\mu$-strongly convex with $\mu=\sum_{i=1}^{m}w_{i}\mu_{i}$.

###### Proof.

By assumption,

$$ $\phi_{i}(y)+\left\langle\nabla\phi_{i}(y),x-y\right\rangle+\frac{\mu_{i}}{2}{ \left\|x-y\right\|^{2}}\leq\phi_{i}(x),\qquad\forall x,y\in\mathbb{R}^{d}.$ (76) $$

This means that

$$ $\sum_{i=1}^{m}w_{i}\left(\phi_{i}(y)+\left\langle\nabla\phi_{i}(y),x-y\right \rangle+\frac{\mu_{i}}{2}{\left\|x-y\right\|^{2}}\right)\leq\sum_{i=1}^{m}w_{i }\phi_{i}(x),\qquad\forall x,y\in\mathbb{R}^{d},$ $$

which is equivalent to

$$ $\phi(y)+\left\langle\nabla\phi(y),x-y\right\rangle+\frac{\sum_{i=1}^{m}w_{i} \mu_{i}}{2}{\left\|x-y\right\|^{2}}\leq\phi(x),\qquad\forall x,y\in\mathbb{R}^ {d},$ (77) $$

So, $\phi$ is $\mu$-strongly convex.
∎

###### Lemma 10 .

###### Proof.

### D.2 L-SVRP

###### Lemma 11 .

Recalling that

$$ $h_{k}:=\nabla f_{\xi_{k}}(w_{k})-\nabla f(w_{k}),$ (78) $$

we can write

$$ $\displaystyle{\rm E}\left[\left.\left\langle h_{k}-\nabla f_{\xi_{k}}(x_{\star }),x_{k}-x_{\star}\right\rangle\right|x_{k},w_{k}\right]$ $\displaystyle=$ $\displaystyle\left\langle{\rm E}\left[\left.h_{k}-\nabla f_{\xi_{k}}(x_{\star} )\right|x_{k},w_{k}\right],x_{k}-x_{\star}\right\rangle$ (79) $\displaystyle\overset{\eqref{eq:h_k-9u980fd-LSVRP}}{=}$ $\displaystyle\left\langle{\rm E}\left[\left.\nabla f_{\xi_{k}}(w_{k})-\nabla f (w_{k})-\nabla f_{\xi_{k}}(x_{\star})\right|x_{k},w_{k}\right],x_{k}-x_{\star}\right\rangle$ $\displaystyle=$ $\displaystyle\left\langle\underbrace{\nabla f(w_{k})-\nabla f(w_{k})-\nabla f( x_{\star})}_{=0},x_{k}-x_{\star}\right\rangle$ $\displaystyle=$ $\displaystyle 0,$ $$

and

$$ $\displaystyle{\rm E}\left[\left.{\left\|h_{k}-\nabla f_{\xi_{k}}(x_{\star}) \right\|^{2}}\right|x_{k},w_{k}\right]$ $\displaystyle\overset{\eqref{eq:h_k-9u980fd-LSVRP}}{=}$ $\displaystyle{\rm E}\left[\left.{\left\|\nabla f_{\xi_{k}}(w_{k})-\nabla f(w_{ k})-\nabla f_{\xi_{k}}(x_{\star})\right\|^{2}}\right|x_{k},w_{k}\right]$ (80) $\displaystyle\overset{\eqref{eq:similarity-1}}{\leq}$ $\displaystyle\delta^{2}{\left\|w_{k}-x_{\star}\right\|^{2}}.$ $$

###### Lemma 12 .

Observe that by the way $w_{k+1}$ is defined, we have

$$ $\displaystyle{\rm E}\left[\left.{\left\|w_{k+1}-x_{\star}\right\|^{2}}\right|x _{k+1},w_{k}\right]=p{\left\|x_{k+1}-x_{\star}\right\|^{2}}+(1-p){\left\|w_{k} -x_{\star}\right\|^{2}}.$ (81) $$

Taking expectation again, and applying the tower property of expectation, we get

$$ $\displaystyle{\rm E}\left[{\left\|w_{k+1}-x_{\star}\right\|^{2}}\right]$ $\displaystyle=$ $\displaystyle{\rm E}\left[{\rm E}\left[\left.{\left\|w_{k+1}-x_{\star}\right\| ^{2}}\right|x_{k+1},w_{k}\right]\right]$ (82) $\displaystyle\overset{\eqref{eq:nuD*7g7f}}{=}$ $\displaystyle p{\rm E}\left[{\left\|x_{k+1}-x_{\star}\right\|^{2}}\right]+(1-p ){\rm E}\left[{\left\|w_{k}-x_{\star}\right\|^{2}}\right],$ $$

which is what we set out to prove.

###### Lemma 11 .

###### Lemma 12 .

## Appendix E Auxiliary facts

###### Fact 1 (Every point is a fixed point) .

Let $\phi:\mathbb{R}^{d}\to\mathbb{R}$ be a differentiable convex function. Then

$$ ${\rm prox}_{\gamma\phi}\left(x+\gamma\nabla\phi(x)\right)=x,\qquad\forall \gamma>0,\quad\forall x\in\mathbb{R}^{d}.$ $$

In particular, if $x_{\star}$ is a minimizer of $\phi$, then ${\rm prox}_{\gamma\phi}\left(x_{\star}\right)=x_{\star}$.

###### Proof.

Pick any $x\in\mathbb{R}^{d}$ and $\gamma>0$.
Evaluating the proximity operator at

$$ $y:=x+\gamma\nabla\phi(x)$ (83) $$

gives

$$ $\displaystyle{\rm prox}_{\gamma\phi}\left(y\right)=\arg\min_{x^{\prime}\in \mathbb{R}^{d}}\left(\phi(x^{\prime})+\frac{1}{2\gamma}{\left\|x^{\prime}-y \right\|^{2}}\right).$ $$

This is a strongly convex minimization problem, and hence the (necessarily unique) minimizer

$$ $x^{\prime}:={\rm prox}_{\gamma\phi}\left(y\right)$ (84) $$

of this problem satisfies the first-order optimality condition

$$ $\nabla\phi(x^{\prime})+\frac{1}{\gamma}\left(x^{\prime}-y\right)=0.$ $$

Note that $x^{\prime}=x$ satisfies this equation, and hence

$$ $x\overset{\eqref{eq:98g9d_98gfd}}{=}{\rm prox}_{\gamma\phi}\left(y\right) \overset{\eqref{eq:b98g(*-98hf}}{=}{\rm prox}_{\gamma\phi}\left(x+\gamma\nabla \phi(x)\right).$ $$

∎

The next statement is (Khaled and Jin, , Fact 4)

###### Fact 2 (Contractivity of the prox) .

If $\phi$ is differentiable and $\mu$-strongly convex, then for all $\gamma>0$ and for any $x,y\in\mathbb{R}^{d}$ we have

$$ $\displaystyle{\left\|{\rm prox}_{\gamma\phi}\left(x\right)-{\rm prox}_{\gamma \phi}\left(y\right)\right\|^{2}}\leq\frac{1}{(1+\gamma\mu)^{2}}{\left\|x-y \right\|^{2}}.$ $$

###### Proof.

This lemma can be seen as a tighter version of (Mishchenko et al., , Lemma 5) though our proof technique is different. Note that $p(x)={\rm prox}_{\gamma h}(x)$ satisfies $\gamma\nabla h(p(x))+\left[p(x)-x\right]=0$, or equivalently $p(x)=x-\gamma\nabla h(p(x))$. Using this we have

$$ $\displaystyle{\left\|p(x)-p(y)\right\|^{2}}$ $\displaystyle={\left\|\left[x-\gamma\nabla h(p(x))\right]-\left[y-\gamma\nabla h (p(y))\right]\right\|^{2}}$ $\displaystyle={\left\|\left[x-y\right]-\gamma\left[\nabla h(p(x))-\nabla h(p(y ))\right]\right\|^{2}}$ $\displaystyle={\left\|x-y\right\|^{2}}+\gamma^{2}{\left\|\nabla h(p(x))-\nabla h (p(y))\right\|^{2}}-2\gamma\left\langle x-y,\nabla h(p(x))-\nabla h(p(y)) \right\rangle.$ (85) $$

Now note that

$$ $\displaystyle\left\langle x-y,\nabla h(p(x))-\nabla h(p(y))\right\rangle$ $\displaystyle=\left\langle p(x)+\gamma\nabla h(p(x))-\left[p(y)+\gamma\nabla h (p(y))\right],\nabla h(p(x))-\nabla h(p(y))\right\rangle$ $\displaystyle=\left\langle p(x)-p(y),\nabla h(p(x))-\nabla h(p(y))\right \rangle+\gamma{\left\|\nabla h(p(x))-\nabla h(p(y))\right\|^{2}}.$ (86) $$

Combining Equations 86 and 85 we get

$$ $\displaystyle\begin{split}{\left\|p(x)-p(y)\right\|^{2}}&={\left\|x-y\right\|^ {2}}+\gamma^{2}{\left\|\nabla h(p(x))-\nabla h(p(y))\right\|^{2}}-2\gamma\left \langle p(x)-p(y),\nabla h(p(x))-\nabla h(p(y))\right\rangle\\ &\qquad-2\gamma^{2}{\left\|\nabla h(p(x))-\nabla h(p(y))\right\|^{2}}\end{split}$ $\displaystyle={\left\|x-y\right\|^{2}}-\gamma^{2}{\left\|\nabla h(p(x))-\nabla h (p(y))\right\|^{2}}-2\gamma\left\langle p(x)-p(y),\nabla h(p(x))-\nabla h(p(y) )\right\rangle.$ (87) $$

Let $D_{h}(u,v)=h(u)-h(v)-\left\langle\nabla h(v),u-v\right\rangle$ be the Bregman divergence associated with $h$ at $u,v$. It is easy to show that

$$ $\displaystyle\left\langle u-v,\nabla h(u)-\nabla h(v)\right\rangle=D_{h}(u,v)+ D_{h}(v,u).$ $$

This is a special case of the three-point identity (Chen and Teboulle, , Lemma 3.1). Using this with $u=p(x)$ and $v=p(y)$ and plugging back into (87) we get

$$ $\displaystyle{\left\|p(x)-p(y)\right\|^{2}}$ $\displaystyle={\left\|x-y\right\|^{2}}-\gamma^{2}{\left\|\nabla h(p(x))-\nabla h (p(y))\right\|^{2}}-2\gamma\left[D_{h}(p(x),p(y))+D_{h}(p(y),p(x))\right].$ $$

Note that because $h$ is strongly convex, we have that $D_{h}(p(y),p(x))\geq\frac{\mu}{2}{\left\|p(y)-p(x)\right\|^{2}}$ and $D_{h}(p(x),p(y))\geq\frac{\mu}{2}{\left\|p(y)-p(x)\right\|^{2}}$, hence

$$ $\displaystyle{\left\|p(x)-p(y)\right\|^{2}}$ $\displaystyle\leq{\left\|x-y\right\|^{2}}-\gamma^{2}{\left\|\nabla h(p(x))- \nabla h(p(y))\right\|^{2}}-2\gamma\mu{\left\|p(x)-p(y)\right\|^{2}}.$ (88) $$

Strong convexity implies that for any two points $u,v$

$$ $\displaystyle{\left\|\nabla h(u)-\nabla h(v)\right\|^{2}}\geq\mu^{2}{\left\|u- v\right\|^{2}},$ $$

see (Nesterov, , Theorem 2.1.10) for a proof. Using this in Equation 88 with $u=p(x)$ and $v=p(y)$ yields

$$ $\displaystyle{\left\|p(x)-p(y)\right\|^{2}}$ $\displaystyle\leq{\left\|x-y\right\|^{2}}-\gamma^{2}\mu^{2}{\left\|p(x)-p(y) \right\|^{2}}-2\gamma\mu{\left\|p(x)-p(y)\right\|^{2}}.$ $$

Rearranging gives

$$ $\displaystyle\left[1+\gamma^{2}\mu^{2}+2\gamma\mu\right]{\left\|p(x)-p(y) \right\|^{2}}\leq{\left\|x-y\right\|^{2}}.$ $$

It remains to notice that $(1+\gamma\mu)^{2}=1+\gamma^{2}\mu^{2}+2\gamma\mu$.
∎

###### Fact 3 (Recurrence) .

Assume that a sequence $\{s_{k}\}_{k\geq 0}$ of positive real numbers for all $k\geq 0$ satisfies

$$ $\displaystyle s_{k+1}\leq as_{k}+b,$ $$

where $0<a<1$ and $b\geq 0$. Then the sequence for all $k\geq 0$ satisfies

$$ $s_{k}\leq a^{k}s_{0}+b\min\left\{k,\frac{1}{1-a}\right\}.$ (89) $$

###### Proof.

Unrolling the recurrence, we get

$$ $\displaystyle s_{k}$ $\displaystyle\leq as_{k-1}+b\leq a\left(as_{k-2}+b\right)+b\leq\cdots\leq a^{k }s_{0}+b\sum_{i=0}^{k-1}a^{i}.$ (90) $$

We can now bound the sum $\sum_{i=0}^{k-1}a^{i}$ in two different ways. First, since $a<1$, we get the estimate

$$ $\displaystyle\sum_{i=0}^{k-1}a^{i}\leq\sum_{i=0}^{k-1}1=k.$ (91) $$

Second, we sum a geometric series

$$ $\displaystyle\sum_{i=0}^{k-1}a^{i}\leq\sum_{i=0}^{\infty}a^{i}=\frac{1}{1-a}.$ (92) $$

Note that either of these bounds can be better. So, we apply the best of these bounds. Substituting Equations 91 and 92 into (90) gives (89).
∎

###### Fact 1 (Every point is a fixed point) .

###### Proof.

###### Fact 2 (Contractivity of the prox) .

###### Proof.

###### Fact 3 (Recurrence) .

###### Proof.

## Peer Reviews

confidence: 5
contribution: 2
initial_rating: 3
presentation: 4
questions: 1) Can you provide  class  of problems satisfying Assumption 5 but not the one in  Khaled and Jin ?
2) The experimental setting is not very interesting since it deals with a function with a LIpschitz continuous gradient (see eq 17).
Can you provide experiments in a setting covered by your assumptions but not by previous existing papers?
soundness: 3
strengths: The main technical contribution of the paper is a generalization of an existing framework from functions with a Lipschitz continuous gradient to the case of differentiable functions. A non negligible contribution in my opinion is the very clear presentation with several  interesting remarks. 

The assumptions and the statements are clear and the proofs technically sound.
summary: The paper deals with variance reduction techniques for stochastic proximal point algorithms, which recently emerged as a robust alternative to stochastic gradient algorithms, primarily due to their reduced sensitivity to step size. The authors develop a versatile stochastic proximal point algorithm, which can encompass both the standard stochastic proximal point algorithm, and various variance 
reduced versions.  For this algorithm, and several specific instances of it,  assuming differentiability and strong convexity of the objective function, the authors derive several convergence rates for the expected objective function values, showing improvements over the standard stochastic proximal point algorithm and available theoretical results. The paper contains a section devoted to numerical experiments elucidating some theoretical aspects.
weaknesses: The main novelty of the paper is the analysis of a unifying algorithm that allows to deal with variance reduced stochastic proximal point methods for an objective function which is  only differentiable and strongly convex. 

Though the main proofs are partially different from the ones used in the related literature, the idea is not new and it is a generalization of the approach proposed in  the papers: 

1) E. Gorbunov, F. Hanzely, and P. Richtarik. A unified theory of sgd: Variance reduction, sampling, quantization and coordinate descent.

2) C. Traore, V. Apidopoulos, S. Salzo, and S. Villa. Variance reduction techniques for stochastic
proximal point algorithms

So the novelty of the contribution is limited. 

In general, there is a tendency in the paper to overstate the contributions, which I did not appreciate. I enclose some examples below: 

1) The title claims "A unified theory of stochastic proximal poin methods without smoothness". This is misleading since the analysis is performed only for smooth (differentiable) functions, while nonsmooth ones are not covered.  With respect to the existing literature, the authors drop the assumption of Lipschitz continuity of the gradient, but in turn assume that the function is strongly convex.

2) the similarity assumption (Assumption 5) is indeed more general than the one in Khaled and Jin, but since it involves unknown quantities it is very unclear how more general it is in practice. Concrete examples of problems satisfying Assumption 5 but not the one in  Khaled and Jin must be provided in order to establish whether this generalization has some kind of impact, or not. 

3) the authors often state that a unifying framework for the study of variance reduction techniques is missing, but this is false, since both the papers mentioned above provides such framework (of course, under different assumptions on the objective function). We find out that the proposed framework s is similar to the ones of the papers 1) and 2) above only on line 174, after a long discussion. 
 
4) In the table 1:


The authors show convergence of 8 algorithms (see Table 1). However: 

Algorithm 1 is the "universal" one including all the others; 

SPPM* and SPPM GC are not interesting from the practical point of view since it requires the gradients at the solution

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: My main consideration is the importance of the Assumptions 5 and 6. Although these assumptions is more general than smooth condition, it is unclear whether they are useful in real-applications. It looks that Assumptions 5 and 6 are proposed only for convergence analysis.
1. Can you provide some applications which satisfy Assumption 5 or 6 and but do not hold smoothness?
2. The experiments on linear regression look not very relevant to this paper, since it can be solved by the algorithms based on smoothness assumption.
soundness: 3
strengths: The presentation is good. The authors provide solid theoretical analysis to support the proposed framework.
summary: This paper provided a framework for analyzing the stochastic proximal point methods for solving strongly convex optimization problem without smoothness. It leads to several new algorithm by introducing the techniques of sampling strategies, variance-reduction and correction.
weaknesses: The motivation of this paper is unclear. The assumption without smoothness used in this paper looks not popular.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: **Q1.** Is the combination of Assumption 6 and strong convexity a weaker condition than the $\nu$-smooth condition? If so, could you provide an example?

**Q2.** Could the proposed SPPM-LC framework be extended or adapted for non-convex optimization, such as weakly convex optimization?

**Q3.** How do the proposed algorithms compare with the proximal point method in [1], which is applicable to weakly convex functions using a diminishing step size?

**References:**
[1] Davis, D., & Drusvyatskiy, D. Stochastic model-based minimization of weakly convex functions. *SIOPT*, 2019.
soundness: 2
strengths: **S1:** This paper is generally well-written and easy to follow. It is well-grounded in theoretical analysis, establishing convergence guarantees without relying on smoothness assumptions.

**S2:** It unifies multiple SPPM variants under a single theoretical framework, making it easier to understand the relationships and convergence behavior across methods.

**S3:** The development of new SPPM variants, such as SPPM with Nonuniform Sampling and SPPM with Arbitrary Sampling, enriches the field by expanding the applicability of SPPM.
summary: The paper presents a unified theoretical framework for stochastic proximal point methods (SPPM) applied to problems that do not require smoothness assumptions but do rely on a strong convexity assumption.

The paper establishes a single convergence theorem applicable to multiple SPPM variants, providing the robust convergence behavior of these methods.

They propose a general algorithm called SPPM-LC (SPPM with Learned Correction) and integrates techniques such as variance reduction and arbitrary sampling through a flexible correction vector. Several new SPPM variants have been considered.

Experiments have been conducted to show some inherent properties of these approaches.
weaknesses: **W1.** This paper appears more like a review article, which may not align well with the scope of ICLR. The authors propose seven algorithms for solving differentiable convex problems. However, these algorithms may lack sufficient novelty or clear performance advantages. Readers might still be uncertain about which algorithm is the best choice.

**W2.** The framework assumes strong convexity, limiting the applicability of the results to problems that are not strongly convex or to non-convex settings.

**W3.** Algorithm 1 uses a constant, sufficiently small step size, which is rarely practical in real applications. This can introduce additional estimation errors, as highlighted in Theorem 1. In practice, a more effective approach would involve a diminishing step size [1]. The authors should at least discuss the comparison between constant and diminishing step sizes.

**W4.** The authors claim to use a weaker assumption (without smoothness); however, in their experiment, they still consider a convex least squares problem that satisfies the smoothness assumption. I think the authors should provide an example that satisfies both strong convexity and Assumption 6 but violates the $\nu$-smoothness assumption.

**W5.** The empirical experiments provided by the authors are limited to strongly convex least squares problems, which may restrict their relevance to practical applications and broader research interests.


**References:**
[1] Davis, D., & Drusvyatskiy, D. Stochastic model-based minimization of weakly convex functions. *SIOPT*, 2019.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: Can the main result in Theorem 1 be strengthened to interpolate between the sublinear rate of the standard SPPMs and the linear rate of SPPMs with gradient correction or variance reduction?
soundness: 2
strengths: The proposed universal SPPM algorithmic framework sounds interesting, and novel as far as the reviewer knows about. The convergence analysis is generally well presented, with sufficient details provided in  a clear and neat way in the appendix sections.
summary: This paper presents a unified algorithmic framework of stochastic proximal point methods that encapsulates several existing SPPM-style algorithms, and some new extensions with variance reduction and arbitrary sampling as well. The proposed universal algorithm comes with a convergence analysis under relatively mild conditions, without imposing individual smoothness assumption on loss functions. A preliminary numerical study is carried out to demonstrate the performance of newly introduced algorithms on some toy linear regression tasks.
weaknesses: The major concern is about the strength of the main result in Theorem 1. While exhibiting contraction rates, the bound there seems non-vanishing with respect to time instance unless $C_1=C_2=0$.  As shown in Sectons B.6-B.8 that such a requirement can be fulfilled by variance reduction variants of SPPM, but not by the standard SPPM (with $h_k \equiv 0$) as shown in Section B.2 (which is not surprising due to the fundamental limit of pure stochastic optimization methods). It will be much more interesting if the convergence bound could be made adaptive to standard SPPM with sublinear rate of convergence. 

The paper also falls short on the experiment side, where only toy linear regression problems are considered with very limited data scale. Since the SPPM-type algorithms have long been acknowledged to be superior to SGD in non-smooth settings, it is more desirable to test the performances of the proposed algorithms in non-smooth learning tasks such as SVMs with hinge loss.

Meta-review: A framework of convergence analysis is presented for a range of variations of SPPM, such as variance reduction and Point-SAGA. While the analysis is sound, reviewers assess that the novelty is limited and the contribution is a little overstated.

Experiments are limited to the simple least squares problem.

Decision: Reject