# Sharper Analysis for Minibatch Stochastic Proximal Point Methods: Stability,
Smoothness, and Deviation

Xiao-Tong Yuan and Ping Li

Cognitive Computing Lab
Baidu Research
No. 10 Xibeiwang East Road, Beijing 100193, China
10900 NE 8th St. Bellevue, Washington 98004, USA
E-mail: {xtyuan1980<EMAIL_ADDRESS>
###### Abstract

The stochastic proximal point (SPP) methods have gained recent attention for
stochastic optimization, with strong convergence guarantees and superior
robustness to the classic stochastic gradient descent (SGD) methods showcased
at little to no cost of computational overhead added. In this article, we
study a minibatch variant of SPP, namely M-SPP, for solving convex composite
risk minimization problems. The core contribution is a set of novel excess
risk bounds of M-SPP derived through the lens of algorithmic stability theory.
Particularly under smoothness and quadratic growth conditions, we show that
M-SPP with minibatch-size $n$ and iteration count $T$ enjoys an in-expectation
fast rate of convergence consisting of an
$\mathcal{O}\left(\frac{1}{T^{2}}\right)$ _bias_ decaying term and an
$\mathcal{O}\left(\frac{1}{nT}\right)$ _variance_ decaying term. In the
small-$n$-large-$T$ setting, this result substantially improves the best known
results of SPP-type approaches by revealing the impact of noise level of model
on convergence rate. In the complementary small-$T$-large-$n$ regime, we
provide a two-phase extension of M-SPP to achieve comparable convergence
rates. Moreover, we derive a near-tight high probability (over the randomness
of data) bound on the parameter estimation error of a sampling-without-
replacement variant of M-SPP. Numerical evidences are provided to support our
theoretical predictions when substantialized to Lasso and logistic regression
models.

##  1 Introduction

We consider the following problem of regularized risk minimization over a
closed convex subset $\mathcal{W}\subseteq\mathbb{R}^{p}$:

$\min_{w\in\mathcal{W}}R(w):=R^{\ell}(w)+r(w),\ \ \text{where }\
R^{\ell}(w):=\mathbb{E}_{z\sim\mathcal{D}}[\ell(w;z)],$ (1)

where $\ell:\mathcal{W}\times\mathcal{Z}\mapsto R^{+}$ is a non-negative
convex loss function whose value $\ell(w;z)$ measures the loss of a
hypothesis, parameterized by $w\in\mathcal{W}$, evaluated over a data sample
$z\in\mathcal{Z}$, $\mathcal{D}$ represents a distribution over $\mathcal{Z}$,
and $r:\mathcal{W}\mapsto\mathbb{R}^{+}$ is a data-independent non-negative
convex function whose value $r(w)$ measures certain complexity of the
hypothesis. We are particularly interested in the situation where the
composite population risk $R$ is strongly convex around its minimizers, though
in this setting the terms $R^{\ell}$ and $r$ are not necessarily required to
be so simultaneously. For an instance, the $\ell_{1}$-norm regularizer
$r(w)=\mu\|w\|_{1}$ or its grouped variants are often used for sparse
generalized linear models learning with quadratic or logistic loss functions
(Van de Geer, 2008; Ravikumar et al., 2009; Negahban et al., 2012).

In statistical machine learning, it is usually assumed that the estimator only
has access to, either as a batch training set or in an online/incremental
manner, a collection $S=\\{z_{i}\\}_{i=1}^{N}$ of i.i.d. random data instances
drawn from $\mathcal{D}$. The goal is to compute a stochastic estimator
$\hat{w}_{S}$ based on the knowledge of $S$, hopefully that it generalizes
well as a near minimizer of the population risk. More precisely, we aim at
deriving a suitable law of large numbers, i.e., a sample size vanishing rate
$\delta_{N}$ so that the _excess risk_ at $\hat{w}_{S}$ satisfies
$R(\hat{w}_{S})-R^{*}\leq\delta_{N}$ in expectation or with high probability
over $S$ where $R^{*}:=\min_{w\in\mathcal{W}}R(w)$ represents the minimal
value of composite risk.

In this work, inspired by the recent remarkable success of the stochastic
proximal point (SPP) algorithms (Patrascu and Necoara, 2017; Asi and Duchi,
2019a, b; Davis and Drusvyatskiy, 2019) and their minibatch extensions (Wang
et al., 2017b; Zhou et al., 2019; Asi et al., 2020), we provide a sharper
generalization performance analysis for a class of minibatch SPP methods for
solving the stochastic composite risk minimization problem (1).

###  1.1 Algorithm and Motivation of Study

_Minibatch Stochastic Proximal Point Algorithm._ Let
$S_{t}=\\{z_{i,t}\\}_{i=1}^{n}$ be a minibatch of $n$ i.i.d. samples drawn
from distribution $\mathcal{D}$ at time instance $t\geq 1$ and denote

$R_{S_{t}}(w):=\frac{1}{n}\sum_{i=1}^{n}\ell(w;z_{i,t})+r(w)$

as the regularized empirical risk over $S_{t}$. We consider the Minibatch
Stochastic Proximal Point (M-SPP) algorithm, as outlined in Algorithm 1, for
composite risk minimization based on a sequence of data minibatches
$S=\\{S_{t}\\}_{t=1}^{T}$. The precision value $\epsilon_{t}$ in the algorithm
quantifies the sub-optimality of $w_{t}$ for solving the inner-loop
regularized ERM over the minibatch $S_{t}$. The M-SPP algorithm is generic and
it encompasses several existing SPP methods as special cases. For example in
the extreme case when $n=1$ and $\epsilon_{t}\equiv 0$ M-SPP reduces to a
composite variant of the standard SPP method (Bertsekas, 2011), as formulated
in (5). In general, the recursion update formulation (2) can be regarded as a
natural composite extension of the existing minibatch stochastic proximal
point methods for statistical estimation (Wang et al., 2017b; Asi et al.,
2020).

Input : Regularization modulus $\\{\gamma_{t}\\}_{t\geq 1}$.

Output : $\bar{w}_{T}$ as a weighted average of $\\{w_{t}\\}_{1\leq t\leq T}$.

Initialization _Specify a value of $w_{0}$. Typically $w_{0}=0$. _

for  _$t=1,2,...,T$_ do

Sample a minibatch
$S_{t}:=\\{z_{i,t}\\}_{i=1}^{n}\overset{\text{i.i.d.}}{\sim}\mathcal{D}^{n}$
and estimate $w_{t}$ satisfying
$F_{t}(w_{t})\leq\min_{w\in\mathcal{W}}\left\\{F_{t}(w):=R_{S_{t}}(w)+\frac{\gamma_{t}}{2}\|w-w_{t-1}\|^{2}\right\\}+\epsilon_{t},$
(2) where $R_{S_{t}}(w):=\frac{1}{n}\sum_{i=1}^{n}\ell(w;z_{i,t})+r(w)$ and
$\epsilon_{t}\geq 0$ measures the sub-optimality of estimation.

end for

Algorithm 1 Minibatch Stochastic Proximal Point (M-SPP)

_Prior results and limitations._ The present study focuses on the
generalization analysis of M-SPP for convex composite risk optimization.
Recently, it has been shown by Asi et al. (2020, Theorem 2) that if the
instantaneous loss functions are strongly convex with respect to the
parameters, then the M-SPP algorithm converges at the rate of
$\mathcal{O}\left(\frac{\log(nT)}{nT}\right)$. Prior to that, Wang et al.
(2017b, Theorem 5) proved an $\mathcal{O}(\frac{1}{nT})$ rate for M-SPP when
the individual loss functions are Lipschitz continuous and strongly convex.
There results, among others for SPP (Patrascu and Necoara, 2017; Davis and
Drusvyatskiy, 2019), commonly require that each instantaneous loss should be
strongly convex which is too stringent to be fulfilled in high-dimensional or
infinite spaces. For an instance, the quadratic loss
$\ell(w;z)=\frac{1}{2}(w^{\top}x-y)^{2}$ over a feature-label pair $z=(x,y)$
is convex but in general not strongly convex, although the population risk
$R^{\ell}(w)=\frac{1}{2}\mathbb{E}(y-w^{\top}x)^{2}$ is strongly convex
provided that the covariance matrix of random feature $x$ is non-degenerate.
In the meanwhile, the Lipschitz-loss assumption made for the analysis (Wang et
al., 2017b, Theorem 5) limits its applicability to smooth losses like
quadratic loss, not to mention an interaction between Lipschitz continuity and
strong convexity (Agarwal et al., 2012; Asi and Duchi, 2019b).

The above mentioned deficiencies of prior results motivate us to investigate
the convergence behavior of M-SPP for composite risk minimization beyond the
setting where each individual loss is strongly convex and Lipschitz
continuous. From the perspective of optimization, smoothness is essential for
establishing strong convergence guarantees for solving the inner-loop strongly
convex risk minimization subproblems in (6), e.g., with variance reduced
stochastic algorithms (Johnson and Zhang, 2013; Xiao and Zhang, 2014) or
communication-efficient distributed optimization algorithms (Shamir et al.,
2014; Zhang and Lin, 2015; Yuan and Li, 2020). Aiming at covering such an
important yet less understood problem regime, we focus our study on analyzing
the convergence behavior of M-SPP when the convex loss functions are smooth
and the risk function exhibits quadratic growth property (see Assumption 2 for
a formal definition).

###  1.2 Our Contributions and Main Results

The main contribution of the present work is a sharper non-asymptotic
convergence analysis of the M-SPP algorithm through the lens of algorithmic
stability theory (Bousquet and Elisseeff, 2002; Feldman and Vondrák, 2018).
Let $W^{*}:=\\{w\in\mathcal{W}:R(w)=R^{*}\\}$ be the set of minimizers of the
composite population risk $R$. We are particularly interested in the regime
where the loss function $\ell$ is convex and smooth but not necessarily
Lipschitz (e.g., quadratic loss), while the population risk $R$ satisfies the
quadratic growth condition, i.e.,
$R(w)-R^{*}\geq\frac{\lambda}{2}\min_{w^{*}\in W^{*}}\|w-w^{*}\|^{2}$,
$\forall w\in\mathcal{W}$, for some $\lambda>0$, which can be satisfied by
strongly convex objectives, and various other statistical estimation problems
(see, e.g., Karimi et al., 2016; Drusvyatskiy and Lewis, 2018). For the family
of $L$-smooth loss functions, with $\gamma_{t}=\mathcal{O}(\lambda\rho t)$ for
an arbitrary scalar $\rho\in(0,0.5]$ and $\epsilon_{t}\equiv 0$, we show in
Theorem 1 that the excess risk at the weighted average output
$\bar{w}_{T}=\frac{2}{T(T+1)}\sum_{t=1}^{T}tw_{t}$ is in expectation upper
bounded by the following bound:

$R(\bar{w}_{T})-R^{*}\lesssim\frac{\rho\left[R(w_{0})-R^{*}\right]}{T^{2}}+\frac{LR^{*}}{\rho\lambda
nT}.$ (3)

In this composite bound, the first _bias_ component associated with initial
gap $R(w_{0})-R^{*}$ has a decaying rate
$\mathcal{O}\left(\frac{1}{T^{2}}\right)$ and the second _variance_ component
associated with $R^{*}$ converges at the rate of
$\mathcal{O}\left(\frac{1}{\lambda nT}\right)$. The variance decaying rate
actually matches the corresponding optimal rates of the SGD-type methods for
strongly convex optimization (Rakhlin et al., 2012; Dieuleveut et al., 2017;
Woodworth and Srebro, 2021). Also, such an
$\mathcal{O}\left(\frac{1}{T^{2}}+\frac{1}{\lambda nT}\right)$ bounds matches
those bounds for SPP (Davis and Drusvyatskiy, 2019) or M-SPP (Wang et al.,
2017b) which are in contrast obtained under a substantially stronger
assumption that each individual loss function should be strongly convex and
Lipschitz as well. In the realizable or near realizable machine learning
regimes where $R^{*}$ equals to or approximates zero, the variance term in (3)
would be sharper than those bounds of Wang et al. (2017b); Davis and
Drusvyatskiy (2019). To our best knowledge, the bound in (3) for smooth and
convex loss functions is _new_ to the SPP-type methods. More generally for
arbitrary convex risk functions, we present in Theorem 3 an
$\mathcal{O}(\frac{1}{\sqrt{nT}})$ excess risk bound for exact M-SPP. Further,
as shown in Theorem 4 and Theorem 5, similar results can be extended to the
inexact M-SPP given that the inner-loop sub-optimality is sufficiently small.

In the regime $T\ll n$ which is of special interest for off-line incremental
learning with large data batches, setting a near-optimal value
$\rho=\sqrt{\frac{T}{n\lambda}}$ in the excess risk bound (3) yields an
$\mathcal{O}\left(\frac{1}{T\sqrt{\lambda nT}}\right)$ rate of convergence.
This rate, in terms of $n$, is substantially slower than the
$\mathcal{O}(\frac{1}{\lambda nT})$ rate available for the previous
small-$n$-large-$T$ setup. In order to address such a deficiency, we propose a
two-phase variant of M-SSP (see Algorithm 2) to boost its performance in the
small-$T$-large-$n$ regime: in the first phase, M-SPP with sufficiently small
minibatch-size is invoked over $S_{1}$ to obtain $w_{1}$, and then initialized
by $w_{1}$ the second phase applies M-SPP to the rest minibatches. Then in
Theorem 2 we show that the in-expectation excess risk at the output of the
second phase can be accelerated to scale as

$R(\bar{w}_{T})-R^{*}\lesssim\frac{L^{2}(R(w_{0})-R^{*})}{\lambda^{2}n^{2}T^{2}}+\frac{LR^{*}}{\lambda
nT},$ (4)

which holds regardless to the mutual strength of minibatch size $n$ and
iteration count $T$.

In addition to the above in-expectation risk bounds, we further derive a high-
probability model estimation error bound of M-SPP based on algorithmic
stability theory. Our deviation analysis is carried out over a sampling-
without-replacement variant of M-SPP (see Algorithm 3). For population risk
with quadratic growth property, up to an additive term on the inner-loop sub-
optimality $\epsilon_{t}$, we establish in Theorem 6 the following deviation
bound on the estimation error $D(\bar{w}_{T},W^{*})$ that holds with
probability at least $1-\delta$ over $S$ while in expectation over the
randomness of sampling:

$D(\bar{w}_{T},W^{*})\lesssim\frac{\sqrt{L\log(1/\delta)}\log(T)}{\lambda\sqrt{nT}}+\sqrt{\frac{\left[R(w_{0})-R^{*}\right]}{\lambda
T^{2}}+\frac{LR^{*}}{\rho\lambda^{2}nT}}.$

When $T=\Omega(n)$, up to the logarithmic factors, this above bound matches
(in terms of the total sample size $N=nT$) the known minimax lower bounds for
statistical estimation even without computational limits (Tsybakov, 2008).

To highlight the core contribution of this work, the following three new
insights into M-SPP make our results distinguished from the best known of SPP-
type methods for convex optimization:

  1. 1.

First and for most, the fast rates in (3) and (4) reveal the impact of noise
level, as quantified by $R^{*}$, to convergence rate which has not been
previously known for SPP-type methods. These bounds are valid for smooth
losses and thus complement the previous ones for Lipschitz losses (Patrascu
and Necoara, 2017; Wang et al., 2017b; Davis and Drusvyatskiy, 2019).

  2. 2.

Second, the risk bounds in (3) and (4) are established under the quadratic
growth condition of population risk. This is substantially weaker than the
instantaneous-loss-wise strong convexity assumption commonly imposed by prior
analysis to achieve the comparable rates for SPP-type methods (Toulis and
Airoldi, 2017; Wang et al., 2017b; Asi et al., 2020).

  3. 3.

Third, we provide a deviation analysis of M-SPP from the viewpoint of uniform
algorithmic stability which to our best knowledge has not yet been addressed
in the previous study on SPP-type methods.

We should emphasize that, while we provide some insights into the numerical
aspects of M-SPP through an empirical study, this work is largely a
theoretical contribution.

###  1.3 Related Work

Our work is situated at the intersection of two lines of machine learning
research: stochastic optimization and algorithmic stability theory, both of
which have been actively studied with a vast body of beautiful and insightful
theoretical results established in literature. We next incompletely review
some representative work that are closely relevant to ours.

_Stochastic optimization._ Stemming from the pioneering work of Robbins and
Monro (1951), stochastic gradient descent (SGD) methods have been extensively
studied to approximately solve a simplified version of the problem (1) with
$r\equiv 0$ (Zhang, 2004; Nemirovski et al., 2009; Rakhlin et al., 2012;
Bottou et al., 2018). For the composite formulation, a vast body of proximal
SGD methods have been developed for efficient optimization in the presence of
potentially non-smooth regularizers (Hu et al., 2009; Duchi et al., 2010;
Ghadimi and Lan, 2012; Lan, 2012; Kulunchakov and Mairal, 2019). To handle the
challenges associated with stepsize selection and numerical instability of SGD
(Nemirovski et al., 2009; Bach and Moulines, 2011), a number of more
sophisticated methods including implicit stochastic/online learning (Crammer
et al., 2006; Kulis and Bartlett, 2010; Toulis et al., 2016; Toulis and
Airoldi, 2017) and stochastic proximal point (SPP) methods (Bertsekas, 2011;
Patrascu and Necoara, 2017; Asi and Duchi, 2019a, b; Davis and Drusvyatskiy,
2019) have recently been investigated for enhancing stability and adaptivity
of stochastic (composite) optimization. For an example, in our considered
composite optimization regime, the iteration procedure of vanilla SPP can be
expressed as the following recursion form for $\ i\geq 1$:

$\hat{w}^{\text{spp}}_{i}:=\mathop{\arg\min}_{w\in\mathcal{W}}\ell(w;z_{i})+r(w)+\frac{\gamma_{i}}{2}\|w-\hat{w}^{\text{spp}}_{i-1}\|^{2},$
(5)

where $z_{i}\sim\mathcal{D}$ is a random data sample, $\gamma_{i}$ is a
regularization modulus and $\|\cdot\|$ stands for the Euclidean norm. In
contrast to standard SGD methods which are simple in per-iteration modeling
but brittle to stepsize choice, the SPP methods are more accurate in objective
approximation which leads to substantially improved stability to the choice of
algorithm hyper-parameters while enjoying optimal guarantees on convergence
(Asi and Duchi, 2019a, b).

An attractive feature of these above (proximal) stochastic optimization
methods is that their convergence guarantees directly apply to the population
risk and the minimax optimal rates of order $\mathcal{O}(\frac{1}{T})$ are
achievable after $T$ rounds of iteration for strongly convex problems
(Nemirovski et al., 2009; Agarwal et al., 2012; Rakhlin et al., 2012). For
large-scale machine learning, the improved memory efficiency is another
practical argument in favor of stochastic over batch optimization methods.
However, due to the sequential processing nature, the stochastic optimization
methods tend to be less efficient for parallelization especially in
distributed computing environment where excessive communication between nodes
would be required for model update (Bottou et al., 2018).

_Empirical risk minimization._ At the opposite end of SGD-type and online
learning, the following defined (composite) empirical risk minimization (ERM,
a.k.a., M-estimation) is another popularly studied formulation for statistical
learning (Lehmann and Casella, 2006):

$\hat{w}^{\text{erm}}_{S}:=\mathop{\arg\min}_{w\in\mathcal{W}}\left\\{R_{S}(w):=\frac{1}{N}\sum_{i=1}^{N}\ell(w;z_{i})+r(w)\right\\}.$

Thanks to the finite-sum structure, a large body of randomized incremental
algorithms with linear rates of convergence have been established for ERM
including SVRG (Johnson and Zhang, 2013; Xiao and Zhang, 2014), SAGA (Defazio
et al., 2014) and Katyusha (Allen-Zhu, 2017), to name a few. From the
perspective of distributed computation, one intrinsic advantage of ERM over
SGD-type methods lies in that it can better explore the statistical
correlation among data samples for designing communication-efficient
distributed optimization algorithms (Jaggi et al., 2014; Shamir et al., 2014;
Zhang and Lin, 2015; Lee et al., 2017). Unlike stochastic optimization
methods, the generalization performances of the batch or incremental
algorithms are by nature controlled by that of ERM (Bottou and Bousquet, 2007)
which has long been studied with a bunch of insightful results available
(Vapnik, 1999; Bartlett et al., 2005; Srebro et al., 2010; Mei et al., 2018).
Particularly for strongly convex risk functions, the
$\mathcal{O}(\frac{1}{N})$ rate of convergence is possible for ERM (Bartlett
et al., 2005; Koltchinskii, 2006; Zhang et al., 2017), though these fast rates
are in general dimensionality-dependent for parametric learning models.

It has been recognized that SGD-type and ERM-type approaches cannot dominate
each other in terms of generalization, runtime, storage and parallelization
efficiency. This motivates a recent trend of trying to propose the so called
stochastic model-based methods that can achieve the best of two worlds. Among
others, a popular paradigm for such a purpose of combination is _minibatch
proximal update_ which in each iteration updates the model via (approximately)
solving a local ERM over a stochastic minibatch (Li et al., 2014; Wang et al.,
2017b; Asi et al., 2020; Deng and Gao, 2021). This strategy can be viewed as a
minibatch extension to the SPP algorithm and it has been shown to attain a
substantially improved trade-off between computation, communication and memory
efficiency for large-scale distributed machine learning (Li et al., 2014; Wang
et al., 2017a). Alternatively, a number of online extensions of the
incremental finite-sum algorithms, such as streaming SVRG (Frostig et al.,
2015) and streaming SAGA (Jothimurugesan et al., 2018), have been proposed for
stochastic optimization with competitive guarantees to ERM but at lower cost
of computation.

_Algorithmic stability and generalization._ Since the seminal work of Bousquet
and Elisseeff (2002), algorithmic stability has been extensively studied with
remarkable success achieved in establishing generalization bounds for strongly
convex ERM estimators (Zhang, 2003; Mukherjee et al., 2006; Shalev-Shwartz et
al., 2010). Particularly, the state-of-the-art risk bounds of strongly convex
ERM are offered by approaches based on the notion of uniform stability
(Feldman and Vondrák, 2018, 2019; Bousquet et al., 2020; Klochkov and
Zhivotovskiy, 2021). It was shown by Hardt et al. (2016) that the solution
obtained via (stochastic) gradient descent is stable for smooth convex or non-
convex loss functions. For non-smooth convex losses, the stability induced
generalization bounds of SGD have been established in expectation (Lei and
Ying, 2020) or deviation (Bassily et al., 2020). For learning with sparsity,
algorithmic stability theory has been employed to derive the generalization
bounds of the popularly used iterative hard thresholding (IHT) algorithm (Yuan
and Li, 2022). Through the lens of uniform algorithmic stability, convergence
rates of M-SPP have been studied for convex (Wang et al., 2017b) and weakly
convex (Deng and Gao, 2021) Lipschitz losses. While sharing a similar spirit
to Wang et al. (2017b); Deng and Gao (2021), our analysis customized for
smooth convex loss functions is considerably different and the resultant fast
rates are of special interest in low-noise statistical settings (Srebro et
al., 2010).

###  1.4 Notation and Paper Organization

_Notation._ The key quantities and notations frequently used in our analysis
are summarized in Table 1.

Notation | Definition
---|---
$n$ | minibatch size
$T$ | round of iteration
$N$ | total number of samples visited, i.e., $N=nT$
$f$ | hypothesis
$\ell$ | loss function
$r$ | regularization term
$R^{\ell}$ | population risk: $R^{\ell}(w):=\mathbb{E}_{(x,y)\sim\mathcal{D}}[\ell(f_{w}(x),y)]$
$R$ | composite population risk: $R(w):=R^{\ell}(w)+r(w)$
$R^{*}$ | the optimal value of composite risk, i.e., $R^{*}:=\min_{w\in\mathcal{W}}R(w)$
$W^{*}$ | the optimal solution set of composite risk, i.e., $W^{*}:=\mathop{\arg\min}_{w\in\mathcal{W}}R(w)$
$S_{t}$ | data minibatch at time instance $t$
$S_{I}$ | The union of data minibatch over $I$, i.e., $S_{I}:=\\{S_{t}\\}_{t\in I}$
$R^{\ell}_{S}$ | empirical risk over $S$, i.e., $R^{\ell}_{S}(w):=\frac{1}{|S|}\sum_{(x,y)\in S}\ell(f_{w}(x,y)$
$R_{S}$ | composite empirical risk over $S$, i.e., $R_{S}(w):=R^{\ell}_{S}(w)+r(w)$
$\epsilon_{t}$ | precision of minibatch risk minimization at time instance $t$
$\|w\|_{1}$ |  $\ell_{1}$-norm of a vector $w$, i.e., $\|w\|_{1}:=\sum_{i}|[w]_{i}|$
$\|w\|$ | Euclidean norm of a vector $w$
$D(w,W^{*})$ | the distance from $w$ to $W^{*}$, i.e., $D(w,W^{*})=\min_{w^{*}\in W^{*}}\|w-w^{*}\|$
$[T]$ | $[T]:=\\{1,...,T\\}$
$\mathbf{1}_{\\{C\\}}$ | the indicator function of the condition $C$
Table 1: Table of notation.

_Organization._ The paper proceeds with the material organized as follows: In
Section 2, we analyze the risk bounds of exact M-SPP with convex and smooth
loss functions and present a two-phase variant to further improve convergence
performance. In Section 3, we extend our analysis to the more realistic
setting where inexact M-SPP iteration is allowed. In Section 4, we study the
high-probability estimation error bounds of M-SPP. A comprehensive comparison
to some closely relevant results is highlighted in Section 5. The numerical
study for theory verification and algorithm evaluation is provided in Section
6. The concluding remarks are made in Section 7. All the proofs of main
results and some additional results on the iteration stability of M-SPP are
relegated to appendix.

##  2 A Sharper Analysis of M-SPP for Smooth Loss

In this section, we analyze the convergence rate of M-SPP for smooth and
convex loss functions using the tools developed in algorithmic stability
theory. In what follows, for the sake of notation simplicity and presentation
clarity of core ideas, we assume for the time being that the inner-loop
composite ERM in the M-SPP iteration procedure (2) has been solved exactly
with $\epsilon_{t}\equiv 0$, i.e.,

$w_{t}=\mathop{\arg\min}_{w\in\mathcal{W}}\left\\{F_{t}(w):=R_{S_{t}}(w)+\frac{\gamma_{t}}{2}\|w-w_{t-1}\|^{2}\right\\}.$
(6)

A full convergence analysis for the inexact variant (i.e., $\epsilon_{t}>0$)
will be presented in the Section 3 via a slightly more involved perturbation
analysis.

###  2.1 Basic Assumptions

We begin by introducing some basic assumptions that will be used in the
analysis to follow. We say a differentiable function
$g:\mathcal{W}\mapsto\mathbb{R}$ is $L$-smooth if $\forall s,t\in\mathbb{R}$,

$\left|g(w)-g(w^{\prime})-\langle\nabla
g(w),w-w^{\prime}\rangle\right|\leq\frac{L}{2}|w-w^{\prime}|^{2}.$

As formally stated in the following assumption, we suppose that the individual
loss functions are convex and $L$-smooth which can be satisfied, e.g., by the
quadratic loss (for regression) and the logistic loss (for prediction).

######  Assumption 1.

The loss function $\ell$ is convex and $L$-smooth with respect to its first
argument. Also, we assume that the regularization term $r$ is convex over
$\mathcal{W}$.

Let us define $D(w,W^{*}):=\min_{w^{*}\in W^{*}}\|w-w^{*}\|$ as the distance
from $w$ to the set $W^{*}$ of minimizers. The next assumption requires that
the population risk has the characterization of quadratic growth away from the
set of minimizers (Anitescu, 2000; Karimi et al., 2016).

######  Assumption 2.

The population risk function $R$ satisfies $R(w)\geq
R^{*}+\frac{\lambda}{2}D^{2}(w,W^{*}),\forall w\in\mathcal{W}$ for some
$\lambda>0$.

Clearly, the quadratic growth property can be implied by the traditional
strong convexity condition (around the minimizers) which is satisfied by a
number of popular learning models including linear and logistic regression,
generalized linear models, smoothed Huber losses, and various other
statistical estimation problems. Particularly, Assumption 2 holds when
$R^{\ell}$ is strongly convex and $r$ is convex. Notice that for risk
functions with quadratic growth property, the prior analysis of M-SPP for
Lipschitz losses (Wang et al., 2017b) is not generally applicable because
Assumption 2 implies that the Lipschitz constant of loss could be arbitrarily
large if the infinite distance $\min_{w^{*}\in
W^{*}}\|w-w^{*}\|\rightarrow\infty$ is allowed.

###  2.2 Main Results

The following theorem is our main result on the in-expectation rate of
convergence of the exact M-SPP with smooth loss and quadratic growth
population risk functions. Recall that $N=nT$ is the total number of data
points visited up to the iteration counter $T$.

######  Theorem 1.

Suppose that Assumptions 1 and 2 hold. Consider $\epsilon_{t}\equiv 0$ and the
weighted average output $\bar{w}_{T}=\frac{2}{T(T+1)}\sum_{t=1}^{T}tw_{t}$ in
Algorithm 1. Let $\rho\in(0,0.5]$ be an arbitrary scalar.

  * (a)

Suppose that $n\geq\frac{64L}{\lambda\rho}$. Set $\gamma_{t}=\frac{\lambda\rho
t}{4}$ for $t\geq 1$. Then for any $T\geq 1$,

$\mathbb{E}\left[R(\bar{w}_{T})-R^{*}\right]\leq\frac{4\rho\left[R(w_{0})-R^{*}\right]}{T^{2}}+\frac{2^{9}L}{\lambda\rho
nT}R^{*}.$

  * (b)

Set $\gamma_{t}=\frac{\lambda\rho t}{4}+\frac{16L}{n}$ for $t\geq 1$. Then for
any $T\geq 1$,

$\mathbb{E}\left[R(\bar{w}_{T})-R^{*}\right]\leq\left(\frac{4\rho}{T^{2}}+\frac{2^{8}L}{\lambda
nT}\right)[R(w_{0})-R^{*}]+\left(\frac{2^{16}L^{2}}{\lambda^{2}\rho^{2}n^{2}T}+\frac{2^{9}L}{\lambda\rho
nT}\right)R^{*},$

###### Proof.

The proof technique is inspired by the uniform stability arguments developed
by Wang et al. (2017b) for Lipschitz and instance-wise strongly convex loss,
with several new elements along developed for handling smooth loss and
quadratic growth of risk function. As a non-trivial ingredient, we show that
it is possible to extend those stability arguments to smooth losses in view of
a classical result from Srebro et al. (2010, Lemma 2.1) that allows the
derivative of a smooth loss to be bounded in terms of its function value. See
Appendix A.1 for a full proof of this result. ∎

A few remarks on Theorem 1 are in order.

######  Remark 1.

In Part (a), the minibatch size is required to be sufficiently large. In this
setting, the excess risk bound consists of two components: the first _bias_
component associated with initial gap $R(w_{0})-R^{*}$ has a decaying rate
$\mathcal{O}(\frac{1}{T^{2}})$ and the second _variance_ component associated
with $R^{*}$ vanishes at a dominate rate of $\mathcal{O}(\frac{1}{\lambda
nT})$. The variance term shows that the convergence rate can be improved in
the low-noise settings where the factor of $R^{*}$ is relatively small.
Extremely in the separable case with $R^{*}=0$, the excess risk bound of
Theorem 1 would scale as fast as $\mathcal{O}(\frac{1}{T^{2}})$.

######  Remark 2.

One disadvantage of the result in Part (a) lie in that the minibatch size is
required to be sufficiently larger than the condition number of the population
risk $R$. Contrastively, the excess risk bound in Part (b) holds for arbitrary
minibatch sizes. The cost, however, is a relatively slower bias decaying term
$\mathcal{O}(\frac{1}{T^{2}}+\frac{1}{\lambda nT})$ which is dominated by
$\mathcal{O}(\frac{1}{\lambda nT})$ in the case of $T\gg n$.

######  Remark 3.

Let $N=nT$ be the total number of data points accessed. When $T\gg n$, the
$\mathcal{O}(\frac{1}{N})$ dominant rates in Theorem 1 match those prior ones
for SPP-type methods (Wang et al., 2017b; Davis and Drusvyatskiy, 2019) which
are, however, obtained under the assumption that each individual loss function
should be Lipschitz continuous and strongly convex. In comparison to the
$\mathcal{O}(\frac{1}{N})$ rate established for SGD with smooth loss (Lei and
Ying, 2020, Theorem 12), our result in Theorem 1 is stronger and less
stringent in the following senses: 1) our bound shows explicitly the impact of
$R^{*}$ which usually represents the noise level of model, and 2) we only
require the population risk to have quadratic growth property while the bound
of Lei and Ying (2020, Theorem 12) not only requires the loss to be Lipschitz
but also assumes the empirical risk to be strongly convex.

Let us further look into the choice of the scalar $\rho$ in Theorem 1. We
focus the discussion on the part (a) and similar observations apply to the
part (b). We distinguish the discussion in the following two complementary
cases regarding the mutual strength of minibatch-size $n$ and round of
iteration $T$:

  * •

Case I: Small-$n$-large-$T$. Suppose that $n=\mathcal{O}(1)$ and
$T\rightarrow\infty$ is allowed. In this case, simply setting $\rho=0.5$
yields the convergence rate of order
$\mathcal{O}\left(\frac{1}{T^{2}}+\frac{1}{\lambda nT}\right)$ in the part
(a).

  * •

Case II: Small-$T$-large-$n$. Suppose that $T=\mathcal{O}(1)$ and
$n\rightarrow\infty$ is allowed. In this setup, given that
$n\geq\frac{4T}{\lambda}$, then with a roughly optimal choice
$\rho=\sqrt{\frac{T}{n\lambda}}$ the excess risk bound in Theorem 1(a) will be
of the order $\mathcal{O}\left(\frac{1}{T\sqrt{\lambda nT}}\right)$, which is
substantially slower than the previous fast rate in Case I. This is intuitive
because M-SPP with large minibatches behaves more like regularized ERM which
is known to exhibit slow rate of convergence even for strongly convex problems
(Shalev-Shwartz et al., 2010; Srebro et al., 2010). Nevertheless, such a
small-$T$-large-$n$ setup is of special interest for off-line incremental
learning with large minibatches and distributed statistical learning (Li et
al., 2014; Wang et al., 2017b; You et al., 2020). We will address this
critical case in the next subsection.

###  2.3 A Two-Phase M-SPP Method

Input : Dataset $S=\\{S_{t}\\}_{t=1}^{T}$ in which
$S_{t}:=\\{z_{i,t}\\}_{i=1}^{n}\overset{\text{i.i.d.}}{\sim}\mathcal{D}^{n}$,
regularization modulus $\\{\gamma_{t}>0\\}_{t\in[T]}$.

Output : $\bar{w}_{T}$ as a weighted average of $\\{w_{t}\\}_{2\leq t\leq T}$.

Initialization _Specify a value of $w_{0}$. Typically $w_{0}=0$._

/* Phase-I */

Divide sample $S_{1}$ into disjoint minibatches of equal size $m$;

Run M-SPP over these minibatches to obtain the output $w_{1}$;

/* Phase-II */

Initialized with $w_{1}$, run M-SPP over data minibatches $\\{S_{t}\\}_{2\leq
t\leq T}$ with $\\{\gamma_{t}\\}_{2\leq t\leq T}$ to obtain the sequence
$\\{w_{t}\\}_{2\leq t\leq T}$.

Algorithm 2 Two-Phase M-SPP (M-SPP-TP)

To remedy the deficiencies mentioned in the previous discussion, we propose a
two-phase variant of M-SSP, as outlined in Algorithm 2, to boost its
performance in the small-$T$-large-$n$ regimes. The procedure can be regarded
as sort of a restarting argument (Nemirovskii and Nesterov, 1985; Renegar and
Grimmer, 2022; Zhou et al., 2022) for M-SPP. More specifically, the Phase-I
serves as an initialization step that invokes M-SPP to a uniform division of
$S_{1}$ with minibatch size $m$ to obtain $w_{1}$. Then starting from $w_{1}$,
the Phase-II just invokes M-SPP to the consequent large minibatches
$\\{S_{t}\\}_{t\geq 2}$ which is suitable for large-scale parallelization if
applicable. The following theorem is a consequence of Theorem 1 to such a two-
phase M-SPP procedure.

######  Theorem 2.

Suppose that Assumptions 1 and 2 hold. Consider $\epsilon_{t}\equiv 0$ for
implementing M-SPP in both Phase-I and Phase-II of Algorithm 2. Consider the
weighted average output $\bar{w}_{T}=\frac{2}{(T-1)(T+2)}\sum_{t=2}^{T}tw_{t}$
in Phase-II.

  * (a)

Suppose that $n\geq\frac{128L}{\lambda}$. Set $m=\frac{128L}{\lambda}$ in
Phase-I and $\gamma_{t}=\frac{\lambda t}{8}$ for implementing M-SPP in both
Phase-I and Phase II. Then for any $T\geq 2$, $\bar{w}_{T}$ satisfies

$\mathbb{E}\left[R(\bar{w}_{T})-R^{*}\right]\lesssim\frac{L^{2}\left[R(w_{0})-R^{*}\right]}{\lambda^{2}n^{2}T^{2}}+\frac{L}{\lambda
nT}R^{*}.$

  * (b)

Set $m=\mathcal{O}(1)$ in Phase-I and $\gamma_{t}=\frac{\lambda
t}{8}+\frac{16L}{n}$ for implementing M-SPP in both Phase-I and Phase-II. Then
for any $T\geq 2$, $\bar{w}_{T}$ satisfies

$\mathbb{E}\left[R(\bar{w}_{T})-R^{*}\right]\lesssim\frac{L^{2}\left[R(w_{0})-R^{*}\right]}{\lambda^{2}nT}+\frac{L^{3}}{\lambda^{3}nT}R^{*}.$

###### Proof.

See Appendix A.2 for a proof of this result. ∎

######  Remark 4.

The part (a) of Theorem 2 suggests that when the minibatch size is
sufficiently large, the excess risk bound of two-phase M-SPP has a bias
decaying term of scale $\mathcal{O}\left(\frac{1}{n^{2}T^{2}}\right)$ and a
variance term that decays at the rate of $\mathcal{O}(\frac{1}{nT})$. The rate
is valid even when the scales of $T$ relatively small, and thus is stronger
than the $\mathcal{O}\left(\frac{1}{T\sqrt{nT}}\right)$ rate implied by
Theorem 1 for the vanilla M-SPP in the small-$T$-large-$n$ regime. It is worth
to mention that both the bias and variance components in our bound for M-SPP
are faster than those derived for strongly convex ERM (Srebro et al., 2010).

######  Remark 5.

The excess risk bound in Part (b) of Theorem 2 is valid for arbitrary
minibatch sizes, but at the cost of a relatively slower
$\mathcal{O}(\frac{1}{nT})$ bias decaying rate.

###  2.4 Results for Arbitrary Convex Risks

We further analyze the proposed M-SPP algorithm when the loss function $\ell$
is convex and smooth, but without requiring that the composite risk $R$ has
quadratic growth property. The following is our main result in such a generic
setting.

######  Theorem 3.

Suppose that Assumption 1 holds. Set
$\gamma_{t}\equiv\gamma\geq\frac{16L}{n}$. Let
$\bar{w}_{T}=\frac{1}{T}\sum_{t=1}^{T}w_{t}$ be the average output of
Algorithm 1. Then

$\mathbb{E}\left[R(\bar{w}_{T})-R^{*}\right]\lesssim\frac{\gamma}{T}D^{2}(w_{0},W^{*})+\frac{L}{\gamma
n}R^{*}.$

Particularly for $\gamma=\sqrt{\frac{T}{n}}+\frac{16L}{n}$, it holds that

$\mathbb{E}\left[R(\bar{w}_{T})-R^{*}\right]\lesssim\left(\frac{1}{\sqrt{nT}}+\frac{L}{nT}\right)D^{2}(w_{0},W^{*})+\frac{L}{\sqrt{nT}}R^{*}.$

###### Proof.

See Appendix A.3 for a proof of this result. ∎

######  Remark 6.

The first bound of Theorem 3 implies that for any $\epsilon\in(0,1)$, by
setting $\gamma=\mathcal{O}\left(\frac{L}{\epsilon n}\right)$,
$R(\bar{w}_{t})$ converges to $(1+\epsilon)R^{*}$ at the rate of
$\mathcal{O}(\frac{1}{nT\epsilon})$. This bound matches the results of Lei and
Ying (2020, Theorem 4) for smooth SGD method. The second bound of Theorem 3
further shows that by setting
$\gamma=\mathcal{O}(\sqrt{\frac{T}{n}}+\frac{L}{n})$, the excess risk of
$\bar{w}_{T}$ decays at the rate of $\mathcal{O}(\frac{1}{\sqrt{nT}})$ for
both bias and variance terms, which matches in order the corresponding bound
derived for Lipschitz-loss (Wang et al., 2017b, Theorem 4). To our knowledge,
such a bias-variance composite rate of convergence is new for SPP-type methods
with convex and smooth loss functions.

Analogous to the robustness analysis of SPP (Asi and Duchi, 2019a, b), we have
also analyzed the iteration stability of M-SPP for convex losses with respect
to the choice of regularization modulus $\gamma_{t}$. The corresponding
results, which can be found in Appendix A.4, confirm that the choice of
$\gamma_{t}$ is insensitive to the gradient scale of loss functions for
generating a non-divergent sequence of estimation errors.

##  3 Perturbation Analysis for Inexact M-SPP

In the preceding section, we have analyzed the convergence rates of M-SPP
under the assumption that the inner-loop proximal ERM subproblems constructed
in its iteration procedure (2) are solved exactly, i.e., $\epsilon_{t}\equiv
0$. To make our analysis more practical, we further provide in this section a
perturbation analysis of M-SPP when the inner-loop proximal ERM subproblems
are only required to be solved approximately up to certain precision
$\epsilon_{t}>0$. As a starting point, we need to impose the following
Lipschitz continuity assumption on the regularization term $r$.

######  Assumption 3.

The regularization term $r$ is Lipschitz continuous over $\mathcal{W}$, i.e.,
$|r(w)-r(w^{\prime})|\leq G\|w-w^{\prime}\|$, $\forall
w,w^{\prime}\in\mathcal{W}$.

For example, the $\ell_{1}$-norm regularizer $r(w)=\mu\|w\|_{1}$ satisfies
this assumption with respect to Euclidean norm as
$|r(w)-r(w^{\prime})|=\mu|\|w\|_{1}-\|w^{\prime}\|_{1}|\leq\mu\|w-w^{\prime}\|_{1}\leq\mu\sqrt{p}\|w-w^{\prime}\|$.

The following theorem is our main result on the rate of convergence of the
inexact M-SPP for composite stochastic convex optimization with smooth losses.

######  Theorem 4.

Suppose that Assumptions 1, 2 and 3 hold. Let $\rho\in(0,1/4]$ be an arbitrary
scalar and set $\gamma_{t}=\frac{\lambda\rho t}{4}$. Suppose that
$n\geq\frac{76L}{\lambda\rho}$. Assume that
$\epsilon_{t}\leq\frac{\epsilon}{nt^{4}}$ for some $\epsilon\in[0,1]$. Then
for any $T\geq 1$, the weighted average output
$\bar{w}_{t}=\frac{2}{T(T+1)}\sum_{t=1}^{T}tw_{t}$ of Algorithm 1 satisfies

$\mathbb{E}\left[R(\bar{w}_{t})-R^{*}\right]\lesssim\frac{\rho}{T^{2}}(R(w_{0})-R^{*})+\frac{L}{\lambda\rho
nT}R^{*}+{\frac{\sqrt{\epsilon}}{T^{2}}\left(\frac{L}{\lambda\rho}+G\sqrt{\frac{1}{\lambda\rho}}\right)}.$

###### Proof.

See Appendix B.1 for a proof of this result. We would like to highlight that
our perturbation analysis for smooth loss is considerably different from that
of Wang et al. (2017b) developed for Lipschitz loss. This is mainly because in
the smooth loss case, the change of loss could no longer be upper bounded by
the change of prediction, and thus we need to make a more careful treatment to
the perturbation caused by inexact minimization of the regularized minibatch
empirical risk. ∎

We provide in order a few remarks on Theorem 4.

######  Remark 7.

Theorem 4 suggests that the excess risk bound of exact M-SPP in the part (a)
of Theorem 1 can be extended to its inexact version, provided that the inner-
loop minibatch ERMs (2) are solved to sufficient accuracy, say,
$\epsilon_{t}\leq\mathcal{O}\left(\frac{1}{nt^{4}}\right)$. Similarly, the
result in the part (b) of Theorem 1 for arbitrary minibatch sizes can also be
extended to the inexact M-SPP, which is omitted to avoid redundancy. Since the
inner-loop minibatch ERMs are strongly convex and the loss functions are
smooth, in average the desired accuracy can be attained in logarithmic time
$\mathcal{O}\left(\log\left(\frac{1}{\epsilon_{t}}\right)\right)$ via
variance-reduced SGD methods (Xiao and Zhang, 2014).

######  Remark 8.

Analogous to the discussions at the end of Section 2.2, by specifying the
choice of $\rho$ we can derive a direct consequent result of Theorem 4 which
more explicitly shows the rate of convergence with respect to $N=nT$. Also for
the two-phase M-SPP, in view of Theorem 4 we can show that the bound in
Theorem 2 can be extended to the inexact setting if the minibatch optimization
is sufficiently accurate. These extensions are more or less straightforward
and thus are omitted.

In the following theorem, we provide an excess risk bound for the inexact
M-SPP when the composite risk $R$ is convex but not necessarily has quadratic
growth property.

######  Theorem 5.

Suppose that Assumptions 1 and 3 hold. Set
$\gamma_{t}\equiv\gamma\geq\frac{19L}{n}$. Assume that
$\epsilon_{t}\leq\min\left\\{\frac{\epsilon}{n^{2}t^{5}},\frac{2G^{2}}{9n^{2}\gamma}\right\\}$
for some $\epsilon\in[0,1]$. Then the average output
$\bar{w}_{T}=\frac{1}{T}\sum_{t=1}^{T}w_{t}$ of Algorithm 1 satisfies

$\mathbb{E}\left[R(\bar{w}_{T})-R^{*}\right]\lesssim\frac{\gamma}{T}D^{2}(w_{0},W^{*})+\frac{L}{\gamma
n}R^{*}+{\left(\frac{L}{\gamma
n}+\frac{\gamma}{LnT}+\frac{G}{\sqrt{\gamma}nT}\right)\sqrt{\epsilon}}.$

Particularly for $\gamma=\sqrt{\frac{T}{n}}+\frac{19L}{n}$, it holds that

$\mathbb{E}\left[R(\bar{w}_{T})-R^{*}\right]\lesssim\left(\frac{1}{\sqrt{nT}}+\frac{L}{nT}\right)D^{2}(w_{0},W^{*})+\frac{L}{\sqrt{nT}}R^{*}+{\left(\frac{L+G}{\sqrt{nT}}+\frac{1}{nT}\right)\sqrt{\epsilon}}.$

###### Proof.

See Appendix B.2 for a proof of this result. ∎

######  Remark 9.

Theorem 5 confirms that the excess risk bounds established in Theorem 3 for
exact M-SPP are tolerant to sufficiently small sub-optimality
$\epsilon_{t}\leq\mathcal{O}(\frac{1}{n^{2}t^{5}})$ of minibatch proximal ERM
subproblems.

##  4 Performance Guarantees with High Probability

In the previous two sections, we have analyzed the excess risk bounds of M-SPP
in expectation. In this section, we move on to study high-probability
guarantees of M-SPP with respect to the randomness of training data, still
under the notion of algorithmic stability. To this end, we first introduce a
variant of M-SPP which carries out the proximal point update via sampling
without replacement over the given data minibatches. We then show that the
output of the proposed algorithm is uniformly stable in expectation over the
randomness of sampling. As a main result of this section, for strongly convex
population risk, we establish a near-optimal high probability (with respect to
data) bound on the estimation error $\|\bar{w}_{t}-w^{*}\|$ that holds in
expectation over the randomness of inner-data sampling. Additionally, we
provide a high-probability generalization bound for arbitrary convex loss.

###  4.1 Sampling Without Replacement M-SPP

Input : Dataset $S=\\{S_{t}\\}_{t=1}^{T}$ in which
$S_{t}:=\\{z_{i,t}\\}_{i=1}^{n}\overset{\text{i.i.d.}}{\sim}\mathcal{D}^{n}$,
regularization modulus $\\{\gamma_{t}>0\\}_{t\in[T]}$.

Output : $\bar{w}_{T}$ as a weighted average of $\\{w_{t}\\}_{1\leq t\leq
T}$..

Initialization _Specify a value of $w_{0}$. Typically $w_{0}=0$._

for  _$t=1,2,...,T$_ do

Uniformly randomly sample an index $\xi_{t}\in[T]$ without replacement.

Estimate $w_{t}$ satisfying
$F_{t}(w_{t})\leq\min_{w\in\mathcal{W}}\left\\{F_{t}(w):=R_{S_{\xi_{t}}}(w)+\frac{\gamma_{t}}{2}\|w-w_{t-1}\|^{2}\right\\}+\epsilon_{t},$
(7) where $\epsilon_{t}\geq 0$ measures the sub-optimality.

end for

Algorithm 3 Sampling Without Replacement M-SPP (M-SPP-SWoR)

Let us consider the M-SPP-SWoR (M-SPP via Sampling Without Replacement)
procedure as outlined in Algorithm 3. Given a set $S$ of $T$ data minibatches,
at each iteration, the algorithm uniformly randomly samples one minibatch from
$S$ without replacement for proximal update. After $T$ rounds of iteration,
all the minibatches are used to update the model. Since this procedure is
merely a random shuffling variant of M-SPP as presented in Algorithm 1, we can
see that all the in-expectation bounds established in the previous sections
for M-SPP directly transfer to M-SPP-SWoR under any implementation of
shuffling. As we will show shortly in the next subsection that such a random
shuffling scheme is beneficial for boosting the on-average algorithmic
stability of M-SPP which then leads to strong high-probability guarantees for
M-SPP-SWoR.

###  4.2 A Uniform Stability Analysis

Let $S=\\{S_{t}\\}_{t\in[T]}$ and $S^{\prime}=\\{S^{\prime}_{t}\\}_{t\in[T]}$
be two sets of data minibatches. We denote by $S_{t}\doteq S^{\prime}_{t}$ if
$S_{t}$ and $S^{\prime}_{t}$ differ in a single data point, and by $S\doteq
S^{\prime}$ if $S$ and $S^{\prime}$ differ in a single minibatch and a single
data point in that minibatch. We introduce the following concept of uniform
stability of M-SPP which substantializes the concept of uniform algorithmic
stability that serves as a powerful tool for analyzing generalization bounds
of statistical estimators and their learning algorithms (Bousquet and
Elisseeff, 2002; Hardt et al., 2016; Feldman and Vondrák, 2019).

######  Definition 1 (Uniform Stability of M-SPP).

The M-SPP algorithm is said to be $\varrho$-uniformly stable with respect to a
mapping $h:\mathcal{W}\mapsto\mathbb{R}^{q}$ if
$\left\|h(\bar{w}_{T})-h(\bar{w}^{\prime}_{T})\right\|\leq\varrho$ for any
pair of data sets $S\doteq S^{\prime}$.

The following result gives a uniform stability (with respect to identical
mapping) bound of the vanilla M-SPP (Algorithm 1) that holds
deterministically, and a corresponding bound for M-SPP-SWoR (Algorithm 3) that
holds in expectation over the randomness of minibatch sampling.

######  Proposition 1.

Suppose that Assumption 1 holds and the loss function is bounded such that
$0\leq\ell(y^{\prime},y)\leq M$ for all $y,y^{\prime}$. Let
$S=\\{S_{t}\\}_{t\in[T]}$ and $S^{\prime}=\\{S^{\prime}_{t}\\}_{t\in[T]}$ be
two sets of data minibatches satisfying $S\doteq S^{\prime}$. Then

  * (a)

The weighted average output $\bar{w}_{T}$ and $\bar{w}^{\prime}_{T}$
respectively generated by M-SPP (Algorithm 1) over $S$ and $S^{\prime}$
satisfy

$\sup_{S,S^{\prime}}\|\bar{w}_{T}-\bar{w}^{\prime}_{T}\|\leq\frac{4\sqrt{2LM}}{n\min_{t\in[T]}\gamma_{t}}+\sum_{t=1}^{T}2\sqrt{\frac{2\epsilon_{t}}{\gamma_{t}}}.$

  * (b)

The weighted average output $\bar{w}_{T}$ and $\bar{w}^{\prime}_{T}$
respectively generated by M-SPP-SWoR (Algorithm 3) over $S$ and $S^{\prime}$
satisfy

$\sup_{S,S^{\prime}}\mathbb{E}_{\xi_{[T]}}\left[\|\bar{w}_{T}-\bar{w}^{\prime}_{T}\|\right]\leq\sum_{t=1}^{T}\left\\{\frac{4\sqrt{2LM}}{nT\gamma_{t}}+2\sqrt{\frac{2\epsilon_{t}}{\gamma_{t}}}\right\\}.$

###### Proof.

See Appendix C.1 for a proof. ∎

######  Remark 10.

Suppose that the sub-optimality $\\{\epsilon_{t}\\}_{t\in[T]}$ are
sufficiently small. If setting $\gamma_{t}=\mathcal{O}(t)$ as used for
population risks with quadratic growth property, then Proposition 1 shows that
M-SPP is $\mathcal{O}\big{(}\frac{1}{n}\big{)}$-uniformly stable, while in
expectation over the randomness of without-replacement sampling, M-SPP-SWoR
has an much improved uniform stability parameter scaling as
$\mathcal{O}\big{(}\frac{\log(T)}{nT}\big{)}$. If setting
$\gamma_{t}\equiv\sqrt{\frac{T}{n}}$ as used for generic convex loss, then
M-SPP will be $\mathcal{O}\big{(}\frac{1}{\sqrt{nT}}\big{)}$-uniformly stable
while M-SPP-SWoR has an identical uniform stability parameter in expectation
over sampling.

In the following theorem, based on the uniform stability bounds in Proposition
1, we derive an upper bound on the estimation error $D(\bar{w}_{T},W^{*})$ of
M-SPP-SWoR that holds with high probability over data distribution while in
expectation over randomly sampling the minibatches for update.

######  Theorem 6.

Suppose that Assumptions 1, 2, 3 hold and the loss function $\ell$ is bounded
in the interval $(0,M]$. Let $\rho\in(0,1/4]$ be an arbitrary scalar and set
$\gamma_{t}=\frac{\lambda\rho t}{4}$. Suppose that
$n\geq\frac{76L}{\lambda\rho}$. Assume that
$\epsilon_{t}\leq\min\left\\{\frac{\epsilon}{nt^{4}},\frac{LM}{\lambda\rho
n^{2}T^{2}t}\right\\}$ for some $\epsilon\in[0,1]$. Then with probability at
least $1-\delta$ over $S$, the weighted average output $\bar{w}_{T}$ of M-SPP-
SWoR (Algorithm 3) satisfies

$\displaystyle\mathbb{E}_{\xi_{[T]}}\left[D(\bar{w}_{T},W^{*})\right]$
$\displaystyle\lesssim$
$\displaystyle\frac{\sqrt{LM\log(1/\delta)}\log(T)}{\lambda\rho\sqrt{nT}}+\sqrt{\frac{\rho\left[R(w_{0})-R^{*}\right]}{\lambda
T^{2}}+\frac{L}{\lambda^{2}\rho nT}R^{*}+\frac{\sqrt{\epsilon}}{\lambda
T^{2}}\left(\frac{L}{\lambda\rho}+G\sqrt{\frac{1}{\lambda\rho}}\right)}.$

###### Proof.

See Appendix C.2 for a proof of this result. ∎

######  Remark 11.

We comment on the optimality of the bound in Theorem 6. Consider
$\rho=\mathcal{O}(1)$. The first term of scale
$\mathcal{O}\big{(}\frac{\sqrt{\log(1/\delta)}\log(T)}{\sqrt{nT}}\big{)}$
represents the overhead of getting generalization with high probability over
data. The second term matches the corresponding in-expectation estimation
error bound in Theorem 4, which matches the known optimal rates for strongly
convex SGD (Rakhlin et al., 2012; Dieuleveut et al., 2017). In view of the
minimax lower bounds for statistical estimation (Tsybakov, 2008), the
estimation error bound established in Theorem 6 is near-optimal for strongly
convex risk minimization.

Finally, we provide a high-probability generalization bound of M-SPP for
arbitrary convex population risk functions.

######  Theorem 7.

Suppose that Assumptions 1 and 3 hold and the loss function $\ell$ is bounded
in the interval $[0,M]$. Set $\gamma_{t}\equiv\sqrt{\frac{T}{n}}$. Assume that
$\epsilon_{t}\leq\frac{LM}{4nT^{2}\sqrt{nT}}$. Then with probability at least
$1-\delta$ over $S$, the average output
$\bar{w}_{T}=\frac{1}{T}\sum_{t=1}^{T}w_{t}$ of M-SPP (Algorithm 1) satisfies

$\left|R(\bar{w}_{T})-R_{S}(\bar{w}_{T})\right|\lesssim\frac{(LM+G\sqrt{LM})\log(N)\log(1/\delta)}{\sqrt{nT}}+M\sqrt{\frac{\log\left(1/\delta\right)}{nT}}.$

###### Proof.

See Appendix C.3 for a proof of this result. ∎

We remark in passing that using similar uniform stability argument, the high-
probability generalization bound in Theorem 7 can be shown to hold for convex
and non-smooth loss functions as well. We omit the detailed analysis as it is
out of the scope of this paper focusing on smooth losses.

##  5 Comparison with Prior Methods

_Comparison with M-SPP and SPP methods._ The M-SPP algorithm considered in
this article is a minibatch extension of the SPP methods. The convergence
analysis of SPP has received recent wide attention in stochastic optimization
community. Specially for finite-sum optimization over $N$ data points, an
incremental SPP method was proposed and analyzed in (Bertsekas, 2011). For
learning with linear prediction models and strongly convex Lipschitz-loss,
(Toulis et al., 2016) established a set of $\mathcal{O}(\frac{1}{N^{\gamma}})$
rates of convergence for SPP with suitable $\gamma\in(0.5,1]$, where $N$ is
the iteration counter. For arbitrary convex loss functions, the non-asymptotic
convergence performance of SPP was studied with
$\mathcal{O}(\frac{1}{\sqrt{N}})$ rate obtained for Lipschitz losses (Patrascu
and Necoara, 2017; Davis and Drusvyatskiy, 2019), $\mathcal{O}(\frac{1}{N})$
for strongly convex and Lipschitz (Davis and Drusvyatskiy, 2019) or smooth
(Patrascu and Necoara, 2017) losses, or
$\mathcal{O}\left(\frac{\log(N)}{N}\right)$ rate for strongly convex non-
smooth losses (Asi and Duchi, 2019b). Recently, it has been shown that the
$\mathcal{O}\left(\frac{\log(N)}{N}\right)$ rate also extends to M-SPP with
strongly convex losses (Asi et al., 2020). The asymptotic and non-asymptotic
behaviors of SPP for weakly convex losses (e.g., composite of convex loss with
smooth map) have been studied for stochastic optimization with (Duchi and
Ruan, 2018) or without (Davis and Drusvyatskiy, 2019) composite structures.
Among others, our work is most closely related to the minibatch proximal
update method developed for communication-efficient distributed optimization
(Wang et al., 2017b). Similarly from the viewpoint of algorithmic stability,
the $\mathcal{O}(\frac{1}{N^{\gamma}})$ rates were established for that method
for Lipschitz-loss with arbitrary convexity ($\gamma=0.5$) or strong convexity
($\gamma=1$). In comparison to these prior results, our convergence results
for M-SPP are new in the following aspects:

  * •

The convergence rates are derived for smooth losses and they explicitly show
the impact of noise level of a statistical model, as encoded in $R^{*}$, to
convergence performance which has not been previously known for SPP-type
methods.

  * •

The $\mathcal{O}(N^{-1})$ fast rate attained in this article is valid for
population risks with quadratic growth property, without requiring each
instantaneous loss to be strongly convex.

  * •

We provide a near-optimal model estimation error bound of a sampling-without-
replacement variant of M-SPP that holds with high probability over the
randomness of data while in expectation over the randomness of sampling.

_Comparison with SGD and ERM._ Similar to those in Theorem 1 and Theorem 3,
the bias-variance composite rates have been known for accelerated SGD for
least squares regression (Dieuleveut et al., 2017), or minibatch SGD (M-SGD)
for generic convex and smooth learning problems (Woodworth and Srebro, 2021).
While the results are of similar flavor, we came to the path in a distinct
algorithmic framework using quite different proof techniques. Particularly, in
contrast to Woodworth and Srebro (2021), our analysis neither uses the
knowledge of model scale which is typically inaccessible in real problems, nor
relies on the restarting arguments for strongly convex problems. Also for SGD
with smooth loss functions, a fast rate of $\mathcal{O}(\frac{1}{N})$ has
recently been established via stability theory in the ideally clean case where
the optimal population risk is zero (Lei and Ying, 2020, Theorem 4). With
$\gamma=\mathcal{O}(\frac{1}{n})$, the first bound of our Theorem 3 matches
that bound in the context of M-SPP. For strongly convex problems, our results
in Theorem 1 are stronger than (Lei and Ying, 2020, Theorem 12) in the sense
that the formers (ours) only require the population risk to have quadratic
growth property while the latter requires the loss to be Lipschitz and the
empirical risk to be strongly convex. Finally, for convex ERM, similar
composite risk bounds have been established by Srebro et al. (2010); Zhang et
al. (2017) under somewhat more stringent conditions such as bounded domain of
interest and huge sample with $N\gg p$.

Table 2 summaries a comparison of the risk bounds obtained in this work to
several prior ones for (M-)SPP, (M-)SGD and ERM.

Method | Literature | Risk Bound | Conditions
---|---|---|---
Loss | $R$ | $R_{S}$
M-SPP | Asi et al. (2020) | $\mathcal{O}\left(\frac{\log(N)}{N}\right)$ | s.cvx | — | —
Wang et al. (2017b) | $\mathcal{O}\left(\frac{1}{N}\right)$ | Lip & s.cvx | — | —
Theorem 1 (our work) |  |  $\mathcal{O}\left(\frac{1}{T^{2}}+\frac{R^{*}}{N}\right)$ or
---
$\mathcal{O}\left(\frac{1}{T^{2}}+\frac{1+R^{*}}{N}\right)$
sm & cvx | qg | —
Theorem 3 (our work) |  |  $\mathcal{O}\left(\frac{1}{N}+R^{*}\right)$ or
---
$\mathcal{O}\left(\frac{1+R^{*}}{\sqrt{N}}\right)$
sm & cvx | — | —
SPP | Asi and Duchi (2019b) | $\mathcal{O}\left(\frac{\log(N)}{N}\right)$ | s.cvx | — | —
Patrascu and Necoara (2017) | $\mathcal{O}\left(\frac{1}{N}\right)$ | sm & s.cvx | — |
Davis and Drusvyatskiy (2019) | $\mathcal{O}\left(\frac{1}{N^{2}}+\frac{1}{N}\right)$ | Lip & s.cvx | — |
M-SGD | Woodworth and Srebro (2021) | $\mathcal{O}\left(\frac{1}{T^{2}}+\frac{1}{N}+\sqrt{\frac{R^{*}}{N}}\right)$ | sm & cvx | — | —
$\mathcal{O}\left(e^{-T}+\frac{R^{*}}{N}\right)$ | sm & cvx | qg | —
| Dieuleveut et al. (2017) | $\mathcal{O}\left(\frac{1}{N^{2}}+\frac{R^{*}}{N}\right)$ | quadratic | s.cvx | —
SGD | Lei and Ying (2020) |  |  $\mathcal{O}\left(\frac{1}{N}+R^{*}\right)$ or
---
$\mathcal{O}\left(\frac{1+R^{*}}{\sqrt{N}}\right)$
sm & cvx | — | s.cvx
Rakhlin et al. (2012) | $\mathcal{O}\left(\frac{1}{N}\right)$ |  | Lip &
---
sm & cvx
s.cvx | —
ERM | Zhang et al. (2017) |  |  $\mathcal{O}\left(\frac{p}{N}+\frac{R^{*}}{N}\right)$ or
---
$\mathcal{O}\left(\frac{1}{N^{2}}+\frac{R^{*}}{N}\right)$
for $N\gtrsim p$
sm & cvx |  | Lip
---
& s.cvx
—
Srebro et al. (2010) | $\mathcal{O}\left(\frac{1}{N}+\sqrt{\frac{R^{*}}{N}}\right)$ | sm & cvx | — | —
Table 2: Comparison of our risk bounds to some prior results for M-SPP and SPP
as well as for SGD and ERM. Recall that $T$ is the iteration count and $N$ is
the total number of samples accessed. All the listed bounds hold in
expectation. Here we have used the following abbreviations: cvx (convex),
s.cvx (strongly convex), Lip (Lipschitz continuous), sm (smooth), qg
(quadratic growth).

##  6 Experiments

We carry out a set of numerical study to demonstrate the convergence
performance of minibatch stochastic proximal point methods in (composite)
statistical learning problems, to answer the following 3 questions associated
with the key theory and algorithms established in this article:

  * •

Question 1: _How the size of minibatch and noise level of a statistical
learning model affect the convergence speed of M-SPP for smooth loss
function?_ This question is mainly about verifying Theorem 1 and Theorem 5,
and it is answered through a simulation study on Lasso estimation in Section
6.1.

  * •

Question 2: _Can the two-phase variant of M-SPP improve over M-SPP in the
small- $T$-large-$n$ setting?_ The simulation results presented in Section 6.1
also answer this question related to the verification of Theorem 2.

  * •

Question 3: _How M-SPP(-TP) methods compare with M-SGD in convergence
performance?_ The real-data experimental results on logistic regression tasks
in Section 6.2 answer this question about algorithm comparison.

###  6.1 Simulation Study

We first provide a simulation study to verify our theoretical results for
smooth losses when substantialize to the widely used Lasso regression model
(Wainwright, 2009) with quadratic loss function
$\ell(f_{w}(x),y)=\frac{1}{2}(y-w^{\top}x)^{2}$ and $r(f_{w})=\mu\|w\|_{1}$
where $\mu$ is the $\ell_{1}$-penalty modulus. Given a model parameter
$\bar{w}\in\mathbb{R}^{p}$ and a feature point $x\in\mathbb{R}^{p}$ drawn from
standard Gaussian distribution $\mathcal{N}(0,I_{p\times p})$, the responses
$y$ is generated according to a linear model $y=\bar{w}^{\top}x+\varepsilon$
with a random Gaussian noise $\varepsilon\sim\mathcal{N}(0,\sigma^{2})$. In
this case, the population risk function can be expressed in a close form as

$\displaystyle
R(w)=\frac{1}{2}\|w-\bar{w}\|^{2}+\frac{\sigma^{2}}{2}+\mu\|w\|_{1}.$

Given a set of $T$ random $n$-minibatches
$\left\\{S_{t}=\\{x_{i,t},y_{i,t}\\}_{i\in[n]}\right\\}_{t\in[T]}$ drawn from
the above data distribution, we aim at evaluating the convergence performance
of M-SPP towards the minimizer of $R$ which can be expressed as

$\displaystyle w^{*}=(\bar{w}-\mu)_{+}-(-\bar{w}-\mu)_{+},$

where $(\cdot)_{+}$ is an element-wise function that preserves the positive
parts of a vector.

We test with $p=5000$ and $N=nT=100p$, and consider a well-specified sparse
regression model where the true parameter vector $\bar{w}$ is $\bar{k}$-sparse
with $\bar{k}=0.2p$ and its non-zero entries are sampled from a zero-mean
Gaussian distribution. We set $\mu=10^{-3}$ and initialize $w^{(0)}=0$. The
inner-loop minibatch proximal Lasso subproblems are optimized via a standard
proximal gradient descent method, using either of the following two
termination criteria: 1) the difference between consecutive objective values
is below $10^{-3}$ and 2) the iteration step reaches $1000$.

(a) Results under varying $T$.  (b) Results under varying $\sigma$.  (c) M-SPP
_versus_ M-SPP-TP

Figure 1: Simulation study on Lasso regression: Convergence performances of
M-SPP and M-SPP-TP. The y-axis represents the logarithmic scale of excess
risk.

The following two experimental setups are considered for theory verification:

  * •

We fix the noise level $\sigma=0.1$ and study the impact of varying
$T\in\\{10,20,100,500\\}$ on the convergence of M-SPP. Figure 1(a) shows the
evolving curves of excess risk as functions of sample size, in a semi-log
layout with y-axis representing the logarithmic scale of excess risk. From
this set of curves we can observe a clear trend that in the early stage, M-SPP
converges faster when the total number of minibatches is relatively large
(say, $T\in\\{20,100\\}$). This is consistent with the prediction of Theorem 1
about the impact of $T$ and $n$ on convergence rates. While in the final
stage, relatively slower convergence behavior is exhibited under relatively
larger $T$ (say, $T\in\\{100,500\\}$). This observation can be explained by
the inexact analysis in Theorem 4 which shows that to guarantee the desired
convergence rate, the inner-loop proximal ERM update needs to be extremely
accurate when $T$ is relatively large. Therefore, the question raised in
Question 1 on the impact of minibatch size on convergence rate is answered by
this group of results.

Also in this setup, we have compared M-SPP and its two-phase variant M-SPP-TP
for $T\in\\{5,10\\}$. The related results are shown in Figure 1(c), which
indicate that M-SPP-TP significantly improves the convergence of M-SPP in the
small-$T$-large-$n$ cases. This observation supports the result of Theorem 2
and answers Question 2 affirmatively.

  * •

We fix $T=50$ and study the impact of varying noise level
$\sigma\in\\{0.1,1,5\\}$ on the convergence performance. The results are shown
in Figure 1(b). From this group of results we can see that faster convergence
speed is attained at relatively smaller noise level $\sigma$, while the speed
becomes insensitive to noise level when $\sigma$ is sufficiently small (say,
$\sigma\leq 1$). This is consistent with the predication by Theorem 1, keeping
in mind the fact that
$R^{*}=\frac{1}{2}\|w^{*}-\bar{w}\|^{2}+\frac{\sigma^{2}}{2}+\mu\|w^{*}\|_{1}\leq\|\bar{w}\|^{2}+\frac{1}{2}\sigma^{2}$.
The question raised in Question 1 on the impact of noise level on convergence
performance is answered by this group of results.

###  6.2 Experiment on Real Data

We further compare our methods with M-SGD for binary prediction problems using
the logistic loss $\ell(w^{\top}x,y)=\log(1+\exp(-yw^{\top}x))$. Here the
M-SGD method is implemented by an SGD solver from SGDLibrary (Kasai, 2017).
For M-SPP and M-SPP-TP, the The inner-loop minibatch proximal ERMs are solved
by the same SGD solver applied with a fixed SGD-batch-size $10$ and a single
epoch of data processing. We initialize $w^{(0)}=0$ for all the considered
methods.

We use two public data sets for evaluation: the gisette data (Guyon et al.,
2004) with $p=5000,N=6000$ and the covtype.binary data (Collobert et al.,
2001) with $p=54,N=581,012$ 111Both data sets are available at
https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/.. For each data set,
we use half of the samples as training set and the rest as test set. We are
interested in the impact of minibatch-size $n$ on the prediction performance
of model measured by test error. All the considered stochastic algorithms are
executed with 10 epochs of data processing, and thus the overall number of
minibatches is $T=N/n\times 10$. We replicate each experiment $10$ times over
random split of data and report the results in mean-value along with error
bar.

(a) $n=N/5$ (b) $n=N/20$ (c) $n=N/100$

Figure 2: Real-data results on logistic regression: Test error convergence
comparison on gisette under varying minibatch size.

In Figure 2, we show the evolving curves (error bar shaded in color) of test
error with respect to the number of minibatches accessed on gisette, under
varying minibatch size $n\in\\{\frac{N}{5},\frac{N}{20},\frac{N}{100}\\}$.
From this set of curves we can observe that:

  * •

Under the same minibatch size, M-SPP and M-SPP-TP converge faster and stabler
than M-SGD, especially when the minibatch size is relatively large (see Figure
2(a)). This is as expected because when minibatch size becomes large, M-SGD
approaches to gradient descent method while M-SPP approaches ERMs. This
answers Question 3 raised at the beginning of the experiment section.

  * •

M-SPP-TP exhibits sharper convergence behavior than M-SPP at the early stage
of iteration, especially when the minibatch-size is relatively large. This is
consistent with our theoretical results in Theorem 1 and Theorem 2.

Figure 3 shows the corresponding results on covtype under
$n\in\left\\{\frac{N}{20},\frac{N}{100},\frac{N}{1000}\right\\}$. From this
set of results we once again see that M-SPP and M-SPP-TP consistently
outperform M-SGD under the same minibatch size, and M-SPP-TP converges faster
than M-SPP under relatively larger minibatch size (say, $n=\frac{N}{20}$).

(a) $n=N/20$ (b) $n=N/100$ (c) $n=N/1000$

Figure 3: Real-data results on logistic regression: Test error convergence
comparison on covtype.binary under varying minibatch size.

##  7 Conclusions and Future Prospects

In this article, we presented an improved convergence analysis for the
minibatch stochastic proximal point methods with smooth and convex losses.
Under the quadratic growth condition on population risk, we showed that M-SPP
with minibatch-size $n$ and iteration count $T$ converges at a composite rate
consisting of an $\mathcal{O}(\frac{1}{T^{2}})$ _bias_ decaying component and
an $\mathcal{O}(\frac{1}{N})$ _variance_ decaying component. In the
small-$n$-large-$T$ case, this result substantially improves the prior
relevant results of SPP-type approaches which typically require each
instantaneous loss to be Lipschitz and strongly convex. Complementally in the
small-$T$-large-$n$ setting, we provide a two-phase acceleration of M-SPP
which improves the $\mathcal{O}(\frac{1}{T^{2}})$ bias decaying rate to
$\mathcal{O}\left(\frac{\log(N)}{N^{2}}\right)$. Perhaps the most interesting
theoretical finding is that the (dominant) variance decaying term has a factor
dependence on the minimal value of population risk, justifying the sharper
convergence behavior of M-SPP in low-noise statistical setting as backed up by
our numerical evidence. In addition to the in-expectation risk bounds, we have
also derived a near-optimal parameter estimation error bound for a random
shuffling variant of M-SPP that holds with high probability over data
distribution and in expectation over the random shuffling. To conclude, our
theory lays a novel and stronger foundation for understanding the convex M-SPP
style algorithms that have gained recent significant attention, both in theory
and practice, for large-scale machine learning (Li et al., 2014; Wang et al.,
2017a; Asi et al., 2020).

There are several key prospects for future investigation of our theory:

  * •

It is still open to derive near-optimal exponential excess risk bounds for
M-SPP that apply to the (suffix) average or last of iterates over training
data.

  * •

Inspired by the recent progresses made towards understanding M-SPP with
momentum acceleration (Deng and Gao, 2021; Chadha et al., 2022), it is
interesting to provide momentum and weakly-convex extensions of our theory for
smooth loss functions.

  * •

Last but not least, we expect that the theory developed in this article can be
extended to the setup of non-parametric learning with minibatch stochastic
proximal point methods.

## Acknowledgements

The authors sincerely thank the anonymous referees for their constructive
comments. The work of Xiao-Tong Yuan is also funded in part by the National
Key Research and Development Program of China under Grant No. 2018AAA0100400
and in part by the Natural Science Foundation of China (NSFC) under Grant
No.U21B2049, No.61876090 and No.61936005.

##  Appendix A Proofs for the Results in Section 2

In this section, we present the technical proofs for the main results stated
in Section 2.

###  A.1 Proof of Theorem 1

Here we prove Theorem 1 as restated below for convenience. See 1

We first present the following lemma which will be used in the proof. It can
be viewed as a straightforward extension of the prior result (Wang et al.,
2017b, Lemma 1) to the setup of composite minimization. A proof is included
here for the sake of completeness.

######  Lemma 1.

Assume that the loss function $\ell$ is convex with respect to its first
argument and the regularization function $r$ is convex. Then for any
$w\in\mathcal{W}$, we have

$R_{S_{t}}(w_{t})-R_{S_{t}}(w)\leq\frac{\gamma_{t}}{2}\left(\|w-w_{t-1}\|^{2}-\|w-w_{t}\|^{2}-\|w_{t}-w_{t-1}\|^{2}\right).$

###### Proof.

Since $\ell$ and $r$ are both convex, $R_{S_{t}}$ is convex over
$\mathcal{W}$. The optimality of $w_{t}$ implies that for any
$w\in\mathcal{W}$ and $\eta\in(0,1)$

$\displaystyle R_{S_{t}}(w_{t})+\frac{\gamma_{t}}{2}\|w_{t}-w_{t-1}\|^{2}\leq
R_{S_{t}}((1-\eta)w_{t}+\eta w)+\frac{\gamma_{t}}{2}\|(1-\eta)w_{t}+\eta
w-w_{t-1}\|^{2}$ $\displaystyle\leq$
$\displaystyle(1-\eta)R_{S_{t}}(w_{t})+\eta
R_{S_{t}}(w)+\frac{\gamma_{t}}{2}\left[(1-\eta)\|w_{t}-w_{t-1}\|^{2}+\eta\|w-w_{t-1}\|^{2}-\eta(1-\eta)\|w-w_{t}\|^{2}\right],$

where in the last inequality we have used the definition of the norm
$\|\cdot\|$. Rearranging both sides of the above inequality yields

$\eta(R_{S_{t}}(w_{t})-R_{S_{t}}(w))\leq\frac{\eta\gamma_{t}}{2}\left[\|w-w_{t-1}\|^{2}-(1-\eta)\|w-w_{t}\|^{2}-\|w_{t}-w_{t-1}\|^{2}\right],$

which then implies (keep in mind that $\eta>0$)

$R_{S_{t}}(w_{t})-R_{S_{t}}(w)\leq\frac{\gamma_{t}}{2}\left[\|w-w_{t-1}\|^{2}-(1-\eta)\|w-w_{t}\|^{2}-\|w_{t}-w_{t-1}\|^{2}\right].$

Limiting $\eta\rightarrow 0^{+}$ in the above inequality yields the desired
bound. ∎

The following boundedness result for smooth function is due to Srebro et al.
(2010, Lemma 3.1).

######  Lemma 2.

If $g$ is non-negative and $L$-smooth, then $\|\nabla
g(w)\|\leq\sqrt{2Lg(w)}$.

Let $\\{\mathcal{F}_{t}\\}_{t\geq 1}$ be the filtration generated by the
iterates $\\{w_{t}\\}_{t\geq 1}$ as
$\mathcal{F}_{t}=\sigma\left(w_{1},w_{2},...,w_{t}\right)$. With Lemma 1 and
Lemma 2 in place, we can further establish the following key lemma that plays
a fundamental role in proving Theorem 1.

######  Lemma 3.

Suppose that the Assumptions 1 holds. Set $\gamma_{t}\geq\frac{16L}{n}$. Then
we have

$\mathbb{E}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]\leq\gamma_{t}\left(D^{2}(w_{t-1},W^{*})-\mathbb{E}\left[D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]\right)+\frac{16L}{\gamma_{t}n}R^{*}.$

###### Proof.

Let us consider a sample set $S_{t}^{(i)}$ which is identical to $S_{t}$
except that one of the $z_{i,t}$ is replaced by another random sample
$z^{\prime}_{i,t}$. Denote

$w^{(i)}_{t}=\mathop{\arg\min}_{w\in\mathcal{W}}\left\\{F^{(i)}_{t}(w):=R_{S^{(i)}_{t}}(w)+\frac{\gamma_{t}}{2}\|w-w_{t-1}\|^{2}\right\\},$

where $R_{S^{(i)}_{t}}(w):=\frac{1}{n}\left(\sum_{j\neq
i}\ell(w;z_{j,t})+\ell(w;z^{\prime}_{i,t})\right)+r(w)$. Then we can show that

$\displaystyle F_{t}(w^{(i)}_{t})-F_{t}(w_{t})$ $\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{j\neq
i}\left(\ell(w^{(i)}_{t};z_{j,t})-\ell(w_{t};z_{j,t})\right)+\frac{1}{n}\left(\ell(w^{(i)}_{t};z_{i,t})-\ell(w_{t};z_{i,t})\right)$
$\displaystyle+r(w^{(i)}_{t})-r(w_{t})+\frac{\gamma_{t}}{2}\|w^{(i)}_{t}-w_{t-1}\|^{2}-\frac{\gamma_{t}}{2}\|w_{t}-w_{t-1}\|^{2}$
$\displaystyle=$ $\displaystyle
F^{(i)}_{t}(w^{(i)}_{t})-F^{(i)}_{t}(w_{t})+\frac{1}{n}\left(\ell(w^{(i)}_{t};z_{i,t})-\ell(w_{t};z_{i,t})\right)-\frac{1}{n}\left(\ell(w^{(i)}_{t};z^{\prime}_{i,t})-\ell(w_{t};z^{\prime}_{i,t})\right)$
$\displaystyle\leq$
$\displaystyle\frac{1}{n}\left|\ell(w^{(i)}_{t};z_{i,t})-\ell(w_{t};z_{i,t})\right|+\frac{1}{n}\left|\ell(w^{(i)}_{t};z^{\prime}_{i,t})-\ell(w_{t};z^{\prime}_{i,t})\right|$
$\displaystyle\overset{\zeta_{1}}{\leq}$
$\displaystyle\frac{\|\nabla\ell(w^{(i)}_{t};z_{i,t})\|+\|\nabla\ell(w_{t};z^{\prime}_{i,t})\|}{n}\|w^{(i)}_{t}-w_{t}\|$
$\displaystyle\overset{\zeta_{2}}{\leq}$
$\displaystyle\frac{\sqrt{2L\ell(w^{(i)}_{t};z_{i,t})}+\sqrt{2L\ell(w_{t};z^{\prime}_{i,t})}}{n}\|w^{(i)}_{t}-w_{t}\|,$

where “$\zeta_{1}$” is due to the convexity of loss and in “$\zeta_{2}$”we
have used Lemma 2. The bound in Lemma 1 implies

$F_{t}(w^{(i)}_{t})-F_{t}(w_{t})\geq\frac{\gamma_{t}}{2}\|w^{(i)}_{t}-w_{t}\|^{2}.$

Combining the preceding two inequalities yields

$\frac{\gamma_{t}}{2}\|w^{(i)}_{t}-w_{t}\|\leq\frac{\sqrt{2L\ell(w^{(i)}_{t};z_{i,t})}+\sqrt{2L\ell(w_{t};z^{\prime}_{i,t})}}{n},$

which immediately gives

$\|w^{(i)}_{t}-w_{t}\|\leq\frac{2\left(\sqrt{2L\ell(w^{(i)}_{t};z_{i,t})}+\sqrt{2L\ell(w_{t};z^{\prime}_{i,t})}\right)}{\gamma_{t}n}.$
(8)

Let us now consider the following population risk and empirical risk over
$S_{t}$ with respect to the loss function $\ell$:

$R^{\ell}(w):=\mathbb{E}_{(x,y)\sim\mathcal{D}}[\ell(w;z)],\quad
R^{\ell}_{S_{t}}(w):=\frac{1}{n}\sum_{i=1}^{n}\ell(w;z_{i,t}).$

Since $S_{t}$ and $S_{t}^{(i)}$ are both i.i.d. samples of the data
distribution. It follows that

$\displaystyle\mathbb{E}_{S_{t}}\left[R^{\ell}(w_{t})\mid\mathcal{F}_{t-1}\right]=\mathbb{E}_{S_{t}\cup\\{z^{\prime}_{i,t}\\}}\left[\ell(w_{t};z^{\prime}_{i,t})\mid\mathcal{F}_{t-1}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{S_{t}^{(i)}}\left[R^{\ell}(w^{(i)}_{t})\mid\mathcal{F}_{t-1}\right]=\mathbb{E}_{S_{t}^{(i)}\cup\\{z_{i,t}\\}}\left[\ell(w^{(i)}_{t};z_{i,t})\mid\mathcal{F}_{t-1}\right]$
.

Since the above holds for all $i=1,...,n$, we can further show that

$\displaystyle\mathbb{E}_{S_{t}}\left[R^{\ell}(w_{t})\mid\mathcal{F}_{t-1}\right]$
(9) $\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{S_{t}^{(i)}\cup\\{z_{i,t}\\}}\left[\ell(w^{(i)}_{t};z_{i,t})\mid\mathcal{F}_{t-1}\right]=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{S_{t}\cup\\{z^{\prime}_{i,t}\\}}\left[\ell(w^{(i)}_{t};z_{i,t})\mid\mathcal{F}_{t-1}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{S_{t}\cup\\{z^{\prime}_{i,t}\\}}\left[\ell(w_{t};z^{\prime}_{i,t})\mid\mathcal{F}_{t-1}\right]=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{S^{(i)}_{t}\cup\\{z_{i,t}\\}}\left[\ell(w_{t};z^{\prime}_{i,t})\mid\mathcal{F}_{t-1}\right].$

Regarding the empirical case, we find that

$\displaystyle\mathbb{E}_{S_{t}}\left[R^{\ell}_{S_{t}}(w_{t})\mid\mathcal{F}_{t-1}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{S_{t}}\left[\ell(w_{t};z_{i,t})\mid\mathcal{F}_{t-1}\right]=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{S_{t}\cup\\{z^{\prime}_{i,t}\\}}\left[\ell(w_{t};z_{i,t})\mid\mathcal{F}_{t-1}\right].$

Combining the preceding two equalities gives that

$\displaystyle\left|\mathbb{E}_{S_{t}}\left[R(w_{t})-R_{S_{t}}(w_{t})\mid\mathcal{F}_{t-1}\right]\right|$
$\displaystyle=$
$\displaystyle\left|\mathbb{E}_{S_{t}}\left[R^{\ell}(w_{t})-R^{\ell}_{S_{t}}(w_{t})\mid\mathcal{F}_{t-1}\right]\right|$
$\displaystyle=$
$\displaystyle\left|\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{S_{t}\cup\\{z^{\prime}_{i,t}\\}}\left[\ell(w^{(i)}_{t};z_{i,t})-\ell(w_{t};z_{i,t})\mid\mathcal{F}_{t-1}\right]\right|$
$\displaystyle\leq$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{S_{t}\cup\\{z^{\prime}_{i,t}\\}}\left[\left|\ell(w^{(i)}_{t};z_{i,t})-\ell(w_{t};z_{i,t})\right|\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{S_{t}\cup\\{z^{\prime}_{i,t}\\}}\left[\sqrt{2L\ell(w^{(i)}_{t};z_{i,t})}\|w^{(i)}_{t}-w_{t}\|\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\overset{\eqref{inequat:proof_exp_risk_smooth_key_1}}{\leq}$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{S^{(i)}_{t}\cup\\{z_{i,t}\\}}\left[\frac{4L\ell(w^{(i)}_{t};z_{i,t})}{\gamma_{t}n}+\frac{4L\sqrt{\ell(w^{(i)}_{t};z_{i,t})\ell(w_{t};z^{\prime}_{i,t})}}{\gamma_{t}n}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\overset{\zeta_{1}}{\leq}$
$\displaystyle\left(\frac{L}{\gamma_{t}n}\right)\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{S_{t}\cup\\{z^{\prime}_{i,t}\\}}\left[6\ell(w^{(i)}_{t};z_{i,t})+2\ell(w_{t};z^{\prime}_{i,t})\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\overset{\eqref{inequat:proof_exp_risk_smooth_key_2}}{=}$
$\displaystyle\frac{8L}{\gamma_{t}n}\mathbb{E}_{S_{t}}\left[R^{\ell}(w_{t})\mid\mathcal{F}_{t-1}\right]\leq\frac{8L}{\gamma_{t}n}\mathbb{E}_{S_{t}}\left[R(w_{t})\mid\mathcal{F}_{t-1}\right],$

where in “$\zeta_{1}$” we have used the fact $a^{2}+b^{2}\geq 2ab$ and the
last inequality is due to the fact $r\geq 0$.

Let us now denote $w^{*}_{t}=\mathop{\arg\min}_{w\in W^{*}}\|w-w_{t}\|$.
Conditioned on $\mathcal{F}_{t-1}$, taking expectation on both sides of the
bound in Lemma 1 for $w=w^{*}_{t-1}$ yields

$\displaystyle\mathbb{E}_{S_{t}}\left[R_{S_{t}}(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$
$\displaystyle\frac{\gamma_{t}}{2}\mathbb{E}_{S_{t}}\left[\|w^{*}_{t-1}-w_{t-1}\|^{2}-\|w^{*}_{t-1}-w_{t}\|^{2}-\|w_{t}-w_{t-1}\|^{2}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$
$\displaystyle\frac{\gamma_{t}}{2}\left(\|w^{*}_{t-1}-w_{t-1}\|^{2}-\mathbb{E}_{S_{t}}\left[\|w^{*}_{t}-w_{t}\|^{2}\mid\mathcal{F}_{t-1}\right]\right).$

Combining the preceding two inequalities yields

$\displaystyle\mathbb{E}_{S_{t}}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{S_{t}}\left[R(w_{t})-R_{S_{t}}(w_{t})+R_{S_{t}}(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$
$\displaystyle\left|\mathbb{E}_{S_{t}}\left[R(w_{t})-R_{S_{t}}(w_{t})\mid\mathcal{F}_{t-1}\right]\right|+\mathbb{E}_{S_{t}}\left[R_{S_{t}}(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$
$\displaystyle\frac{\gamma_{t}}{2}\left(\|w^{*}_{t-1}-w_{t-1}\|^{2}-\mathbb{E}_{S_{t}}\left[\|w^{*}_{t}-w_{t}\|^{2}\mid\mathcal{F}_{t-1}\right]\right)+\frac{8L}{\gamma_{t}n}\mathbb{E}_{S_{t}}\left[R(w_{t})\mid\mathcal{F}_{t-1}\right]$
$\displaystyle=$
$\displaystyle\frac{\gamma_{t}}{2}\left(\|w_{t-1}^{*}-w_{t-1}\|^{2}-\mathbb{E}_{S_{t}}\left[\|w^{*}_{t}-w_{t}\|^{2}\mid\mathcal{F}_{t-1}\right]\right)+\frac{8L}{\gamma_{t}n}\mathbb{E}_{S_{t}}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]+\frac{8L}{\gamma_{t}n}\mathbb{E}_{S_{t}}\left[R^{*}\right]$
$\displaystyle\leq$
$\displaystyle\frac{\gamma_{t}}{2}\left(\|w^{*}_{t-1}-w_{t-1}\|^{2}-\mathbb{E}_{S_{t}}\left[\|w^{*}_{t}-w_{t}\|^{2}\mid\mathcal{F}_{t-1}\right]\right)+\frac{1}{2}\mathbb{E}_{S_{t}}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]+\frac{8L}{\gamma_{t}n}R^{*},$

where in the last inequality we have used the condition
$\gamma_{t}\geq\frac{52L}{n}$. After rearranging the terms in the above
inequality we obtain

$\displaystyle\mathbb{E}_{S_{t}}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]\leq$
$\displaystyle\gamma_{t}\left(\|w^{*}_{t-1}-w_{t-1}\|^{2}-\mathbb{E}_{S_{t}}\left[\|w^{*}_{t}-w_{t}\|^{2}\mid\mathcal{F}_{t-1}\right]\right)+\frac{16L}{\gamma_{t}n}R^{*}$
$\displaystyle=$
$\displaystyle\gamma_{t}\left(D^{2}(w_{t-1},W^{*})-\mathbb{E}_{S_{t}}\left[D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]\right)+\frac{16L}{\gamma_{t}n}R^{*}.$

This implies the desired bound. ∎

The following lemma is a direct consequence of Lemma 3.

######  Lemma 4.

Suppose that the Assumptions 1 holds. Set $\gamma_{t}\geq\frac{16L}{n}$. Then
the following holds for all $t\geq 1$:

$\mathbb{E}\left[D^{2}(w_{t},W^{*})\right]\leq
D^{2}(w_{0},W^{*})+\sum_{\tau=1}^{t}\frac{16L}{\gamma_{\tau}^{2}n}R^{*}.$

###### Proof.

Since $R(w_{t})\geq R^{*}$ and $\gamma_{t}\geq\frac{52L}{n}$, the bound in
Lemma 3 immediately implies that

$\mathbb{E}_{S_{t}}\left[D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]\leq
D^{2}(w_{t-1},W^{*})+\frac{16L}{\gamma^{2}_{t}n}R^{*}.$ (10)

By unfolding the above recurrent from time instance $t$ to zero we obtain that
for all $t\geq 1$,

$\mathbb{E}\left[D^{2}(w_{t},W^{*})\right]\leq
D^{2}(w_{0},W^{*})+\sum_{\tau=1}^{t}\frac{16L}{\gamma_{\tau}^{2}n}R^{*}.$

This proves the desired bound. ∎

With all these lemmas in place, we are now ready to prove the main result in
Theorem 1.

###### Proof of Theorem 1.

Part (a): Note that the condition on $n$ implies $\gamma_{t}=\frac{\lambda\rho
t}{4}\geq\frac{\lambda\rho}{4}\geq\frac{16L}{n}$. Applying Lemma 3 along with
the condition $R(w_{t})-R^{*}\geq\frac{\lambda}{2}D^{2}(w_{t},W^{*})$ yields

$\displaystyle(1-\rho)\mathbb{E}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$
$\displaystyle\gamma_{t}D^{2}(w_{t-1},W^{*})-\left(\gamma_{t}+\frac{\lambda\rho}{2}\right)\mathbb{E}\left[D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]+\frac{2^{4}L}{\gamma_{t}n}R^{*}$
$\displaystyle\leq$ $\displaystyle\frac{\lambda\rho
t}{4}D^{2}(w_{t-1},W^{*})-\frac{\lambda\rho(t+2)}{4}\mathbb{E}\left[D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]+\frac{2^{6}L}{\lambda\rho
tn}R^{*}$ $\displaystyle\leq$ $\displaystyle\frac{\lambda\rho
t}{4}D^{2}(w_{t-1},W^{*})-\frac{\lambda\rho(t+2)}{4}\mathbb{E}\left[D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]+\frac{2^{7}L}{\lambda\rho(t+1)n}R^{*},$

where in the last inequality we have used $\frac{1}{t}\leq\frac{2}{t+1}$ for
$t\geq 1$. The above inequality implies

$\displaystyle t\mathbb{E}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$
$\displaystyle(t+1)\mathbb{E}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$ $\displaystyle\frac{\lambda\rho
t(t+1)}{4(1-\rho)}D^{2}(w_{t-1},W^{*})-\frac{\lambda\rho(t+1)(t+2)}{4(1-\rho)}\mathbb{E}\left[D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]+\frac{2^{7}L}{\lambda
n\rho(1-\rho)}R^{*}.$

Then based on the law of total expectation and after proper rearrangement we
obtain

$\displaystyle t\mathbb{E}\left[R(w_{t})-R^{*}\right]$ (11)
$\displaystyle\leq$ $\displaystyle\frac{\lambda\rho
t(t+1)}{4(1-\rho)}\mathbb{E}\left[D^{2}(w_{t-1},W^{*})\right]-\frac{\lambda\rho(t+1)(t+2)}{4(1-\rho)}\mathbb{E}\left[D^{2}(w_{t},W^{*})\right]+\frac{2^{7}L}{\lambda
n\rho(1-\rho)}R^{*}.$

By summing the above inequality from $t=1,...,T$ and after normalization we
obtain

$\displaystyle\frac{2}{T(T+1)}\sum_{t=1}^{T}t\mathbb{E}\left[R(w_{t})-R^{*}\right]\leq$
$\displaystyle\frac{\lambda\rho}{T(T+1)(1-\rho)}D^{2}(w_{0},W^{*})+\frac{2^{8}L}{\lambda\rho(1-\rho)(T+1)n}R^{*}$
$\displaystyle\leq$
$\displaystyle\frac{2\lambda\rho}{T(T+1)}D^{2}(w_{0},W^{*})+\frac{2^{9}L}{\lambda\rho(T+1)n}R^{*},$

where in the last inequality we have used $\rho\leq 0.5$. Consider the
weighted output $\bar{w}_{T}=\frac{2}{T(T+1)}\sum_{t=1}^{T}tw_{t}$. In view of
the above inequality and the convexity and quadratic growth property of the
risk function $R$ we have

$\mathbb{E}\left[R(\bar{w}_{T})-R^{*}\right]\leq\frac{4\rho\left[R(w_{0})-R^{*}\right]}{T(T+1)}+\frac{2^{9}L}{\lambda\rho
n(T+1)}R^{*},$

which then implies the desired bound in part (a).

Part (b): Note that $\gamma_{t}=\frac{\lambda\rho
t}{4}+\frac{16L}{n}\geq\frac{16L}{n}$ for all $t\geq 1$. According to Lemma 4
we have the following holds for all $t\geq 1$:

$\displaystyle\mathbb{E}\left[D^{2}(w_{t},W^{*})\right]$ (12)
$\displaystyle\leq$ $\displaystyle
D^{2}(w_{0},W^{*})+\sum_{\tau=1}^{t}\frac{16L}{\gamma_{\tau}^{2}n}R^{*}\leq
D^{2}(w_{0},W^{*})+\frac{2^{8}L}{\lambda^{2}\rho^{2}n}R^{*}\sum_{\tau=1}^{t}\frac{1}{\tau^{2}}\leq
D^{2}(w_{0},W^{*})+\frac{2^{9}L}{\lambda^{2}\rho^{2}n}R^{*}.$

Similar to the argument in part (a), applying Lemma 3 along with the quadratic
growth condition $R(w_{t})-R^{*}\geq\frac{\lambda}{2}D^{2}(w_{t},W^{*})$ and
$\rho\leq 0.5$ yields

$\displaystyle\frac{1}{2}\mathbb{E}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$
$\displaystyle(1-\rho)\mathbb{E}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$
$\displaystyle\gamma_{t}D^{2}(w_{t-1},W^{*})-\left(\gamma_{t}+\frac{\lambda\rho}{2}\right)\mathbb{E}\left[D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]+\frac{2^{4}L}{\gamma_{t}n}R^{*}$
$\displaystyle\leq$ $\displaystyle\frac{\lambda\rho
t}{4}D^{2}(w_{t-1},W^{*})-\frac{\lambda\rho(t+2)}{4}\mathbb{E}\left[D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]$
$\displaystyle+\frac{16L}{n}\left(D^{2}(w_{t-1},W^{*})-\mathbb{E}\left[D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]\right)+\frac{2^{6}L}{\lambda\rho
tn}R^{*},$

where in the second inequality we have used $\gamma_{t}\geq\frac{52L}{n}$, and
in the last inequality we have used $\gamma_{t}\geq\frac{\lambda\rho t}{4}$.
Then based on the law of total expectation and after proper rearrangement we
have

$\displaystyle\mathbb{E}\left[R(w_{t})-R^{*}\right]$ $\displaystyle\leq$
$\displaystyle\frac{\lambda\rho
t}{2}\mathbb{E}\left[D^{2}(w_{t-1},W^{*})\right]-\frac{\lambda\rho(t+2)}{2}\mathbb{E}\left[D^{2}(w_{t},W^{*})\right]$
$\displaystyle+\frac{2^{5}L}{n}\left(\mathbb{E}\left[D^{2}(w_{t-1},W^{*})\right]-\mathbb{E}\left[D^{2}(w_{t},W^{*})\right]\right)+\frac{2^{7}L}{\lambda
tn\rho}R^{*},$

which implies that

$\displaystyle t\mathbb{E}\left[R(w_{t})-R^{*}\right]$ $\displaystyle\leq$
$\displaystyle(t+1)\mathbb{E}\left[R(w_{t})-R^{*}\right]$ $\displaystyle\leq$
$\displaystyle\frac{\lambda\rho
t(t+1)}{2}\mathbb{E}\left[D^{2}(w_{t-1},W^{*})\right]-\frac{\lambda\rho(t+1)(t+2)}{2}\mathbb{E}\left[D^{2}(w_{t},W^{*})\right]$
$\displaystyle+\frac{2^{5}L(t+1)}{n}\left(\mathbb{E}\left[D^{2}(w_{t-1},W^{*})\right]-\mathbb{E}\left[D^{2}(w_{t},W^{*})\right]\right)+\frac{2^{7}L(t+1)}{\lambda
tn\rho}R^{*}$ $\displaystyle\leq$ $\displaystyle\frac{\lambda\rho
t(t+1)}{2}\mathbb{E}\left[D^{2}(w_{t-1},W^{*})\right]-\frac{\lambda\rho(t+1)(t+2)}{2}\mathbb{E}\left[D^{2}(w_{t},W^{*})\right]$
$\displaystyle+\frac{2^{6}Lt}{n}\left(\mathbb{E}\left[D^{2}(w_{t-1},W^{*})\right]-\mathbb{E}\left[D^{2}(w_{t},W^{*})\right]\right)+\frac{2^{8}L}{\lambda
n\rho}R^{*},$

where in the last inequality we have used the fact $t+1\leq 2t$ as $t\geq 1$.
By summing the above inequality from $t=1,...,T$ and after normalization we
obtain

$\displaystyle\frac{2}{T(T+1)}\sum_{t=1}^{T}t\mathbb{E}\left[R(w_{t})-R^{*}\right]$
$\displaystyle\leq$
$\displaystyle\frac{2\lambda\rho}{T(T+1)}D^{2}(w_{0},W^{*})+\frac{2^{7}L}{nT(T+1)}\sum_{t=1}^{T}D^{2}(w_{t-1},W^{*})+\frac{2^{9}L}{\lambda\rho(T+1)n}R^{*}$
$\displaystyle\leq$
$\displaystyle\frac{2\lambda\rho}{T(T+1)}D^{2}(w_{0},W^{*})+\frac{2^{7}L}{nT(T+1)}\sum_{t=1}^{T}\left(D^{2}(w_{0},W^{*})+\frac{2^{9}L}{\lambda^{2}\rho^{2}n}R^{*}\right)+\frac{2^{9}L}{\lambda\rho(T+1)n}R^{*}$
$\displaystyle=$
$\displaystyle\left(\frac{2\lambda\rho}{T(T+1)}+\frac{2^{7}L}{n(T+1)}\right)D^{2}(w_{0},W^{*})+\left(\frac{2^{16}L^{2}}{\lambda^{2}\rho^{2}n^{2}(T+1)}+\frac{2^{9}L}{\lambda\rho
n(T+1)}\right)R^{*},$

where in the last inequality we have used (12). Using the convexity and
quadratic growth property in the above inequality yields

$\mathbb{E}\left[R(\bar{w}_{T})-R^{*}\right]\leq\left(\frac{4\rho}{T(T+1)}+\frac{2^{8}L}{\lambda
n(T+1)}\right)[R(w_{0})-R^{*}]+\left(\frac{2^{16}L^{2}}{\lambda^{2}\rho^{2}n^{2}(T+1)}+\frac{2^{9}L}{\lambda\rho
n(T+1)}\right)R^{*},$

which then implies the desired bound in part (b). The proof is concluded. ∎

###  A.2 Proof of Theorem 2

In this subsection we prove Theorem 2 which is restated below. See 2

###### Proof.

Part (a): In Phase-I, by invoking the first part of Theorem 1 with $\rho=1/2$
and $T=n/m\geq 1$ (with slight abuse of notation) we get immediately that

$\mathbb{E}_{S_{1}}\left[R(w_{1})-R^{*}\right]\leq\frac{2m^{2}\left[R(w_{0})-R^{*}\right]}{n^{2}}+\frac{2^{10}L}{\lambda
n}R^{*}.$ (13)

In Phase-II, conditioned on $\mathcal{F}_{1}$, summing the recursion form (11)
from $t=2,...,T$ with $\rho=1/2$ and proper normalization yields

$\displaystyle\frac{2}{(T-1)(T+2)}\sum_{t=2}^{T}t\mathbb{E}_{S_{2:t}}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{1}\right]$
$\displaystyle\leq$ $\displaystyle\frac{6\lambda
D^{2}(w_{1},W^{*})}{(T-1)(T+2)}+\frac{2^{10}L}{\lambda
n(T+2)}R^{*}\leq\frac{3\left(R(w_{1})-R^{*}\right)}{(T-1)(T+2)}+\frac{2^{10}L}{\lambda
n(T+2)}R^{*},$

where in the last inequality we have used the quadratic growth property.
Consider the weighted average output
$\bar{w}_{T}=\frac{2}{(T-1)(T+2)}\sum_{t=2}^{T}tw_{t}$. Based on the above
inequality and law of total expectation we must have

$\displaystyle\mathbb{E}\left[R(\bar{w}_{T})-R^{*}\right]\leq$
$\displaystyle\frac{6\mathbb{E}_{S_{1}}\left[R(w_{1})-R^{*}\right]}{(T-1)(T+2)}+\frac{2^{10}L}{\lambda
n(T+2)}R^{*}$ $\displaystyle\leq$
$\displaystyle\frac{6\mathbb{E}_{S_{1}}\left[R(w_{1})-R^{*}\right]}{T^{2}}+\frac{2^{102}L}{\lambda
nT}R^{*}$ $\displaystyle\leq$
$\displaystyle\frac{12m^{2}\left[R(w_{0})-R^{*}\right]}{n^{2}T^{2}}+\frac{2^{13}L}{\lambda
nT}R^{*}$ $\displaystyle\leq$
$\displaystyle\frac{2^{22}L^{2}\left[R(w_{0})-R^{*}\right]}{\lambda^{2}n^{2}T^{2}}+\frac{2^{13}L}{\lambda
nT}R^{*},$

where we have used the fact $T\geq 2$ in multiple places and in the last but
one step we have used (13). This immediately implies the desired bound in Part
(a).

Part (b): In Phase-I, by applying second part of Theorem 1 (with $\rho=1/2$
and $T=n/m\geq 1$) and preserving the leading terms we obtain that

$\displaystyle\mathbb{E}_{S_{1}}\left[R(w_{1})-R^{*}\right]\lesssim$
$\displaystyle\left(\frac{m^{2}}{n^{2}}+\frac{L}{\lambda
n}\right)[R(w_{0})-R^{*}]+\left(\frac{L^{2}}{\lambda^{2}mn}+\frac{L}{\lambda
n}\right)R^{*}$ (14) $\displaystyle\lesssim$ $\displaystyle\frac{L}{\lambda
n}[R(w_{0})-R^{*}]+\frac{L^{2}}{\lambda^{2}n}R^{*}.$

In Phase-II, based on the proof argument of the part (b) of Theorem 1 we can
show that the weighted average output
$\bar{w}_{T}=\frac{2}{(T-1)(T+2)}\sum_{t=2}^{T}tw_{t}$ satisfies

$\displaystyle\mathbb{E}\left[R(\bar{w}_{T})-R^{*}\right]\lesssim$
$\displaystyle\left(\frac{1}{T^{2}}+\frac{L}{\lambda
nT}\right)\mathbb{E}_{S_{1}}\left[R(w_{1})-R^{*}\right]+\left(\frac{L^{2}}{\lambda^{2}n^{2}T}+\frac{L}{\lambda
nT}\right)R^{*}$ $\displaystyle\lesssim$ $\displaystyle\left(\frac{L}{\lambda
nT^{2}}+\frac{L^{2}}{\lambda^{2}n^{2}T}\right)[R(w_{0})-R^{*}]+\left(\frac{L^{3}}{\lambda^{3}n^{2}T}+\frac{L^{2}}{\lambda^{2}nT}\right)R^{*}$
$\displaystyle\lesssim$
$\displaystyle\frac{L^{2}}{\lambda^{2}nT}[R(w_{0})-R^{*}]+\frac{L^{3}}{\lambda^{3}nT}R^{*},$

where in the second step we have used (14). This proves the desired bound in
Part (b). ∎

###  A.3 Proof of Theorem 3

In this subsection, we prove Theorem 3 as following restated. See 3

###### Proof.

Since $\gamma_{t}\equiv\gamma\geq\frac{16L}{n}$, the bound in Lemma 3 is
valid. Based on law of total expectation and by summing that inequality from
$t=1,...,T$ with proper normalization we obtain

$\displaystyle\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\left[R(w_{t})-R^{*}\right]\leq\frac{\gamma}{T}D^{2}(w_{0},W^{*})+\frac{16L}{\gamma
n}R^{*}.$

Consider $\bar{w}_{T}=\frac{1}{T}\sum_{t=1}^{T}w_{t}$. In view of the above
inequality and convexity of $R$ we have

$\mathbb{E}\left[R(\bar{w}_{T})-R^{*}\right]\leq\frac{\gamma}{T}D^{2}(w_{0},W^{*})+\frac{16L}{\gamma
n}R^{*}.$

This proves the first desired bound. The second bound follows immediately by
substituting $\gamma=\sqrt{\frac{T}{n}}+\frac{16L}{n}>\frac{16L}{n}$ into the
above bound. The proof is concluded. ∎

###  A.4 On the (Iteration) Stability of M-SPP

In this appendix subsection, we further provide a sensitivity analysis of
M-SPP to the choice of regularization modulus $\\{\gamma_{t}\\}_{t\geq 1}$,
under the following notion of iteration stability essentially introduced by
Asi and Duchi (2019a, b).

######  Definition 2.

A stochastic optimization algorithm generating iterates $\\{w_{t}\\}_{t\geq
1}$ for minimizing the population risk $R(w)$ is staid to be stable if

$\sup_{t\geq 1}D(w_{t},W^{*})<\infty,\ \ \ \text{with probability 1}.$

Before presenting the main results on the iteration stability of M-SPP, we
first recall the Robbins-Siegmund nonnegative almost supermartingale
convergence lemma which is typically used for establishing the stability and
convergence of stochastic optimization methods including SPP (Asi and Duchi,
2019b).

######  Lemma 5 (Robbins and Siegmund (1971)).

Consider four sequences of nonnegative random variables
$\\{U_{t}\\},\\{V_{t}\\},\\{\alpha_{t}\\},\\{\beta_{t}\\}$ that are measurable
over a filtration $\\{\mathcal{F}_{t}\\}_{t\geq 0}$. Suppose that
$\sum_{t}\alpha_{t}<\infty$, $\sum_{t}\beta_{t}<\infty$, and

$\mathbb{E}[U_{t+1}\mid\mathcal{F}_{t}]\leq(1+\alpha_{t})U_{t}+\beta_{t}-V_{t}.$

Then there exits $U_{\infty}$ such that $U_{t}\xrightarrow{a.s.}U_{\infty}$
and $\sum_{t}V_{t}<\infty$ with probability 1.

The following proposition shows that the sequence of estimation error
$\\{\|w_{t}-w^{*}\|\\}$ is non-divergent in expectation and it converges to
some finite value and is bounded with probability $1$.

######  Proposition 2.

Suppose that the Assumptions 1 holds. Assume that
$\gamma_{t}\geq\frac{16L}{n}$ and $\sum_{t\geq 1}L\gamma_{t}^{-2}<\infty$.
Then we have the following hold:

  * (a)

$\mathbb{E}\left[D(w_{t},W^{*})\right]<\infty$;

  * (b)

$D(w_{t},W^{*})$ converges to some finite value and $\sup_{t\geq
1}D(w_{t},W^{*})<\infty$ with probability $1$.

###### Proof.

Applying Lemma 4 yields that for all $t\geq 1$

$\mathbb{E}\left[D^{2}(w_{t},W^{*})\right]\lesssim
D^{2}(w_{0},W^{*})+\sum_{\tau=1}^{t}\frac{L}{\gamma_{\tau}^{2}n}R^{*}<\infty,$

where we have used the given conditions on $\gamma_{t}$. This proves the part
(a). To show the part (b), invoking Lemma 5 with $\alpha_{t}=V_{t}\equiv 0$
and $\beta_{t}=\frac{16L}{\gamma^{2}_{t}n}R^{*}$ to (10) yields
$D(w_{t},W^{*})$ converges to some finite value and thus $\sup_{t\geq
1}D(w_{t},W^{*})<\infty$ almost surely. ∎

######  Remark 12.

Proposition 2 shows that in contrast to minibatch SGD, the choice of
$\gamma_{t}$ in M-SPP is insensitive to the gradient scale of loss functions
for generating a non-divergent sequence of estimation errors.

##  Appendix B Proofs for the Results in Section 3

In this section, we present the technical proofs for the main results stated
in Section 3.

###  B.1 Proof of Theorem 4

In this subsection, we prove Theorem 4 which is restated below. See 4

_Preliminaries._ In what follows, we denote by
$\tilde{w}_{t}:=\mathop{\arg\min}_{w\in\mathcal{W}}F_{t}(w)$ the exact
solution of the inner-loop minibatch ERM optimization, which plays the same
role as $w_{t}$ in Section 2. We first present the following lemma that upper
bounds the discrepancy between the inexact minimizer $w_{t}$ and the exact
minimizer $\tilde{w}_{t}$.

######  Lemma 6.

Assume that the loss function $\ell$ is convex with respect to its first
argument and $r$ is convex. Then for any $w\in\mathcal{W}$, we have

$\|w_{t}-\tilde{w}_{t}\|\leq\sqrt{\frac{2\epsilon_{t}}{\gamma_{t}}}.$

###### Proof.

Using arguments identical to those of Lemma 1 we can show that for all
$w\in\mathcal{W}$,

$R_{S_{t}}(\tilde{w}_{t})-R_{S_{t}}(w)\leq\frac{\gamma_{t}}{2}\left(\|w-w_{t-1}\|^{2}-\|w-\tilde{w}_{t}\|^{2}-\|\tilde{w}_{t}-w_{t-1}\|^{2}\right).$
(15)

Setting $w=w_{t}$ in the above yields

$\frac{\gamma_{t}}{2}\|w_{t}-\tilde{w}_{t}\|^{2}\leq
F_{t}(w_{t})-F_{t}(\tilde{w}_{t})\leq\epsilon_{t},$

which directly implies
$\|w_{t}-\tilde{w}_{t}\|\leq\sqrt{2\epsilon_{t}/\gamma_{t}}$. This proves the
second desired bound. ∎

The following lemma as an extension of Lemma 3 to the inexact setting.

######  Lemma 7.

Suppose that the Assumptions 1, 2 and 3 hold. Assume that
$\gamma_{t}\geq\frac{19L}{n}$. Then the following bound holds for any
$\rho\in(0,1)$:

$\displaystyle\mathbb{E}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]\leq$
$\displaystyle\gamma_{t}\left(D^{2}(w_{t-1},W^{*})-\mathbb{E}\left[\left(1-\frac{\rho\lambda}{2\gamma_{t}}\right)D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]\right)$
$\displaystyle+\frac{19L}{\gamma_{t}n}R^{*}+\left(3n+\frac{4\gamma_{t}}{\rho\lambda}\right)\epsilon_{t}+3G\sqrt{\frac{2\epsilon_{t}}{\gamma_{t}}}.$

Alternatively, for any $w^{*}\in W^{*}$, under Assumptions 1 and 3 we have

$\displaystyle\mathbb{E}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]\leq$
$\displaystyle\gamma_{t}\left(\|w_{t-1}-w^{*}\|^{2}-\mathbb{E}\left[\|w_{t}-w^{*}\|^{2}\mid\mathcal{F}_{t-1}\right]\right)+\frac{19L}{\gamma_{t}n}R^{*}$
$\displaystyle+3n\epsilon_{t}+\left(2\sqrt{2\gamma_{t}}\mathbb{E}\left[\|w_{t}-w^{*}\|\mid\mathcal{F}_{t-1}\right]+3G\sqrt{\frac{2}{\gamma_{t}}}\right)\sqrt{\epsilon_{t}}.$

###### Proof.

Let us decompose $\mathbb{E}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]$
into the following three terms:

$\displaystyle\mathbb{E}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle=$
$\displaystyle\underbrace{\mathbb{E}\left[R(w_{t})-R(\tilde{w}_{t})\mid\mathcal{F}_{t-1}\right]}_{A}+\underbrace{\mathbb{E}\left[R(\tilde{w}_{t})-R_{S_{t}}(\tilde{w}_{t})\mid\mathcal{F}_{t-1}\right]}_{B}+\underbrace{\mathbb{E}\left[R_{S_{t}}(\tilde{w}_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]}_{C}.$

We next bound these three terms respectively. To bound the term $A$, we can
show that

$\displaystyle|A|:=|\mathbb{E}\left[R(w_{t})-R(\tilde{w}_{t})\mid\mathcal{F}_{t-1}\right]|$
$\displaystyle=$
$\displaystyle\left|\mathbb{E}\left[R^{\ell}(w_{t})-R^{\ell}(\tilde{w}_{t})\mid\mathcal{F}_{t-1}\right]+\mathbb{E}\left[r(w_{t})-r(\tilde{w}_{t})\right]\mid\mathcal{F}_{t-1}\right|$
$\displaystyle\leq$
$\displaystyle\mathbb{E}\left[\mathbb{E}_{z}|\ell(w_{t};z)-\ell(\tilde{w}_{t};z)|\mid\mathcal{F}_{t-1}\right]+\mathbb{E}\left[|r(w_{t})-r(\tilde{w}_{t})|\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\overset{\zeta_{1}}{\leq}$
$\displaystyle\mathbb{E}\left[\mathbb{E}_{z}\left[\sqrt{2L\ell(w_{t};z)}\|w_{t}-\tilde{w}_{t}\|\right]\mid\mathcal{F}_{t-1}\right]+\mathbb{E}\left[G\|w_{t}-\tilde{w}_{t}\|\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$
$\displaystyle\mathbb{E}\left[\mathbb{E}_{z}\left[\frac{L}{\gamma_{t}n}\ell(w_{t};z)+\frac{\gamma_{t}n}{2}\|w_{t}-\tilde{w}_{t}\|^{2}\right]\mid\mathcal{F}_{t-1}\right]+\mathbb{E}\left[G\|w_{t}-\tilde{w}_{t}\|\mid\mathcal{F}_{t-1}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\frac{L}{\gamma_{t}n}R^{\ell}(w_{t})\mid\mathcal{F}_{t-1}\right]+\mathbb{E}_{S_{t}}\left[\frac{\gamma_{t}n}{2}\|w_{t}-\tilde{w}_{t}\|^{2}+G\|w_{t}-\tilde{w}_{t}\|\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$
$\displaystyle\mathbb{E}\left[\frac{L}{\gamma_{t}n}R(w_{t})\mid\mathcal{F}_{t-1}\right]+n\epsilon_{t}+G\sqrt{\frac{2\epsilon_{t}}{\gamma_{t}}},$

where in “$\zeta_{1}$” we have used the convexity of loss and Lemma 2 and the
Assumption 3 and in the last inequality we have used $r>0$ and the
perturbation bound of Lemma 6.

To bound the term $B$, using about the same proof arguments as for Lemma 3 we
can show that

$\displaystyle B:=$
$\displaystyle\mathbb{E}\left[R(\tilde{w}_{t})-R_{S_{t}}(\tilde{w}_{t})\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$
$\displaystyle\frac{8L}{\gamma_{t}n}\mathbb{E}\left[R(\tilde{w}_{t})\mid\mathcal{F}_{t-1}\right]$
$\displaystyle=$
$\displaystyle\frac{8L}{\gamma_{t}n}\mathbb{E}\left[R(\tilde{w}_{t})-R(w_{t})\right]+\frac{8L}{\gamma_{t}n}\mathbb{E}\left[R(w_{t})\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$
$\displaystyle\frac{1}{2}|A|+\frac{8L}{\gamma_{t}n}\mathbb{E}\left[R(w_{t})\mid\mathcal{F}_{t-1}\right],$

where we have used the condition on minibatch size $\gamma_{t}$.

To bound the term $C$, based on the definition of $\tilde{w}_{t}$ and by
invoking Lemma 1 with $w=w^{*}_{t-1}$ we can verify that

$\displaystyle C:=$
$\displaystyle\mathbb{E}\left[R_{S_{t}}(\tilde{w}_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$
$\displaystyle\frac{\gamma_{t}}{2}\mathbb{E}\left[\|w^{*}_{t-1}-w_{t-1}\|^{2}-\|w^{*}_{t-1}-\tilde{w}_{t}\|^{2}-\|\tilde{w}_{t}-w_{t-1}\|^{2}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$
$\displaystyle\frac{\gamma_{t}}{2}\mathbb{E}\left[\|w^{*}_{t-1}-w_{t-1}\|^{2}-\|w^{*}_{t-1}-\tilde{w}_{t}\|^{2}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle=$
$\displaystyle\frac{\gamma_{t}}{2}\left(\|w^{*}_{t-1}-w_{t-1}\|^{2}-\mathbb{E}\left[\|w^{*}_{t-1}-w_{t}+w_{t}-\tilde{w}_{t}\|^{2}\mid\mathcal{F}_{t-1}\right]\right)$
$\displaystyle=$
$\displaystyle\frac{\gamma_{t}}{2}\left(\|w^{*}_{t-1}-w_{t-1}\|^{2}-\mathbb{E}\left[\|w^{*}_{t-1}-w_{t}\|^{2}+2\langle
w^{*}_{t-1}-w_{t},w_{t}-\tilde{w}_{t}\rangle+\|w_{t}-\tilde{w}_{t}\|^{2}\mid\mathcal{F}_{t-1}\right]\right)$
$\displaystyle\leq$
$\displaystyle\frac{\gamma_{t}}{2}\left(\|w^{*}_{t-1}-w_{t-1}\|^{2}-\mathbb{E}\left[\left(1-\frac{\rho\lambda}{2\gamma_{t}}\right)\|w^{*}_{t-1}-w_{t}\|^{2}-\frac{2\gamma_{t}}{\rho\lambda}\|w_{t}-\tilde{w}_{t}\|^{2}\mid\mathcal{F}_{t-1}\right]\right)$
$\displaystyle\leq$
$\displaystyle\frac{\gamma_{t}}{2}\left(\|w^{*}_{t-1}-w_{t-1}\|^{2}-\mathbb{E}\left[\left(1-\frac{\rho\lambda}{2\gamma_{t}}\right)\|w^{*}_{t}-w_{t}\|^{2}\mid\mathcal{F}_{t-1}\right]\right)+\frac{2\gamma_{t}\epsilon_{t}}{\rho\lambda}$
$\displaystyle=$
$\displaystyle\frac{\gamma_{t}}{2}\left(D^{2}(w_{t-1},W^{*})-\mathbb{E}\left[\left(1-\frac{\rho\lambda}{2\gamma_{t}}\right)D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]\right)+\frac{2\gamma_{t}\epsilon_{t}}{\rho\lambda}.$

Combining the above three bounds yields

$\displaystyle\mathbb{E}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]=A+B+C$
$\displaystyle\leq$
$\displaystyle\frac{3}{2}|A|+\frac{8L}{\gamma_{t}n}\mathbb{E}\left[R(w_{t})\mid\mathcal{F}_{t-1}\right]$
$\displaystyle+\frac{\gamma_{t}}{2}\left(D^{2}(w_{t-1},W^{*})-\mathbb{E}\left[\left(1-\frac{\rho\lambda}{2\gamma_{t}}\right)D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]\right)+\frac{2\gamma_{t}\epsilon_{t}}{\rho\lambda}$
$\displaystyle\leq$
$\displaystyle\mathbb{E}\left[\frac{3L}{2\gamma_{t}n}R(w_{t})\mid\mathcal{F}_{t-1}\right]+\frac{3n}{2}\epsilon_{t}+\frac{3G}{2}\sqrt{\frac{2\epsilon_{t}}{\gamma_{t}}}+\frac{8L}{\gamma_{t}n}\mathbb{E}\left[R(w_{t})\mid\mathcal{F}_{t-1}\right]$
$\displaystyle+\frac{\gamma_{t}}{2}\left(D^{2}(w_{t-1},W^{*})-\mathbb{E}\left[\left(1-\frac{\rho\lambda}{2\gamma_{t}}\right)D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]\right)+\frac{2\gamma_{t}\epsilon_{t}}{\rho\lambda}$
$\displaystyle\leq$
$\displaystyle\frac{\gamma_{t}}{2}\left(D^{2}(w_{t-1},W^{*})-\mathbb{E}\left[\left(1-\frac{\rho\lambda}{2\gamma_{t}}\right)D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]\right)+\frac{9.5L}{\gamma_{t}n}\mathbb{E}\left[R(w_{t})\mid\mathcal{F}_{t-1}\right]$
$\displaystyle+\left(\frac{3n}{2}+\frac{2\gamma_{t}}{\rho\lambda}\right)\epsilon_{t}+\frac{3G}{2}\sqrt{\frac{2\epsilon_{t}}{\gamma_{t}}}$
$\displaystyle=$
$\displaystyle\frac{\gamma_{t}}{2}\left(D^{2}(w_{t-1},W^{*})-\mathbb{E}\left[\left(1-\frac{\rho\lambda}{2\gamma_{t}}\right)D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]\right)+\frac{9.5L}{\gamma_{t}n}\mathbb{E}\left[R^{*}\right]+\frac{9.5L}{\gamma_{t}n}\mathbb{E}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle+\left(\frac{3n}{2}+\frac{2\gamma_{t}}{\rho\lambda}\right)\epsilon_{t}+\frac{3G}{2}\sqrt{\frac{2\epsilon_{t}}{\gamma_{t}}}$
$\displaystyle\leq$
$\displaystyle\frac{\gamma_{t}}{2}\left(D^{2}(w_{t-1},W^{*})-\mathbb{E}\left[\left(1-\frac{\rho\lambda}{2\gamma_{t}}\right)D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]\right)+\frac{9.5L}{\gamma_{t}n}\mathbb{E}\left[R^{*}\right]+\frac{1}{2}\mathbb{E}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle+\left(\frac{3n}{2}+\frac{2\gamma_{t}}{\rho\lambda}\right)\epsilon_{t}+\frac{3G}{2}\sqrt{\frac{2\epsilon_{t}}{\gamma_{t}}},$

where in the last inequality we have used the condition
$\gamma_{t}\geq\frac{19L}{n}$. After rearranging the terms in the above
inequality we obtain the first desired bound.

To derive the second bound, for any fixed $w^{*}\in W^{*}$, we note that the
term $C$ can be alternatively bounded as

$\displaystyle C\leq$
$\displaystyle\frac{\gamma_{t}}{2}\left(\|w^{*}-w_{t-1}\|^{2}-\mathbb{E}\left[\|w^{*}-w_{t}\|^{2}+2\langle
w^{*}-w_{t},w_{t}-\tilde{w}_{t}\rangle+\|w_{t}-\tilde{w}_{t}\|^{2}\mid\mathcal{F}_{t-1}\right]\right)$
$\displaystyle\leq$
$\displaystyle\frac{\gamma_{t}}{2}\left(\|w^{*}-w_{t-1}\|^{2}-\mathbb{E}\left[\|w^{*}-w_{t}\|^{2}-2\|w_{t}-w^{*}\|\|w_{t}-\tilde{w}_{t}\|\mid\mathcal{F}_{t-1}\right]\right)$
$\displaystyle\leq$
$\displaystyle\frac{\gamma_{t}}{2}\left(\|w^{*}-w_{t-1}\|^{2}-\mathbb{E}\left[\|w^{*}-w_{t}\|^{2}\mid\mathcal{F}_{t-1}\right]\right)+\sqrt{2\gamma_{t}\epsilon_{t}}\mathbb{E}\left[\|w^{*}-w_{t}\|\mid\mathcal{F}_{t-1}\right].$

Similar to the proof of the first bound, we can derive that

$\displaystyle\mathbb{E}_{S_{t}}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]=A+B+C$
$\displaystyle\leq$
$\displaystyle\frac{3}{2}|A|+\frac{8L}{\gamma_{t}n}\mathbb{E}\left[R(w_{t})\mid\mathcal{F}_{t-1}\right]+\frac{\gamma_{t}}{2}\left(\|w^{*}-w_{t-1}\|^{2}-\mathbb{E}\left[\|w^{*}-w_{t}\|^{2}\mid\mathcal{F}_{t-1}\right]\right)+\sqrt{2\gamma_{t}\epsilon_{t}}\mathbb{E}\left[\|w^{*}-w_{t}\|\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$
$\displaystyle\frac{\gamma_{t}}{2}\left(\|w^{*}-w_{t-1}\|^{2}-\mathbb{E}\left[\|w^{*}-w_{t}\|^{2}\mid\mathcal{F}_{t-1}\right]\right)+\frac{9.5L}{\gamma_{t}n}R^{*}$
$\displaystyle+\frac{1}{2}\mathbb{E}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]+\frac{3n}{2}\epsilon_{t}+\sqrt{2\gamma_{t}\epsilon_{t}}\mathbb{E}\left[\|w^{*}-w_{t}\|\mid\mathcal{F}_{t-1}\right]+\frac{3G}{2}\sqrt{\frac{2\epsilon_{t}}{\gamma_{t}}}.$

After rearranging the terms in the above inequality we obtain the second
desired bound. ∎

With the above preliminary results in hand, we are now in the position to
prove the main result of Theorem 4.

###### Proof of Theorem 4.

Since by assumption $R(w_{t})-R^{*}\geq\frac{\lambda}{2}D^{2}(w_{t},W^{*})$
and $\gamma_{t}=\frac{\lambda\rho
t}{4}\geq\frac{\lambda\rho}{4}\geq\frac{19L}{n}$, based on the first bound in
Lemma 7 we can show that

$\displaystyle(1-2\rho)\mathbb{E}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$
$\displaystyle\gamma_{t}D^{2}(w_{t-1},W^{*})-\left(\gamma_{t}+\frac{\rho\lambda}{2}\right)\mathbb{E}\left[D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]+\frac{19L}{\gamma_{t}n}R^{*}+\left(3n+\frac{4\gamma_{t}}{\rho\lambda}\right)\epsilon_{t}+3G\sqrt{\frac{2\epsilon_{t}}{\gamma_{t}}}$
$\displaystyle\leq$ $\displaystyle\frac{\lambda\rho
t}{4}D^{2}(w_{t-1},W^{*})-\frac{\rho\lambda(t+2)}{4}\mathbb{E}\left[D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]+\frac{76L}{\lambda\rho
nt}R^{*}+\left(3n+t\right)\epsilon_{t}+6G\sqrt{\frac{2\epsilon_{t}}{\lambda\rho
t}}.$

Now suppose that $\epsilon_{t}\leq\frac{\epsilon}{nt^{4}}$ for some
$\epsilon\in[0,1]$. Since $\rho\leq 1/4$, the above implies

$\displaystyle\mathbb{E}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$ $\displaystyle\frac{\lambda\rho
t}{2}D^{2}(w_{t-1},W^{*})-\frac{\rho\lambda(t+2)}{2}\mathbb{E}\left[D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]+\frac{152L}{\lambda\rho
nt}R^{*}+\left(\frac{6}{t^{4}}+\frac{2}{t^{3}}+12G\sqrt{\frac{2}{\lambda\rho
t^{5}}}\right)\sqrt{\epsilon}.$

The above inequality then implies

$\displaystyle
t\mathbb{E}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]\leq$
$\displaystyle(t+1)\mathbb{E}\left[R(w_{t})-R^{*}\mid\mathcal{F}_{t-1}\right]$
$\displaystyle\leq$ $\displaystyle\frac{\lambda\rho
t(t+1)}{2}D^{2}(w_{t-1},W^{*})-\frac{\lambda\rho(t+1)(t+2)}{2}\mathbb{E}\left[D^{2}(w_{t},W^{*})\mid\mathcal{F}_{t-1}\right]$
$\displaystyle+\frac{304L}{\lambda\rho
n}R^{*}+\left(\frac{12}{t^{3}}+\frac{4}{t^{2}}+\frac{24G}{t}\sqrt{\frac{2}{\lambda\rho
t}}\right)\sqrt{\epsilon},$

where we have used the fact $\frac{t+1}{t}\leq 2$ for $t\geq 1$. In view of
the law of total expectation, summing the above inequality from $t=1,...,T$
with natural normalization yields

$\displaystyle\frac{2}{T(T+1)}\sum_{t=1}^{T}t\mathbb{E}\left[R(w_{t})-R^{*}\right]$
$\displaystyle\leq$
$\displaystyle\frac{2\lambda\rho}{T(T+1)}D^{2}(w_{0},W^{*})+\frac{608L}{\lambda\rho(T+1)n}R^{*}+\frac{\sqrt{\epsilon}}{T(T+1)}\left(64+192G\sqrt{\frac{2}{\lambda\rho}}\right)$
$\displaystyle\leq$
$\displaystyle\frac{4\rho}{T(T+1)}(R(w_{0})-R^{*})+\frac{608L}{\lambda\rho(T+1)n}R^{*}+\frac{\sqrt{\epsilon}}{T(T+1)}\left(64+192G\sqrt{\frac{2}{\lambda\rho}}\right),$

which then immediately leads to the desired bound. The proof is concluded. ∎

###  B.2 Proof of Theorem 5

In this subsection, we prove Theorem 5 as following restated. See 5

The following lemma, which can be proved by induction (see, e.g., Schmidt et
al., 2011), will be used to prove the main result.

######  Lemma 8.

Assume that the nonnegative sequence $\\{u_{\tau}\\}_{\tau\geq 1}$ satisfies
the following recursion for all $t\geq 1$:

$u^{2}_{t}\leq S_{t}+\sum_{\tau=1}^{t}\alpha_{\tau}u_{\tau},$

with $\\{S_{\tau}\\}_{\tau\geq 1}$ an increasing sequence, $S_{0}\geq
u_{0}^{2}$ and $\alpha_{\tau}\geq 0$ for all $\tau$. Then, the following bound
holds for all $t\geq 1$:

$u_{t}\leq\sqrt{S_{t}}+\sum_{\tau=1}^{t}\alpha_{\tau}.$

The following lemma gives an upper bound on the expected estimation error
$\mathbb{E}\left[\|w^{*}_{0}-w_{t}\|\right]$.

######  Lemma 9.

Under the conditions of Theorem 5, the following bound holds for all $t\geq
1$:

$\mathbb{E}\left[\|w_{t}-w^{*}_{0}\|\right]\leq\|w_{0}-w^{*}_{0}\|+\sqrt{\frac{t}{\gamma}R^{*}}+\frac{6tG}{\gamma}.$

###### Proof.

Recall that $w^{*}_{0}=\mathop{\arg\min}_{w\in W^{*}}\|w_{0}-w\|$. Since
$\gamma_{t}\equiv\gamma\geq\frac{19L}{n}$, the second bound in Lemma 7 is
valid. For any $t\in[T]$, by summing that inequality with $w^{*}=w^{*}_{0}$
from $\tau=1,...,t$ we obtain

$\displaystyle\sum_{\tau=1}^{t}\mathbb{E}\left[R(w_{\tau})-R^{*}\right]+\gamma\mathbb{E}\left[\|w_{t}-w^{*}_{0}\|^{2}\right]$
(16) $\displaystyle\leq$
$\displaystyle\gamma\|w_{0}-w^{*}_{0}\|^{2}+\frac{19L}{\gamma
n}tR^{*}+3n\sum_{\tau=1}^{t}\epsilon_{\tau}+\sum_{\tau=1}^{t}\left(2\sqrt{2\gamma}\mathbb{E}\left[\|w^{*}_{0}-w_{\tau}\|\right]+3G\sqrt{\frac{2}{\gamma}}\right)\sqrt{\epsilon_{\tau}}.$

Dropping the non-negative term
$\sum_{\tau=1}^{t}\mathbb{E}_{S_{[\tau]}}\left[R(w_{\tau})-R^{*}\right]$ from
the above inequality yields

$\displaystyle\underbrace{\mathbb{E}\left[\|w_{t}-w^{*}_{0}\|^{2}\right]}_{u_{t}^{2}}$
$\displaystyle\leq$
$\displaystyle\|w_{0}-w^{*}_{0}\|^{2}+\frac{19L}{\gamma^{2}n}tR^{*}+\frac{3n}{\gamma}\sum_{\tau=1}^{t}\epsilon_{\tau}+\sum_{\tau=1}^{t}\left(2\sqrt{\frac{2}{\gamma}}\mathbb{E}\left[\|w^{*}_{0}-w_{\tau}\|\right]+3G\frac{\sqrt{2}}{\gamma\sqrt{\gamma}}\right)\sqrt{\epsilon_{\tau}}$
$\displaystyle\overset{\zeta_{1}}{\leq}$
$\displaystyle\|w_{0}-w^{*}_{0}\|^{2}+\frac{t}{\gamma}R^{*}+\sum_{\tau=1}^{t}\left(\frac{3n}{\gamma}\epsilon_{\tau}+\frac{3G\sqrt{2}}{\gamma\sqrt{\gamma}}\sqrt{\epsilon_{\tau}}\right)+\sum_{\tau=1}^{t}\left(2\sqrt{\frac{2\epsilon_{\tau}}{\gamma}}\sqrt{\mathbb{E}\left[\|w^{*}_{0}-w_{\tau}\|^{2}\right]}\right)$
$\displaystyle\leq$
$\displaystyle\underbrace{\|w_{0}-w^{*}_{0}\|^{2}+\frac{t}{\gamma}R^{*}+\sum_{\tau=1}^{t}\frac{4G\sqrt{2\epsilon_{\tau}}}{\gamma\sqrt{\gamma}}}_{S_{t}}+\sum_{\tau=1}^{t}\left(\underbrace{2\sqrt{\frac{2\epsilon_{\tau}}{\gamma}}}_{\alpha_{\tau}}\underbrace{\sqrt{\mathbb{E}\left[\|w^{*}_{0}-w_{\tau}\|^{2}\right]}}_{u_{\tau}}\right),$

where in “$\zeta_{1}$” we have used $\gamma\geq\frac{19L}{n}$ and the basic
inequality $\mathbb{E}^{2}[X]\leq\mathbb{E}[X^{2}]$, and in the last
inequality we have used the condition
$\epsilon_{\tau}\leq\frac{2G^{2}}{9n^{2}\gamma}$ for all $\tau\geq 1$. By
invoking Lemma 8 to the above recursion form we can derive that for all $t\geq
1$,

$\displaystyle\sqrt{\mathbb{E}\left[\|w_{t}-w^{*}_{0}\|^{2}\right]}\leq$
$\displaystyle\sqrt{\|w_{0}-w^{*}_{0}\|^{2}+\frac{t}{\gamma}R^{*}+\sum_{\tau=1}^{t}\frac{4G\sqrt{2\epsilon_{\tau}}}{\gamma}}+\sum_{\tau=1}^{t}2\sqrt{\frac{2\epsilon_{\tau}}{\gamma}}$
$\displaystyle\leq$
$\displaystyle\|w_{0}-w^{*}_{0}\|+\sqrt{\frac{t}{\gamma}R^{*}}+\sum_{\tau=1}^{t}\sqrt{\frac{4G\sqrt{2\epsilon_{\tau}}}{\gamma\sqrt{\gamma}}}+\sum_{\tau=1}^{t}2\sqrt{\frac{2\epsilon_{\tau}}{\gamma}}$
$\displaystyle\leq$
$\displaystyle\|w_{0}-w^{*}_{0}\|+\sqrt{\frac{t}{\gamma}R^{*}}+\frac{6Gt}{\gamma},$

where the last inequality is due to the condition
$\epsilon_{\tau}\leq\frac{2G^{2}}{9\gamma}$ for all $\tau\geq 1$. The above
inequality then directly implies the desired bound for all $t\in[T]$. ∎

Now we are ready to prove the main result of Theorem 5.

###### Proof of Theorem 5.

Dropping non-negative term $\gamma\mathbb{E}\left[\|w_{t}-w^{*}\|^{2}\right]$
in (16) followed by natural normalization yields

$\begin{aligned}
&\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\left[R(w_{t})-R^{*}\right]\\\
\leq&\frac{\gamma}{T}\|w_{0}-w^{*}_{0}\|^{2}+\frac{19L}{\gamma
n}R^{*}+\frac{3n}{T}\sum_{t=1}^{T}\epsilon_{t}+\frac{1}{T}\sum_{t=1}^{T}\left(2\sqrt{2\gamma}\mathbb{E}\left[\|w_{t}-w^{*}_{0}\|\right]+3G\sqrt{\frac{2}{\gamma}}\right)\sqrt{\epsilon_{t}}\\\
\overset{\zeta_{1}}{\leq}&\frac{\gamma}{T}\|w_{0}-w^{*}_{0}\|^{2}+\frac{19L}{\gamma
n}R^{*}+\frac{3n}{T}\sum_{t=1}^{T}\epsilon_{t}\\\
&+\frac{1}{T}\sum_{t=1}^{T}\left(2\sqrt{2}\left(\sqrt{\gamma}\|w_{0}-w^{*}_{0}\|+\sqrt{tR^{*}}+\frac{6Gt}{\sqrt{\gamma}}\right)+3G\sqrt{\frac{2}{\gamma}}\right)\sqrt{\epsilon_{t}}\\\
\overset{\zeta_{2}}{\leq}&\frac{\gamma}{T}\|w_{0}-w^{*}_{0}\|^{2}+\frac{19L}{\gamma
n}R^{*}+\frac{1}{T}\sum_{t=1}^{T}\left(3n\epsilon_{t}+2\sqrt{2\gamma\epsilon_{t}}\|w_{0}-w^{*}_{0}\|+2\sqrt{2tR^{*}\epsilon_{t}}+\frac{15\sqrt{2\epsilon_{t}}Gt}{\sqrt{\gamma}}\right)\\\
\overset{\zeta_{3}}{\leq}&\frac{\gamma}{T}\|w_{0}-w_{0}^{*}\|^{2}+\frac{19L}{\gamma
n}R^{*}\\\
&+\frac{1}{T}\sum_{t=1}^{T}\left(3n\epsilon_{t}+\frac{\gamma\|w_{0}-w^{*}_{0}\|^{2}}{t^{2}}+2t^{2}\epsilon_{t}+\frac{2LR^{*}}{\gamma
n}+\frac{\gamma
nt\epsilon_{t}}{L}+\frac{15\sqrt{2\epsilon_{t}}Gt}{\sqrt{\gamma}}\right)\\\
\leq&\frac{3\gamma}{T}\|w_{0}-w_{0}^{*}\|^{2}+\frac{21L}{\gamma
n}R^{*}+\frac{1}{T}\sum_{t=1}^{T}\left(3n\epsilon_{t}+2t^{2}\epsilon_{t}+\frac{\gamma
nt\epsilon_{t}}{L}+\frac{15\sqrt{2\epsilon_{t}}Gt}{\sqrt{\gamma}}\right),\end{aligned},$

where “$\zeta_{1}$”follows from Lemma 9, “$\zeta_{2}$” is due to $t\geq 1$ and
“$\zeta_{3}$” is due to $ab\leq(a^{2}+b^{2})/2$. Now consider
$\epsilon_{t}\leq\frac{\epsilon}{n^{2}t^{5}}$ for some $\epsilon\in[0,1]$.
Then it follows from the preceding inequality that

$\displaystyle\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\left[R(w_{t})-R^{*}\right]$
$\displaystyle\leq$
$\displaystyle\frac{3\gamma}{T}\|w_{0}-w_{0}^{*}\|^{2}+\frac{21L}{\gamma
n}R^{*}+\frac{1}{T}\sum_{t=1}^{T}\left(\frac{3}{nt^{5}}+\frac{2}{nt^{3}}+\frac{\gamma}{nLt^{4}}+\frac{15\sqrt{2}G}{nt^{1.5}\sqrt{\gamma}}\right)\sqrt{\epsilon}$
$\displaystyle\leq$
$\displaystyle\frac{3\gamma}{T}\|w_{0}-w_{0}^{*}\|^{2}+\frac{21L}{\gamma
n}R^{*}+\frac{1}{T}\left(\frac{6}{n}+\frac{4}{n}+\frac{2\gamma}{nL}+\frac{45\sqrt{2}G}{n\sqrt{\gamma}}\right)\sqrt{\epsilon}.$

Let $\bar{w}_{T}=\frac{1}{T}\sum_{t=1}^{T}w_{t}$. Combined with the convexity
of $R$, the above inequality implies

$\mathbb{E}\left[R(\bar{w}_{T})-R^{*}\right]\lesssim\frac{\gamma}{T}D^{2}(w_{0},W^{*})+\frac{L}{\gamma
n}R^{*}+\left(\frac{1}{nT}+\frac{\gamma}{LnT}+\frac{G}{nT\sqrt{\gamma}}\right)\sqrt{\epsilon}.$

This proves the first bound. Substituting
$\gamma=\sqrt{\frac{T}{n}}+\frac{19L}{n}>\frac{19L}{n}$ into the above bound
and preserving the leading terms yields the following second desired bound:

$\mathbb{E}\left[R(\bar{w}_{T})-R^{*}\right]\lesssim\left(\frac{1}{\sqrt{nT}}+\frac{L}{nT}\right)D^{2}(w_{0},W^{*})+\frac{L}{\sqrt{nT}}R^{*}+\left(\frac{L+G}{\sqrt{nT}}+\frac{1}{nT}\right)\sqrt{\epsilon}.$

The proof is concluded. ∎

##  Appendix C Proofs for the Results in Section 4

In this section, we present the proofs for the high probability estimation
error bounds stated in Section 4.

###  C.1 Proof of Proposition 1

In this subsection, we prove Proposition 1 as below restated . See 1

We first need to show the following preliminary result which is about the
expansion property of M-SPP update when performed over identical or different
minibatches.

######  Lemma 10.

Suppose that Assumptions 1 holds and the loss function $\ell$ is bounded in
the interval $[0,M]$. From $w_{0}=w^{\prime}_{0}$, let us define the sequences
$\\{w_{t}\\}_{t\in[T]}$ and $\\{w^{\prime}_{t}\\}_{t\in[T]}$ that are
respectively generated over $\\{S_{t}\\}_{t\in[T]}$ and
$\\{S^{\prime}_{t}\\}_{t\in[T]}$ according to

$\displaystyle F_{t}(w_{t})\leq$
$\displaystyle\min_{w\in\mathcal{W}}\left\\{F_{t}(w):=R_{S_{t}}(w)+\frac{\gamma_{t}}{2}\|w-w_{t-1}\|^{2}\right\\}+\epsilon_{t},$
$\displaystyle F^{\prime}_{t}(w^{\prime}_{t})\leq$
$\displaystyle\min_{w\in\mathcal{W}}\left\\{F^{\prime}_{t}(w):=R_{S^{\prime}_{t}}(w)+\frac{\gamma_{t}}{2}\|w-w^{\prime}_{t-1}\|^{2}\right\\}+\epsilon_{t}.$

Assume that either $S_{t}=S^{\prime}_{t}$ or $S_{t}\doteq S^{\prime}_{t}$ for
all $t\in[T]$. Let $\beta_{t}=\mathbf{1}_{\left\\{S_{t}\neq
S^{\prime}_{t}\right\\}}$. Then the following bound holds for all $t\in[T]$,

$\|w_{t}-w^{\prime}_{t}\|\leq\sum_{\tau=1}^{t}\left\\{\beta_{\tau}\frac{4\sqrt{2LM}}{n\gamma_{\tau}}+2\sqrt{\frac{2\epsilon_{\tau}}{\gamma_{\tau}}}\right\\}.$

###### Proof.

Let $w^{*}_{t}=\mathop{\arg\min}_{w}F_{t}(w)$ and
$w^{\prime*}_{t}=\mathop{\arg\min}_{w}F^{\prime}_{t}(w)$. It follows from
Lemma 1 that

$\displaystyle R_{S_{t}}(w^{*}_{t})-R_{S_{t}}(w^{\prime*}_{t})\leq$
$\displaystyle\frac{\gamma_{t}}{2}\left(\|w^{\prime*}_{t}-w_{t-1}\|^{2}-\|w^{\prime*}_{t}-w^{*}_{t}\|^{2}-\|w^{*}_{t}-w_{t-1}\|^{2}\right)$
$\displaystyle
R_{S^{\prime}_{t}}(w^{\prime*}_{t})-R_{S^{\prime}_{t}}(w^{*}_{t})\leq$
$\displaystyle\frac{\gamma_{t}}{2}\left(\|w^{*}_{t}-w^{\prime}_{t-1}\|^{2}-\|w^{\prime*}_{t}-w^{*}_{t}\|^{2}-\|w^{\prime*}_{t}-w^{\prime}_{t-1}\|^{2}\right).$

Summing both sides of the above two inequalities yields

$\begin{aligned}
&R_{S_{t}}(w^{*}_{t})-R_{S_{t}}(w^{\prime*}_{t})+R_{S^{\prime}_{t}}(w^{\prime*}_{t})-R_{S^{\prime}_{t}}(w^{*}_{t})\\\
\leq&\frac{\gamma_{t}}{2}\left(\|w^{\prime*}_{t}-w_{t-1}\|^{2}-\|w^{*}_{t}-w_{t-1}\|^{2}+\|w^{*}_{t}-w^{\prime}_{t-1}\|^{2}-\|w^{\prime*}_{t}-w^{\prime}_{t-1}\|^{2}-2\|w^{\prime*}_{t}-w^{*}_{t}\|^{2}\right)\\\
=&\frac{\gamma_{t}}{2}\left(2\langle
w^{*}_{t}-w^{\prime*}_{t},w_{t-1}-w^{\prime}_{t-1}\rangle-2\|w^{\prime*}_{t}-w^{*}_{t}\|^{2}\right)\\\
\leq&\frac{\gamma_{t}}{2}\left(\|w_{t-1}-w^{\prime}_{t-1}\|^{2}-\|w^{*}_{t}-w^{\prime*}_{t}\|^{2}\right)\end{aligned}.$

We need to distinguish the following two complementary cases.

Case I: $S_{t}=S^{\prime}_{t}$. In this case, the previous inequality
immediately leads to

$\|w^{*}_{t}-w^{\prime*}_{t}\|\leq\|w_{t-1}-w^{\prime}_{t-1}\|.$

By using triangle inequality and Lemma 6 we obtain

$\|w_{t}-w^{\prime}_{t}\|\leq\|w_{t}-w^{*}_{t}\|+\|w^{*}_{t}-w^{\prime*}_{t}\|+\|w^{\prime}_{t}-w^{\prime*}_{t}\|\leq\|w_{t-1}-w^{\prime}_{t-1}\|+2\sqrt{\frac{2\epsilon_{t}}{\gamma_{t}}}.$
(17)

Case II: $S_{t}$ and $S^{\prime}_{t}$ differ in a single element. In this
case, we have

$\displaystyle\|w^{*}_{t}-w^{\prime*}_{t}\|^{2}$ $\displaystyle\leq$
$\displaystyle\|w_{t-1}-w^{\prime}_{t-1}\|^{2}+\frac{2}{\gamma_{t}}\left(R_{S_{t}}(w^{\prime*}_{t})-R_{S_{t}}(w^{*}_{t})+R_{S^{\prime}_{t}}(w^{*}_{t})-R_{S^{\prime}_{t}}(w^{\prime*}_{t})\right)$
$\displaystyle=$
$\displaystyle\|w_{t-1}-w^{\prime}_{t-1}\|^{2}+\frac{2}{\gamma_{t}}\left(R^{\ell}_{S_{t}}(w^{\prime*}_{t})-R^{\ell}_{S_{t}}(w^{*}_{t})+R^{\ell}_{S^{\prime}_{t}}(w^{*}_{t})-R^{\ell}_{S^{\prime}_{t}}(w^{\prime*}_{t})\right)$
$\displaystyle=$
$\displaystyle\|w_{t-1}-w^{\prime}_{t-1}\|^{2}+\frac{2}{\gamma_{t}}\left(\frac{1}{|S_{t}|}\sum_{z\in
S_{t}}(\ell(w^{\prime*}_{t};z)-\ell(w^{*}_{t};z)+\frac{1}{|S^{\prime}_{t}|}\sum_{z\in
S^{\prime}_{t}}(\ell(w^{*}_{t};z)-\ell(w^{\prime*}_{t};z))\right)$
$\displaystyle\leq$
$\displaystyle\|w_{t-1}-w^{\prime}_{t-1}\|^{2}+\frac{4\sqrt{2LM}}{n\gamma_{t}}\|w^{*}_{t}-w^{\prime*}_{t}\|.$

where in the last inequality we have used $\ell(\cdot;\cdot)$ is
$\sqrt{2LM}$-Lipschitz with respect to its first argument which is implied by
Lemma 2, and $S_{t}$ and $S^{\prime}_{t}$ differ in a single element as well.
Since $x^{2}\leq y^{2}+ax$ implies $x\leq y+a$ for all $x,y,a>0$, we can
derive from the above that

$\|w^{*}_{t}-w^{\prime*}_{t}\|\leq\|w_{t-1}-w^{\prime}_{t-1}\|+\frac{4\sqrt{2LM}}{n\gamma_{t}}.$

Then based on triangle inequality and Lemma 6 we have

$\displaystyle\|w_{t}-w^{\prime}_{t}\|\leq\|w_{t}-w^{*}_{t}\|+\|w^{*}_{t}-w^{\prime*}_{t}\|+\|w^{\prime}_{t}-w^{\prime*}_{t}\|\leq\|w_{t-1}-w^{\prime}_{t-1}\|+\frac{4\sqrt{2LM}}{n\gamma_{t}}+2\sqrt{\frac{2\epsilon_{t}}{\gamma_{t}}}.$
(18)

Let $\beta_{t}=\mathbf{1}_{\left\\{S_{t}\neq S^{\prime}_{t}\right\\}}$ where
$\mathbf{1}_{\\{C\\}}$ is the indicator function of the condition $C$. Based
on the recursion forms (17) and (18) and the condition $w_{0}=w^{\prime}_{0}$
we can show that for all $t\in[T]$

$\|w_{t}-w^{\prime}_{t}\|\leq\sum_{\tau=1}^{t}\left\\{\frac{4\beta_{\tau}\sqrt{2LM}}{n\gamma_{\tau}}+2\sqrt{\frac{2\epsilon_{\tau}}{\gamma_{\tau}}}\right\\},$

which is the desired bound. ∎

Now we are in the position to prove the main result in Proposition 1.

###### Proof of Proposition 1.

Consider a fixed pair of minibatch sets $S\doteq S^{\prime}$.

Part (a): Let $\\{w_{t}\\}_{t\in[T]}$ and $\\{w^{\prime}_{t}\\}_{t\in[T]}$ be
two solution sequences that are respectively generated over
$\\{S_{t}\\}_{t\in[T]}$ and $\\{S^{\prime}_{t}\\}_{t\in[T]}$ by Algorithm 1.
At each time instance $t$, define random variable
$\beta_{t}:=\mathbf{1}_{\left\\{S_{t}\neq S^{\prime}_{t}\right\\}}$. Since by
assumption $S$ and $S^{\prime}$ differ only in a single minibatch, there must
exist one and only one $t\in[T]$ such that $\beta_{t}=1$ and $\beta_{j}=0$ for
all $j\in[T],j\neq t$. Then in the worst case of $\beta_{\tau}=1$ for
$\tau=\mathop{\arg\min}_{i\in[t]}\gamma_{i}$, it follows from Lemma 10 that
for all $t\in[T]$,

$\|w_{t}-w^{\prime}_{t}\|\leq\frac{4\sqrt{2LM}}{n\min_{i\in[t]}\gamma_{i}}+\sum_{i=1}^{t}2\sqrt{\frac{2\epsilon_{i}}{\gamma_{i}}}\leq\frac{4\sqrt{2LM}}{n\min_{i\in[T]}\gamma_{i}}+\sum_{i=1}^{T}2\sqrt{\frac{2\epsilon_{i}}{\gamma_{i}}}.$

Then the convex combination nature of $\bar{w}_{T}$ and $\bar{w}^{\prime}_{T}$
implies that

$\|\bar{w}_{T}-\bar{w}^{\prime}_{T}\|\leq\frac{\sum_{t}\gamma_{t}\|w_{t}-w^{\prime}_{t}\|}{\sum_{t}\gamma_{t}}\leq\frac{4\sqrt{2LM}}{n\min_{t\in[T]}\gamma_{t}}+\sum_{t=1}^{T}2\sqrt{\frac{2\epsilon_{t}}{\gamma_{t}}}.$

The desired result follows immediately as the above bound holds for any pair
$\\{S,S^{\prime}\\}$.

Part (b): Recall that $\\{\xi_{t}\\}_{t\in[T]}$ are the uniform random indices
for iteratively selecting data minibatches from $S$ and $S^{\prime}$. Let
$\\{w_{t}\\}_{t\in[T]}$ and $\\{w^{\prime}_{t}\\}_{t\in[T]}$ be two solution
sequences that are respectively generated over $\\{S_{\xi_{t}}\\}_{t\in[T]}$
and $\\{S^{\prime}_{\xi_{t}}\\}_{t\in[T]}$ by Algorithm 3. Define random
variable $\beta_{t}:=\mathbf{1}_{\left\\{S_{\xi_{t}}\neq
S^{\prime}_{\xi_{t}}\right\\}}$. Since by assumption $S$ and $S^{\prime}$
differ only in a single minibatch, under without-replacement sampling scheme,
there must exist one and only one $t\in[T]$ such that $\beta_{t}=1$ and
$\beta_{j}=0$ for all $j\in[T],j\neq t$. Let us define the event
$\mathcal{E}_{t}:=\\{\beta_{t}=1\text{ and }\beta_{j\neq t,j\in[T]}=0\\}$ for
all $t\in[T]$. Then the uniform randomness of $\xi_{t}$ implies that

$R\left(\mathcal{E}_{t}\right)=\frac{1}{T},\quad t\in[T].$

Given $t\in[T]$, suppose that $\mathcal{E}_{\tau}$ occurs for some
$\tau\in[t]$. Then it follows from Lemma 10 that

$\|w_{t}-w^{\prime}_{t}\|\leq\frac{4\sqrt{2LM}}{n\gamma_{\tau}}+\sum_{i=1}^{t}2\sqrt{\frac{2\epsilon_{i}}{\gamma_{i}}}.$

Suppose that $\mathcal{E}_{\tau}$ occurs for some
$\tau\in\\{t+1,t+2,...,T\\}$, again it follows from Lemma 10 that

$\|w_{t}-w^{\prime}_{t}\|\leq\sum_{i=1}^{t}2\sqrt{\frac{2\epsilon_{i}}{\gamma_{i}}}.$

Then we have

$\displaystyle\mathbb{E}_{\xi_{[t]}}\left[\|w_{t}-w^{\prime}_{t}\|\right]=$
$\displaystyle\sum_{\tau=1}^{T}R\left(\mathcal{E}_{\tau}\right)\left[\|w_{t}-w^{\prime}_{t}\|\mid\mathcal{E}_{\tau}\right]$
$\displaystyle\leq$
$\displaystyle\sum_{\tau=1}^{t}\left\\{\frac{4\sqrt{2LM}}{nT\gamma_{t}}+\sum_{i=1}^{t}\frac{2}{T}\sqrt{\frac{2\epsilon_{i}}{\gamma_{i}}}\right\\}+\sum_{\tau=t+1}^{T}\left\\{\sum_{i=1}^{t}\frac{2}{T}\sqrt{\frac{2\epsilon_{t}}{\gamma_{t}}}\right\\}$
$\displaystyle=$
$\displaystyle\sum_{\tau=1}^{t}\left\\{\frac{4\sqrt{2LM}}{nT\gamma_{\tau}}+2\sqrt{\frac{2\epsilon_{\tau}}{\gamma_{\tau}}}\right\\}\leq\sum_{t=1}^{T}\left\\{\frac{4\sqrt{2LM}}{nT\gamma_{t}}+2\sqrt{\frac{2\epsilon_{t}}{\gamma_{t}}}\right\\}.$

It follows that

$\mathbb{E}_{\xi_{[T]}}\left[\|\bar{w}_{T}-\bar{w}^{\prime}_{T}\|\right]\leq\frac{\sum_{t}\gamma_{t}\mathbb{E}_{\xi_{[t]}}\left[\|w_{t}-w^{\prime}_{t}\|\right]}{\sum_{t}\gamma_{t}}\leq\sum_{t=1}^{T}\left\\{\frac{4\sqrt{2LM}}{nT\gamma_{t}}+2\sqrt{\frac{2\epsilon_{t}}{\gamma_{t}}}\right\\}.$

The desired result follows immediately as the above bound holds for any pair
$\\{S,S^{\prime}\\}$. ∎

###  C.2 Proof of Theorem 6

In this subsection, we prove Theorem 6 that is restated below. See 6

To show this result, we need to use the following restated McDiarmid’s
inequality (McDiarmid, 1989) which is also known as bounded difference
inequality.

######  Lemma 11 (McDiarmid’s/Bounded differences inequality).

Let $X_{1},X_{2},...,X_{N}$ be independent random variables valued in
$\mathcal{X}$. Suppose that the function $h:\mathcal{X}^{N}\mapsto\mathbb{R}$
satisfies the bounded differences property, i.e., the following inequality
holds for any $i\in[N]$ and any $x_{1},...,x_{N},x^{\prime}_{i}$:

$|h(x_{1},...,x_{i-1},x_{i},x_{i+1},...,x_{N})-h(x_{1},...,x_{i-1},x^{\prime}_{i},x_{i+1},...,x_{N})|\leq
c_{i}.$

Then for any $\varepsilon>0$,

$\mathbb{P}\left(h(X_{1},...,X_{N})-\mathbb{E}\left[h(X_{1},...,X_{N})\right]\geq\varepsilon\right)\leq\exp\left(-\frac{2\varepsilon^{2}}{\sum_{i=1}^{N}c_{i}^{2}}\right).$

Now we are ready to prove Theorem 6.

###### Proof of Theorem 6.

Let $S=\\{S_{t}\\}_{t\in[T]}$ and $S^{\prime}=\\{S^{\prime}_{t}\\}_{t\in[T]}$
be two sets of data minibatches such that $S\doteq S^{\prime}$. Then according
to Proposition 1 the weighted average output $\bar{w}_{T}$ and
$\bar{w}^{\prime}_{T}$ respectively generated by Algorithm 3 over $S$ and
$S^{\prime}$ satisfy

$\sup_{S,S^{\prime}}\mathbb{E}_{\xi_{[T]}}\left[\|\bar{w}_{T}-\bar{w}^{\prime}_{T}\|\right]\leq\sum_{t=1}^{T}\left\\{\frac{4\sqrt{2LM}}{nT\gamma_{t}}+2\sqrt{\frac{2\epsilon_{t}}{\gamma_{t}}}\right\\}\leq\sum_{t=1}^{T}\left\\{\frac{5\sqrt{2LM}}{nT\gamma_{t}}\right\\}\leq\frac{20\sqrt{2LM}(1+\log(T))}{\lambda\rho
nT},$

where in the last but one inequality we have used the condition
$\epsilon_{t}\leq\frac{LM}{4n^{2}T^{2}\gamma_{t}}=\frac{LM}{\lambda\rho