# Stochastic optimization with arbitrary recurrent data sampling

William G. Powell  William G. Powell, Department of Mathematics, University of
Wisconsin - Madison, WI, 53706, USA<EMAIL_ADDRESS>and  Hanbaek Lyu
Hanbaek Lyu, Department of Mathematics, University of Wisconsin - Madison, WI,
53706, USA<EMAIL_ADDRESS>
###### Abstract.

For obtaining optimal first-order convergence guarantee for stochastic
optimization, it is necessary to use a recurrent data sampling algorithm that
samples every data point with sufficient frequency. Most commonly used data
sampling algorithms (e.g., i.i.d., MCMC, random reshuffling) are indeed
recurrent under mild assumptions. In this work, we show that for a particular
class of stochastic optimization algorithms, we do not need any other property
(e.g., independence, exponential mixing, and reshuffling) than recurrence in
data sampling algorithms to guarantee the optimal rate of first-order
convergence. Namely, using regularized versions of Minimization by Incremental
Surrogate Optimization (MISO), we show that for non-convex and possibly non-
smooth objective functions, the expected optimality gap converges at an
optimal rate $O(n^{-1/2})$ under general recurrent sampling schemes.
Furthermore, the implied constant depends explicitly on the ‘speed of
recurrence’, measured by the expected amount of time to visit a given data
point either averaged (‘target time’) or supremized (‘hitting time’) over the
current location. We demonstrate theoretically and empirically that
convergence can be accelerated by selecting sampling algorithms that cover the
data set most effectively. We discuss applications of our general framework to
decentralized optimization and distributed non-negative matrix factorization.

††footnotetext: codes available at https://github.com/wgraysonp/RMISO

##  1\. Introduction

In this paper, we consider the minimization of a non-convex continuous
function $f:\mathbb{R}^{p}\to\mathbb{R}$:

$\displaystyle\boldsymbol{\theta}^{*}\in\operatorname*{arg\,min}_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\Big{\\{}f(\boldsymbol{\theta}):=\mathbb{E}_{v\sim\pi}[f^{v}(\boldsymbol{\theta})]\Big{\\}}$
(1)

where $\pi$ is a probability distribution supported on a finite set
$\mathcal{V}$ and $\boldsymbol{\Theta}\subseteq\mathbb{R}^{p}$ is a convex,
but not necessarily compact feasible set. Here $\boldsymbol{\theta}$
represents the parameters of a model to be optimized, $\mathcal{V}$ is a
collection of data, and $f^{v}(\boldsymbol{\theta})$ is the loss incurred
using parameter $\boldsymbol{\theta}$ with respect to data $v$. We aim to
solve this problem by developing an algorithm that produces iterative
parameter updates $\boldsymbol{\theta}_{n}$ given only access to an arbitrary
sequence of data samples $(v_{n})_{n\geq 1}$. In order to reach a first-order
stationary point of (1) for general objectives, it is necessary to use a
sampling algorithm that is recurrent, meaning that every data point is sampled
infinitely often with ‘sufficient frequency’. Note that recurrence is
satisfied by many common sampling algorithms such as i.i.d., (irreducible)
Markovian, and cyclic [4] and random-reshuffling [57]. The main question we
ask in this work is the following:

$\bullet$ Is there any class of stochastic optimization algorithms for which
recurrent sampling is enough to obtain optimal first-order convergence
guarantee for (1)?

In this paper, we show that for a class of suitable extensions of stochastic
optimization algorithms known as _Minimization by Incremental Surrogate
Optimization_ (MISO) [36], no additional property of a data sampling algorithm
(e.g., independence, exponential mixing, reshuffling) other than recurrence is
needed in order to guarantee convergence to first-order stationary points.
Furthermore, we show that the rate of convergence depends crucially on the
average hitting time between data points.

MISO [36] was originally developed to solve finite sum problems under i.i.d
sampling and proceeds by repeatedly minimizing a surrogate of the empirical
loss function. In [36] it is shown that for MISO the expected objective
optimality gap
$\mathbb{E}[f(\boldsymbol{\theta}_{n})-f(\boldsymbol{\theta}^{*})]$ decays at
rate $O(1/n)$ when the objective function is convex and exponentially fast
when it is strongly convex, just as batch gradient descent does [3]. In this
work, we adapt MISO for solving the problem (1) and provide convergence rate
analysis in the _non-convex_ setting with general recurrent data sampling. We
call our algorithm _Regularized Minimization by Incremental Surrogate
Optimization_ (RMISO), which takes the following form:

$\displaystyle\textbf{(RMISO)}:\begin{cases}\text{Sample }v_{n}\text{ from the
conditional distribution }\pi(\cdot|\mathcal{F}_{n-1})\\\
g_{n}^{v_{n}}\leftarrow\text{Convex majorizing surrogate of $f^{v_{n}}$ at
$\boldsymbol{\theta}_{n-1}$; $g_{n}^{v}=g_{n-1}^{v}$ for $v\neq v_{n}$}\\\
\boldsymbol{\theta}_{n}\in\operatorname*{arg\,min}_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left(\bar{g}_{n}(\boldsymbol{\theta})+\Psi(\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert):=\sum_{v\in\mathcal{V}}g_{n}^{v}(\boldsymbol{\theta})\pi(v)+\Psi(\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert)\right).\end{cases}$
(2)

The algorithm maintains a list of majorizing surrogate functions for each data
point $v$. At each step, a new data sample $v_{n}$ is drawn according to a
recursive data sampling algorithm. We then find a new majorizing convex
surrogate $g_{n}^{v_{n}}$ that is tight at the current parameter
$\boldsymbol{\theta}_{n-1}$. All other surrogates are unchanged. Then the new
parameter $\boldsymbol{\theta}_{n}$ is found by minimizing the empirical mean
of the current surrogates plus a regularization term
$\Psi(\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert)$ that
penalizes large values of
$\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert$. The original
MISO [36] is recovered by omitting this regularization term. The particular
choice of the regularization term is crucial for the success of the analysis
under the general recurrent data sampling setting.

MISO and RMISO may be compared with _Stochastic Averaged Gradient_ (SAG) [53]
as both store the most recent information computed using the data $v$ and
output new parameter updates $\boldsymbol{\theta}_{n}$ depending on an average
of this information over $\mathcal{V}$. Recently in [15], it was shown that
for non-convex objectives SAG produces iterates such that the expected
gradient norm $\mathbb{E}\lVert\nabla f(\boldsymbol{\theta}_{n})\rVert$ decays
at rate $O(n^{-1/2})$ under Markovian sampling. We analyze two extensions of
MISO: one that matches this rate exactly and another up to logarithmic
factors. In comparison, the expected gradient norm converges at rate
$O(n^{-1/4})$ for other stochastic first-order methods such as Stochastic
Gradient Descent (SGD) [55, 2, 15, 25], see Table 1 for a summary. Therefore,
while (R)MISO requires additional storage that scales linearly with
$|\mathcal{V}|$, they enjoy much faster convergence rates.

We remark that the additional storage used by RMISO is the primary reason
simple recurrent sampling is enough to guarantee convergence, and why we
expect the same should hold for SAG. Because RMISO maintains a surrogate for
each $f^{v}$, updating each surrogate frequently enough ensures that the
_averaged_ surrogate $\bar{g}_{n}$ remains a good approximation of the full
objective $f$ at the most recent iterate. In this way, the trajectory of
iterates produced by RMISO stays close to that which would result from
minimizing a fully updated surrogate of $f$ at each step. In contrast,
recurrence alone may not be enough for other first-order methods like SGD.
Without the additional memory, more assumptions may be required to ensure that
stochastic gradients are sufficiently concentrated around their mean so that
the trajectory of iterates stays close to that of gradient descent.

Applications of our work include distributed optimization over networked
systems where $\mathcal{V}$ forms the vertex set of a connected graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ and each vertex $v$ stores some data.
Prior work [22, 21, 51, 32, 42, 15] studies the performance of various
optimization algorithms in this setting with the assumption that data is only
exchanged between neighboring vertices in order to minimize potentially
expensive communication costs. In particular, they assume that the sequence
$(v_{n})_{n\geq 1}$ is generated by a Markov chain on the graph $\mathcal{G}$.
In this setting, we find both theoretically and empirically that convergence
of our algorithm can be accelerated by choosing sampling schemes that
guarantee a higher frequency of visits to each $v\in\mathcal{V}$. Such schemes
may be non-Markovian or not aperiodic and so not guaranteed to converge to a
stationary distribution.

###  1.1. Contribution

Our algorithms and analysis consider three cases: 1) convex surrogates with
constant proximal regularization; 2) convex surrogates with dynamic proximal
regularization; 3) convex surrogates with diminishing radius restriction. We
briefly summarize our results in the following bullet points.

  * MISO with convex surrogates and constant proximal regularization

    * •

MISO with constant quadratic proximal regularization is developed.

    * •

Convergence rates of $O(n^{-1/2})$ are shown for this setting with optimal
dependence on constants.

  * MISO with convex surrogates and dynamic proximal regularization

    * •

MISO with dynamic proximal regularization is developed. Inspired by the
dynamic step size used for SAG in [15], the regularization parameter is
adaptive to the state of the sampling process.

    * •

The same convergence rate results as in the previous case hold.

    * •

Asymptotic convergence of stationarity measures in expectation is proved.

  * MISO with convex surrogates and diminishing radius restriction

    * •

MISO with diminishing search radius restriction is developed, where averaged
surrogates are minimized within a diminishing radius.

    * •

Convergence rates of $O(n^{-1/2}\log n)$ are shown.

    * •

Almost sure convergence to stationarity is proved.

In Section 5.2, we experimentally validate our algorithm for the tasks of non-
negative matrix factorization and binary classification. Empirically, we find
that our method is robust to data heterogeneity as it produces stable iterate
trajectories while still maintaining fast convergence. This property can be
expected as most implementations of our general framework produce updates
based on averages of previous gradients _and_ previous iterates. In this
sense, our algorithm shares characteristics with the ’double-averaging PSGD’
scheme first investigated by Nesterov and Shikhman in [44]. As the authors of
this paper discuss, a common strategy for handling possibly unstable sequences
of iterates is to consider the subsequence that produces the best value of the
objective function seen so far. However, this strategy fails if the value of
the objective function cannot easily be computed. See [44, Sect. 4.3] for an
example of such a situation. In such a case, a more stable algorithm may be
advantageous.

###  1.2. Related Work

In [36], MISO with i.i.d data is studied. It is shown that when the objective
function $f$ is convex, MISO decreases the function value optimality gap at
rate $O(1/n)$ and exponentially fast when $f$ is strongly convex. It is also
shown that when $f$ is non-convex, the iterates generated by MISO converge to
the set of stationary points of $f$ over a convex constraint set, but no
convergence rate analysis is given.

Convergence rate analysis for MISO in the unconstrained non-convex setting was
given in [48] where the authors analyzed a mini-batch version of the algorithm
and obtained convergence rates of $O(n^{-1/2})$. The authors of [26] study a
doubly stochastic extension of MISO in which the surrogate functions are given
as the expectation of a class of random functions with respect to an
underlying distribution. They propose a method in which the surrogates are
estimated via Monte-Carlo approximation and obtain and also obtain convergence
rates of $O(n^{-1/2})$ for non-convex objectives in the constrained setting.
However, both papers assume i.i.d sampling. In contrast, our work focuses on a
much broader class of possibly dependent sampling algorithms.

A large amount of recent work has been devoted to the study of first-order
optimization methods under Markovian sampling (e.g.[8, 7, 56, 17, 2, 34]). In
the non-convex case, [55, Thm. 2] shows that the expected gradient norm
$\mathbb{E}[\lVert\nabla f(\boldsymbol{\theta}_{n})\rVert]$ decays at rate
$\tilde{O}(n^{-1/4})$. First-order optimization with Markovian data has been
further studied in a number of works such as [8, 7, 56, 17, 2, 34] and those
especially concerned with finite sum problems include [21, 22, 42, 15, 54]. In
particular, [22, Thm. 1] shows asymptotic convergence of projected subgradient
descent to the optimal objective value. The authors of [54] consider a variant
of SGD with an adaptive stepsize. They prove that their algorithm requires at
most $O(\varepsilon^{-\frac{4}{2-\nu}})$ iterations to find
$\boldsymbol{\theta}_{n}$ satisfying $\mathbb{E}[\lVert\nabla
f(\boldsymbol{\theta}_{n})\rVert]=O(\varepsilon)$ where $\nu\in(0,1]$ is a
parameter depending on the sparsity of the gradients. As mentioned in the
introduction, [15] shows that for SAG, the expected gradient norm
$\mathbb{E}[\lVert\nabla f(\boldsymbol{\theta}_{n})\rVert]$ decays at rate
$O(n^{-1/2})$ under Markovian sampling schemes.

There has also been a recent focus on proving faster convergence for SGD using
without-replacement sampling methods such as random-reshuffling [16, 57]. This
has been further extended to variance-reduced algorithms [20, 41, 10] and
distributed optimization [39, 18]. New sampling algorithms that aim to improve
over random reshuffling have been suggested in [50, 27, 40]. In particular, in
[29] the authors show that the convergence of SGD can be accelerated provided
a certain concentration inequality holds and propose leveraging this using a
greedy selection strategy.

To obtain our results, we need to adopt a new analytical approach which is
inspired in part by the analysis of SAG in [15]. This strategy differs
significantly from mixing rate arguments used in the analysis of stochastic
optimization methods with Markovian data (e.g [55, 7, 45, 30, 33, 34, 2]). We
believe that these techniques may be of interest in their own right and may
further contribute to analyzing other stochastic optimization methods with
recurrent data sampling.

###  1.3. Organization

In Section 2 we formally define first-order surrogate functions as well as
first return time and last passage time. We then state our main algorithms in
Algorithms 1 and 2. Next, in Section 3 we first define our optimality
conditions and discuss our assumptions. We then give our main results in
Theorems 3.1 and 3.4 as well as Corollaries 3.2 and 3.3. In Section 5 we
discuss specific implementations of our framework and their applications -
distributed matrix factorization and RMISO with prox-linear surrogates. We
provide numerical experiments for the tasks of distributed non-negative matrix
factorization and binary classification in Section 5.2.

The remaining sections are devoted to the convergence analysis of RMISO. In
Section 6 we state and prove some preliminary lemmas used in the proofs of
both Theorems 3.1 and 3.4. Section 7 is devoted to convergence rate analysis
and the proof of Theorem 3.1. We then prove Theorem 3.4 in Section 8. Finally,
Appendix A gives some auxiliary lemmas, and examples of surrogate functions
are provided in Appendix B.

###  1.4. Notation

In this paper, we let $\mathbb{R}^{p}$ denote the ambient space for the
parameter space $\boldsymbol{\Theta}$ equipped with the standard inner product
$\langle\cdot,\cdot\rangle$ and the induced Euclidean norm
$\lVert\cdot\rVert$. For $\boldsymbol{\theta}\in\mathbb{R}^{p}$ and
$\varepsilon>0$, we let $B_{\varepsilon}(\boldsymbol{\theta})$ represent the
closed Euclidean ball of radius $\varepsilon$ centered at
$\boldsymbol{\theta}$. We let $\mathbbm{1}(A)$ be the indicator function of an
event $A$ which takes value $1$ on $A$ and $0$ on $A^{c}$. For a probability
distribution $\mu$ on a finite set $\mathcal{X}$, we let
$\mu_{\text{min}}=\min_{x\in\mathcal{X}}\mu(x)$. We let $a\wedge
b=\min\\{a,b\\}$ for real numbers $a$ and $b$. For a set $\mathcal{X}$ we let
$|\mathcal{X}|$ denote its cardinality. For a random variable $X$ defined on a
probability space $(\Omega,\mathcal{F},\mathbb{P})$ we let $\lVert
X\rVert_{\infty}$ be its $L_{\infty}$ norm on $\Omega$, i.e., $\lVert
X\rVert_{\infty}=\inf\\{t>0:\mathbb{P}(|X|>t)=0\\}$.

##  2\. Algorithm Statement

In this section, we state the two main algorithms used to solve (1). To do
this we start by defining first-order surrogate functions and then define a
few random variables that will be important in both implementation and
analysis. First-order surrogates are defined by

######  Definition 2.1 (First-order surrogates).

A convex function $g:\mathbb{R}^{p}\to\mathbb{R}$ is a first-order surrogate
function of $f$ at $\boldsymbol{\theta}$ if

  * (i)

$g(\boldsymbol{\theta}^{\prime})\geq f(\boldsymbol{\theta}^{\prime})$ holds
for all $\boldsymbol{\theta}^{\prime}\in\boldsymbol{\Theta}$

  * (ii)

the approximation error $h:=g-f$ is differentiable and $\nabla h$ is
$L$-Lipschitz continuous for some $L>0$; moreover $h(\boldsymbol{\theta})=0$
and $\nabla h(\boldsymbol{\theta})=0$.

We denote by $\mathcal{S}_{L}(f,\boldsymbol{\theta})$ the set of all first
order surrogates of $f$ at $\boldsymbol{\theta}$ such that $\nabla h$ is
$L$-Lipschitz.

Algorithms 1 and 2 assume that data $v\in\mathcal{V}$ are sampled according to
an algorithm which generates a stochastic process $(v_{n})_{n\geq 1}$ on state
space $\mathcal{V}$. Certain properties of this process are crucial in our
analysis, especially in proving Lemmas 7.1, 8.2, and 8.4. Below we define the
_first return time_ and _last passage time_.

######  Definition 2.2 (First return time).

For $n\geq 0$ and $v\in\mathcal{V}$, the time to return to data $v$ starting
from time $n$ is defined as

$\displaystyle\tau_{n,v}=\inf\\{j\geq 1:v_{n+j}=v\\}.$ (3)

That is, $\tau_{n,v}$ is the amount of time which one has to wait after time
$n$ for the process to return to $v$. The first return time may be viewed as a
generalization of the return times of a Markov chain. Indeed,
$\tau_{0,v}=\inf\\{n\geq 1:v_{n}=v\\}$ agrees with the classical notion of
return time from the Markov chain literature. This is closely related to the
last passage of time defined below.

######  Definition 2.3 (Last passage time).

For $n\geq 1$ and $v\in\mathcal{V}$ we define the last passage time of $v$
before time $n$ as

$\displaystyle k^{v}(n)=\sup\\{j\leq n:v_{j}=v\\}.$ (4)

If the process has not yet visited $v$, i.e. $\\{j\leq
n:v_{j}=v\\}=\emptyset$, then we set $k^{v}(n)=1$.

Both algorithms (Algorithms 1 and 2) presented in this paper proceed as
follows: 1) Sample a data point $v_{n}\in\mathcal{V}$ using an a priori
sampling algorithm; 2) choose a new surrogate $g_{n}^{v_{n}}$ of the objective
function $f^{v_{n}}$ at $\boldsymbol{\theta}_{n-1}$ and set
$g_{n}^{v}=g_{n-1}^{v}$ for $v\neq v_{n}$; 3) update the average surrogate
$\bar{g}_{n}=\sum_{v\in\mathcal{V}}g_{n}^{v}\pi(v)$; 4) output a minimizer of
the regularized average surrogate
$\bar{g}_{n}(\boldsymbol{\theta})+\Psi_{n}(\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert)$.
Therefore, the last passage time $k^{v}(n)$ appears naturally as it is the
last time the surrogate for data point $v$ has been updated during the
execution of either Algorithm 1 or 2. Thus, $g_{n}^{v}$ is a surrogate of
$f^{v}$ at $\boldsymbol{\theta}_{k^{v}(n)-1}$ and the corresponding surrogate
error at this point $h_{n}^{v}(\boldsymbol{\theta}_{k^{v}(n)-1})$ and its
gradient are equal to zero. We will use this fact crucially in the proof of
the key lemma, Lemma 7.1. Note also that $g_{n}^{v}$ is a surrogate of $f^{v}$
at $\boldsymbol{\theta}_{0}$ if $v$ has not yet been visited, which is why we
set $k^{v}(n)=1$ on this event.

Next, we state our main algorithm with two variants depending on the type of
regularization. Both are extensions of the original MISO algorithm first
proposed in [36]. The use of regularization has the effect of ensuring that
the next iterate $\boldsymbol{\theta}_{n}$ is not too far from the previous
$\boldsymbol{\theta}_{n-1}$. To handle dependent data, many algorithms use
some form of projection or regularization to achieve this property [34, 7,
49]. This allows one to control the bias introduced by dependent sampling
schemes as well as use the broader class of convex surrogates instead of
requiring them to be strongly convex.

Algorithm 1 Incremental Majorization Minimization with Dynamic Proximal
Regularization

1:Input: $\boldsymbol{\theta}_{0}\in\boldsymbol{\Theta}$ (initial estimate);
$N$ (number of iterations); $\rho>0$ (regularization parameter)

2:Option:  Regularization Schedule $\in$ {Dynamic, Constant}

3: Initialize surrogates
$g_{0}^{v}\in\mathcal{S}_{L}(f^{v},\boldsymbol{\theta}_{0})$.

4: for $n=1,\ldots,N$ do

5: sample a data point $v_{n}$

6: choose a surrogate function
$g_{n}^{v_{n}}\in\mathcal{S}_{L}(f^{v_{n}},\boldsymbol{\theta}_{n-1})$;
$g_{n}^{v}=g_{n-1}^{v}$ for all $v\neq v_{n}$

7: $\bar{g}_{n}\leftarrow\sum_{v\in\mathcal{V}}g_{n}^{v}\pi(v)$

8:
$\displaystyle\rho_{n}\leftarrow\begin{cases}\rho+\max_{v\in\mathcal{V}}(n-k^{v}(n))&\textup{
if {Regularization Schedule} $=$ {Dynamic}}\\\ \rho&\textup{if {Regularization
Schedule} $=$ {Constant}}\end{cases}$ (5)

9: update the current estimate $\boldsymbol{\theta}_{n}$:
$\displaystyle\boldsymbol{\theta}_{n}\in\operatorname*{arg\,min}_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\Bigg{[}\bar{g}_{n}(\boldsymbol{\theta})+\frac{\rho_{n}}{2}\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert^{2}\Bigg{]}$
(6)

10:output: $\boldsymbol{\theta}_{N}$

Algorithm 2 Incremental Majorization Minimization with Diminishing Radius

1:Input: $\boldsymbol{\theta}_{0}\in\boldsymbol{\Theta}$ (initial estimate);
$N$(number of iterations); $(r_{n})_{n\geq 1}$(diminishing search radius
constraints)

2: Initialize surrogates
$g_{0}^{v}\in\mathcal{S}_{L}(f^{v},\boldsymbol{\theta}_{0})$.

3: for $n=1,\ldots,N$ do

4: sample a data point $v_{n}$

5: choose a surrogate function
$g_{n}^{v_{n}}\in\mathcal{S}_{L}(f^{v_{n}},\boldsymbol{\theta}_{n-1})$;
$g_{n}^{v}=g_{n-1}^{v}$ for all $v\neq v_{n}$

6: $\bar{g}_{n}\leftarrow\sum_{v\in\mathcal{V}}g_{n}^{v}\pi(v)$

7: update the current estimate $\boldsymbol{\theta}_{n}$:
$\displaystyle\boldsymbol{\theta}_{n}\in\operatorname*{arg\,min}_{\boldsymbol{\theta}\in\boldsymbol{\Theta}\cap
B_{r_{n}}(\boldsymbol{\theta}_{n-1})}\bar{g}_{n}(\boldsymbol{\theta})$ (7)

8:output: $\boldsymbol{\theta}_{N}$

Algorithm 1 uses _Proximal Regularization_ (PR) while Algorithm 2 utilizes a
_Diminishing Radius_ (DR) restriction. The use of DR was first introduced in
[28] to analyze the convergence and complexity of cyclic block majorization-
minimization algorithms for minimizing smooth functions on convex constraint
sets. It may be viewed as a trust-region method where the trust region takes
the form of an Euclidean ball with a diminishing radius $r_{n}$. In contrast
PR penalizes the iterate $\boldsymbol{\theta}_{n}$’s distance from
$\boldsymbol{\theta}_{n-1}$ through the addition of the quadratic term
$\frac{\rho_{n}}{2}\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert^{2}$
to the averaged surrogate. One key difference between PR and DR is that the
use of DR does not change the gradient of the objective while PR does.

Proximal regularization is widely used in the literature [14, 24, 47, 1]. The
multiplicative parameter $\rho_{n}$ may be viewed similarly to the inverse of
the step-size in standard gradient descent. Hence, on the one hand, the use of
a constant regularization parameter may be compared with gradient descent
using a constant step size. On the other hand, the dynamic regularization
increasingly penalizes large values of
$\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert$ as the maximum
amount of time any node $v$ has been left un-visited by the sampling algorithm
(known as the ‘cover time’ in the Markov chain literature) increases. It
therefore adapts the proximal regularization to the stochastic process
$(v_{n})$ helping to ensure that the average surrogate $\bar{g}_{n}$ remains
an accurate approximation of the objective function $f$ near the most recent
iterate.

##  3\. Main Results

###  3.1. Optimality Conditions

In this subsection, we introduce the optimality conditions used in this paper
and related quantities. Here we denote $f$ to be a general objective function
$f:\boldsymbol{\Theta}\to\mathbb{R}$, but elsewhere $f$ will refer to the
objective function in (1) unless otherwise stated.

For a given function $f$ and
$\boldsymbol{\theta}^{*},\boldsymbol{\theta}\in\boldsymbol{\Theta}$, we define
its _directional derivative_ at $\boldsymbol{\theta}^{*}$ in the direction
$\boldsymbol{\theta}-\boldsymbol{\theta}^{*}$ as

$\displaystyle\nabla
f(\boldsymbol{\theta}^{*},\boldsymbol{\theta}-\boldsymbol{\theta}^{*}):=\lim_{\alpha\to
0^{+}}\frac{f(\boldsymbol{\theta}^{*}+\alpha(\boldsymbol{\theta}-\boldsymbol{\theta}^{*}))-f(\boldsymbol{\theta}^{*})}{\alpha}$
(8)

A necessary first order condition for $\boldsymbol{\theta}^{*}$ to be a local
minimum of $f$ is to require $\nabla
f(\boldsymbol{\theta}^{*},\boldsymbol{\theta}-\boldsymbol{\theta}^{*})\geq 0$
for all $\boldsymbol{\theta}\in\boldsymbol{\Theta}$ (see [36]). Thus we say
that $\boldsymbol{\theta}^{*}\in\boldsymbol{\Theta}$ is a _stationary point_
of $f$ over $\boldsymbol{\Theta}$ if

$\displaystyle\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}^{*}\rVert\leq
1}-\nabla
f(\boldsymbol{\theta}^{*},\boldsymbol{\theta}-\boldsymbol{\theta}^{*})\leq 0.$
(9)

If $f$ is differentiable, then $\nabla
f(\boldsymbol{\theta}^{*},\boldsymbol{\theta}-\boldsymbol{\theta}^{*})=\langle\nabla
f(\boldsymbol{\theta}^{*}),\boldsymbol{\theta}-\boldsymbol{\theta}^{*}\rangle$
so the above requirement becomes

$\displaystyle\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}^{*}\rVert\leq
1}\langle-\nabla
f(\boldsymbol{\theta}^{*}),\boldsymbol{\theta}-\boldsymbol{\theta}^{*}\rangle\leq
0.$ (10)

If $\boldsymbol{\Theta}$ is convex, this is equivalent to $-\nabla
f(\boldsymbol{\theta}^{*})$ being in the normal cone of $\boldsymbol{\Theta}$
at $\boldsymbol{\theta}^{*}$. If $\boldsymbol{\theta}^{*}$ is in the interior
of $\boldsymbol{\Theta}$ then it implies that $\lVert\nabla
f(\boldsymbol{\theta}^{*})\rVert=0$.

For iterative algorithms, the first order optimality condition (9) may hardly
be satisfied in a finite number of iterations. A practically important
question is how the worst-case number of iterations required to achieve an
$\varepsilon$-approximate solution scales with the desired precision
$\varepsilon$. We say that $\boldsymbol{\theta}^{*}\in\boldsymbol{\Theta}$ is
an $\varepsilon$-_approximate stationary point_ of $f$ over
$\boldsymbol{\Theta}$ if

$\displaystyle\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}^{*}\rVert\leq
1}-\nabla
f(\boldsymbol{\theta}^{*},\boldsymbol{\theta}-\boldsymbol{\theta}^{*})\leq\varepsilon.$
(11)

This notion of $\varepsilon$-approximate solution is consistent with the
corresponding notion for unconstrained problems. In fact, if $f$ is
differentiable, and if $\boldsymbol{\theta}^{*}$ is distance at least one away
from the boundary $\partial\boldsymbol{\Theta}$, then (11) reduces to
$\lVert\nabla f(\boldsymbol{\theta}^{*})\rVert\leq\varepsilon$. For each
$\varepsilon>0$, we then define the _worst-case iteration complexity_ of an
algorithm for solving (1) as

$\displaystyle N_{\varepsilon}(\boldsymbol{\theta}_{0}):=\inf\\{n\geq
1:\boldsymbol{\theta}_{n}\text{ is an $\varepsilon$-approximate stationary
point of $f$ over $\boldsymbol{\Theta}$}\\},$ (12)

where $(\boldsymbol{\theta}_{n})_{n\geq 0}$ is a sequence of iterates produced
by the algorithm with initial estimate $\boldsymbol{\theta}_{0}$.

###  3.2. Assumptions

In this subsection, we state our assumptions for establishing the main
results. Throughout this paper, we denote by $\mathcal{F}_{n}$ the
$\sigma$-algebra generated by the samples $v_{1},\ldots,v_{n}$ and the
parameters $\boldsymbol{\theta}_{0},\ldots,\boldsymbol{\theta}_{n}$ produced
by Algorithm 1 or 2. With this definition, $(\mathcal{F}_{n})_{n\geq 1}$
defines a filtration.

We first state our main assumption on the sampling scheme.

######  A1Recurrent data sampling.

The sequence $(v_{n})_{n\geq 1}$ of data samples defines a stochastic process
which satisfies the following property: for each $v\in\mathcal{V}$,
$\sup_{n\geq
1}\left\lVert\mathbb{E}[\tau_{n,v}|\mathcal{F}_{n}]\right\rVert_{\infty}<\infty$,
i.e., the expected return time conditioned on $\mathcal{F}_{n}$ is uniformly
bounded.

Assumption A1 states that the data $(v_{n})_{n\geq 1}$ are sampled in such a
way that the expected time between visits to a particular data point is finite
and uniformly bounded. Generalizing the notion of positive recurrence in
Markov chain theory, we say a sampling algorithm is _recurrent_ if A1 is
satisfied. We emphasize that recurrence is the only requirement we make of the
sampling process in order to prove the convergence rate guarantees in Theorem
3.1 and the asymptotic convergence in Theorem 3.4 (We do not assume
independence or Markovian dependence, etc.). We include below a list of some
commonly used recurrent sampling algorithms.

1.

(i.i.d. sampling) Sampling data i.i.d from a fixed distribution is the most
common assumption in the literature [35, 36, 3, 53, 23]. Suppose we sample
$v_{n}$ i.i.d from some distribution $\mu$ on $\mathcal{V}$. Then the
$\tau_{n,v}$ are independent geometric random variables taking values from
$\\{1,2,\dots,\\}$ with success probability $\mu(v)$, so
$\mathbb{E}[\tau_{n,v}|\mathcal{F}_{n}]=1/\mu(v)$. In particular, if $\mu$ is
the uniform distribution
$\mathbb{E}[\tau_{n,v}|\mathcal{F}_{n}]=|\mathcal{V}|$ for all $n$ and $v$. We
point out that it is not necessary for our analysis that $\mu$ agree with the
distribution $\pi$ defining $f$ in (1). Such flexibility may help when using
importance sampling as has been studied in [46, 58, 52].

2.

(MCMC) Markov chain Monte Carlo methods (see e.g. [31, Ch. 3]) produce a
Markov chain $(v_{n})$ on $\mathcal{V}$. If this chain is irreducible then
$\max_{v,w}\mathbb{E}_{w}[\tau_{v}^{+}]$ is finite [31, Lem 1.13] where
$\mathbb{E}_{w}[\cdot]$ is the expectation taken with respect to initial state
$w$ and $\tau_{v}^{+}=\inf\\{n\geq 1:v_{n}=v\\}$ is the first return time to
$v$. For any $n$, $v$, and initial distribution $\nu$, the Markov property
implies
$\mathbb{E}_{\nu}[\tau_{n,v}|\mathcal{F}_{n}]=\mathbb{E}_{v_{n}}[\tau_{0,v}]=\mathbb{E}_{v_{n}}[\tau_{v}^{+}]$.
So any irreducible Markov chain satisfies A1. In contrast to a number of other
works studying stochastic optimization under Markovian sampling [34, 30, 55,
2, 8] we do _not_ require the stationary distribution of this chain to agree
with $\pi$.

3.

(Cyclic sampling) In cyclic sampling one samples data in order according to
some enumeration until the dataset is exhausted. This process is then repeated
until convergence. The authors of [29] show that iteration complexity for SGD
can be improved from $O(\varepsilon^{-4})$ to $O(\varepsilon^{-3})$ using such
methods. To see that A1 holds in this setting, we simply notice that
$\tau_{n,v}\leq|\mathcal{V}|$ for all $n$ and $v$.

4.

(Reshuffling) Random reshuffling is similar to cycling sampling except that
the dataset is randomly permuted at the beginning of each epoch [29]. It was
observed empirically that random reshuffling performs better than i.i.d.
sampling in [6] and further studied in [16, 29]. The authors of [29] show the
same improvement in iteration complexity for SGD as for cyclic sampling. In
this case, A1 is satisfied since $\tau_{n,v}\leq 2|\mathcal{V}|$.

In Markov chain theory, the quantity $\max_{v,w}\mathbb{E}_{w}[\tau_{v}^{+}]$
is commonly denoted $t_{\text{hit}}$. Adapting this to our general setting we
define

$\displaystyle t_{\text{hit}}:=\max_{v\in\mathcal{V}}\,\sup_{n\geq
1}\,\,\left\lVert\mathbb{E}[\tau_{n,v}|\mathcal{F}_{n}]\right\rVert_{\infty}$
(13)

for each $v$ when A1 holds. Continuing the connection with Markov chains, we
also define

$\displaystyle t_{\odot}:=\sup_{n\geq
1}\,\,\left\lVert\sum_{v\in\mathcal{V}}\mathbb{E}[\tau_{n,v}|\mathcal{F}_{n}]\pi(v)\right\rVert_{\infty}$
(14)

where $\pi$ is as in (1). Note that if $v_{n}$ is an irreducible Markov chain
then the Markov property implies

$\displaystyle
t_{\odot}=\max_{w\in\mathcal{V}}\sum_{v\in\mathcal{V}}\mathbb{E}_{w}[\tau_{v}^{+}]\pi(v).$
(15)

This is closely related to the _target time_ defined by

$\displaystyle t_{\odot}^{w}=\sum_{v\in V}\mathbb{E}_{w}[\tau_{v}]\pi(v)$ (16)

with the difference being that here $\tau_{v}=\inf\\{n\geq 0:v_{n}=v\\}$ is
the first hitting time of $v$ rather than the first return time to $v$. The
_random target lemma_ [31, Lem. 10.1] states that if $\pi$ is the stationary
distribution of the Markov chain, then the target time is independent of the
starting state $w$. In this case, the quantities (15) and (16) only differ by
one. This can be seen by first noting that
$\mathbb{E}_{w}[\tau_{v}]=\mathbb{E}_{w}[\tau_{v}^{+}]$ if $v\neq w$. This
leaves only the difference
$\mathbb{E}_{w}[\tau_{w}^{+}]\pi(w)-\mathbb{E}_{w}[\tau_{w}]\pi(w)$. The
second term is equal to zero and the first equals one since if $\pi$ is the
unique stationary distribution for the chain,
$\pi(w)=\frac{1}{\mathbb{E}_{w}[\tau_{w}^{+}]}$ [31, Prop 1.19].

For transitive irreducible Markov chains, $t_{\text{hit}}$ and $t_{\odot}$ are
comparable [31, Cor. 10.3]. However, in other situations $t_{\odot}$ may be
much smaller that $t_{\text{hit}}$. For instance, consider the case that
sampling is done by a simple random walk on a connected graph and $\pi$ is the
stationary distribution. If there is an outlier node with a much lower degree
than the remaining nodes we may have $t_{\odot}\ll t_{\text{hit}}$. An example
of such a graph is displayed in Figure 1 which we refer to as the "lonely"
graph. In the lonely graph, $|\mathcal{V}|-1$ vertices form a clique and the
remaining vertex has degree one. A random walk on this graph has worst-case
hitting time $t_{\text{hit}}=O(|\mathcal{V}|^{2})$: its value is tied to the
lonely vertex with degree one which has low probability of being visited. On
the other hand, $t_{\odot}$ is only $O(|\mathcal{V}|)$ since it only depends
on an _average_ hitting time instead of the worst case. See [31, Ex. 10.4] for
more details.

Figure 1. Lonely graph

We next list some standard requirements on the objectives $f^{v}$.

######  A2Lower-bounded objective and directional derivatives.

For all $v\in\mathcal{V}$ and
$\boldsymbol{\theta},\boldsymbol{\theta}^{\prime}\in\boldsymbol{\Theta}$, the
function $f^{v}$ is bounded below, i.e.
$\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}f^{v}(\boldsymbol{\theta})>-\infty$.
Moreover, the directional derivative

$\displaystyle\nabla
f^{v}(\boldsymbol{\theta},\boldsymbol{\theta}^{\prime}-\boldsymbol{\theta}):=\lim_{\alpha\to
0^{+}}\frac{f^{v}(\boldsymbol{\theta}+\alpha(\boldsymbol{\theta}^{\prime}-\boldsymbol{\theta}))-f^{v}(\boldsymbol{\theta})}{\alpha}$
(17)

exists.

The boundedness assumption in A2 implies that
$\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}f(\boldsymbol{\theta})\geq
M>-\infty$ for some $M$. For the remainder of this paper we will denote
$\Delta_{0}:=\bar{g}_{0}(\boldsymbol{\theta}_{0})-\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}f(\boldsymbol{\theta})$.
It is important to note that if the initial surrogates $g_{0}^{v}$ are in
$\mathcal{S}_{L}(f^{v},\boldsymbol{\theta}_{0})$ (as is the case in both
Algorithms 1 and 2) then
$g_{0}(\boldsymbol{\theta}_{0})=f(\boldsymbol{\theta}_{0})$ so
$\Delta_{0}=f(\boldsymbol{\theta}_{0})-\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}f(\boldsymbol{\theta})$.
The regularity assumption in A2 was used in [36] and is necessary for
analyzing our algorithms using the approximate stationarity measure (11). For
Algorithm 2 we make the following stronger but common assumption which is
crucial to our analysis:

######  A2’.

For each $v\in\mathcal{V}$, the function $f^{v}$ is continuously
differentiable and $\nabla f^{v}$ is $L$-Lipschitz continuous.

For simplicity, we assume that if A2’ holds then the Lipschitz constant of
$f^{v}$ agrees with that of the corresponding approximation error $h^{v}$.

Finally, A3 states that the radii in Algorithm 2 decrease slowly, but not too
slowly. This is analogous to the square-summability of step sizes in gradient
descent.

######  A3Square-summable and non-summable radii.

The sequence $(r_{n})_{n\geq 1}$ is non-increasing,
$\sum_{n=1}^{\infty}r_{n}=\infty$, and $\sum_{n=1}^{\infty}r_{n}^{2}<\infty$.

###  3.3. Statement of main results

In this section and in the remainder of the paper, we let
$h_{n}^{v}(\boldsymbol{\theta}):=g_{n}^{v}(\boldsymbol{\theta})-f^{v}(\boldsymbol{\theta})$
be the surrogate error at $v$ and
$\bar{h}_{n}(\boldsymbol{\theta})=\sum_{v\in\mathcal{V}}h_{n}^{v}(\boldsymbol{\theta})\pi(v)$.
Throughout, let $(\boldsymbol{\theta}_{n})_{n\geq 0}$ be an output of either
Algorithm 1 or 2 with arbitrary initialization $\boldsymbol{\theta}_{0}$ and
assume $g_{n}^{v_{n}}\in\mathcal{S}_{L}(f^{v},\boldsymbol{\theta}_{n-1})$ for
$n\geq 1$ and some $L>0$. We recall the definitions of the following
constants: $t_{\text{hit}}$ and $t_{\odot}$ are as in (13) and (14)
respectively, $\rho$ is the user-defined proximal regularization parameter
when Algorithm 1 is run with Regularization Schedule $=$ Constant and the
deterministic part of the regularization parameter when Algorithm 1 is run
with Regularization Schedule $=$ Dynamic,
$\Delta_{0}=g_{0}(\boldsymbol{\theta}_{0})-\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}f(\boldsymbol{\theta})$
is the initial difference between the value of the averaged surrogate and the
optimal value of $f$. We will consider one of the following cases:

######  C1.

Assume A1-A2. Use Algorithm 1 with Regularization Schedule $=$ Constant

######  C2.

Assume A1 \- A2. Use Algorithm 1 with Regularization Schedule $=$ Dynamic.

######  C3.

Assume A1 \- A3 and A2’. Use Algorithm 2.

Our first main result, Theorem 3.1, gives worst-case upper bounds on the
expected rate of convergence to optimality. For each of the cases C1-C3 we
give rates of convergence for the objective function $f$.

######  Theorem 3.1 (Rate of Convergence to Stationarity).

Algorithms 1 and 2 satisfy the following:

  1. (i)

Assume case C1. Then

$\displaystyle\min_{1\leq n\leq
N}\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}-\nabla
f(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\right]\leq\frac{\sqrt{2\Delta_{0}}\left(\sqrt{\rho}+\frac{Lt_{\odot}}{\sqrt{\rho}}\right)}{\sqrt{N}}.$
(18)

In particular, the optimal choice $\rho=Lt_{\odot}$ results in

$\displaystyle\min_{1\leq n\leq
N}\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}-\nabla
f(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\right]\leq
2\sqrt{\frac{2\Delta_{0}Lt_{\odot}}{N}}.$ (19)

  2. (ii)

Assume case C2. Then

$\displaystyle\min_{1\leq n\leq
N}\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}-\nabla
f(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\right]\leq\frac{\sqrt{2\Delta_{0}}\left(\sqrt{\rho+(2t_{\text{hit}}+1)\log_{2}(4|\mathcal{V}|)}+\frac{Lt_{\odot}}{\sqrt{\rho}}\right)}{\sqrt{N}}.$
(20)

Taking $\rho=Lt_{\odot}$ gives

$\displaystyle\min_{1\leq n\leq
N}\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\leq
1\rVert}-\nabla
f(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\right]\leq\frac{\sqrt{2\Delta_{0}}\left(\sqrt{Lt_{\odot}+(2t_{\text{hit}}+1)\log_{2}(4|\mathcal{V}|)}+\sqrt{Lt_{\odot}}\right)}{\sqrt{N}}.$
(21)

  3. (iii)

Assume case C3. Then

$\displaystyle\min_{1\leq n\leq
N}\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}\langle-\nabla
f(\boldsymbol{\theta}_{n}),\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rangle\right]\leq\frac{\Delta_{0}+\sqrt{\frac{2L}{\pi_{\text{min}}}C_{N}\Delta_{0}}+\Big{(}3+t_{\odot}\Big{)}C_{N}L}{\sum_{n=1}^{N}1\wedge
r_{n+1}},$ (22)

where $C_{N}=\sum_{n=1}^{N}r_{n}^{2}$.

To our knowledge, the rates of convergence given in Theorem 3.1 are entirely
new for first-order algorithms with general recurrent data sampling. In
contrast to the convergence result for SAG in [15] that depends on the hitting
time $t_{\text{hit}}$, Algorithm 1 with constant proximal regularization (case
C1) depends on the possibly much smaller target time $t_{\odot}$.

The dependence on $t_{\odot}$ (as well as $t_{\text{hit}}$ for case C2)
suggest that one can accelerate convergence by choosing an algorithm with the
smallest values of these constants appropriate for the context. If one is free
to choose any sampling scheme, then the optimal order of dependence on
$|\mathcal{V}|$ for $t_{\odot}$ is achieved by sampling the data i.i.d from
$\pi$. Indeed, in this case
$\mathbb{E}[\tau_{n,v}|\mathcal{F}_{n}]=\frac{1}{\pi(v)}$ for all $n$ and $v$
so $t_{\odot}=|\mathcal{V}|$. In contrast,
$t_{\text{hit}}=\frac{1}{\pi_{\text{min}}}$ for i.i.d. sampling which may be
much larger than $t_{\odot}$ if $\pi$ is not uniform. In settings such as
decentralized optimization i.i.d sampling may not be feasible [21, 22].
However, one could opt for cyclic sampling which has
$t_{\text{hit}}=|\mathcal{V}|$ for any distribution $\pi$. It is always true
that $t_{\odot}\leq t_{\text{hit}}$, but the two are comparable with cyclic
sampling. Indeed if $n$ is fixed then for each $1\leq k\leq|\mathcal{V}|$
there is a $v$ with $\mathbb{E}[\tau_{n,v}|\mathcal{F}_{n}]=k$. So $t_{\odot}$
is the largest possible value of

$\displaystyle\sum_{v\in\mathcal{V}}\sigma_{v}\pi(v)$ (23)

where $\sigma$ ranges over all permutations of
$\\{1,2,\ldots,|\mathcal{V}|\\}$. In particular,
$t_{\odot}=\frac{|\mathcal{V}|-1}{2}$ if $\pi$ is uniform.

We also remark that it would not be difficult to adapt our analysis to the
specific setting that the return times $\tau_{n,v}$ are uniformly bounded. In
this case the results of Theorem 3.1 would hold _almost surely_ , meaning the
above inequalities would remain true without the expectation on the left-hand
side. In this case $t_{\odot}$ would be replaced with

$\displaystyle\sup_{n\geq 1}\left\lVert\sum_{v\in
V}\tau_{n,v}\pi(v)\right\rVert_{\infty}$ (24)

and $t_{\text{hit}}$ replaced by $\max_{v\in\mathcal{V}}\sup_{n\geq
1}\lVert\tau_{n,v}\rVert_{\infty}$. This is true for cyclic sampling and
random reshuffling but not for i.i.d. or general Markovian sampling.

Next, we state a corollary of Theorem 3.1 which specializes these results to
the setting of unconstrained nonconvex optimization.

######  Corollary 3.2.

Assume either $\boldsymbol{\Theta}=\mathbb{R}^{p}$ or that there exists
$c\in(0,1]$ so that
$\text{dist}(\boldsymbol{\theta}_{n},\partial\boldsymbol{\Theta})\geq c$ for
all $n\geq 1$.

  1. (i)

Let $(\boldsymbol{\theta}_{n})_{n\geq 0}$ be an output of Algorithm 1. Assume
case C1 or C2. Then for any $N\geq 1$

$\displaystyle\min_{1\leq n\leq N}\mathbb{E}\left[\lVert\nabla
f(\boldsymbol{\theta}_{n})\rVert\right]=O\left(N^{-1/2}\right).$ (25)

  2. (ii)

Let $(\boldsymbol{\theta}_{n})_{n\geq 0}$ be an output of Algorithm 2. Assume
case C3. Then for any $N\geq 1$

$\displaystyle\min_{1\leq n\leq N}\mathbb{E}\left[\lVert\nabla
f(\boldsymbol{\theta}_{n})\rVert\right]=O\left(\left(\sum_{n=1}^{N}1\wedge
r_{n}\right)^{-1}\right).$ (26)

Notice that we may take $r_{n}=\frac{1}{\sqrt{n}\log n}$ in Algorithm 2. Then
Corollary 3.2 implies

$\displaystyle\min_{1\leq n\leq N}\mathbb{E}\left[\lVert\nabla
f(\boldsymbol{\theta}_{n})\rVert\right]=O\left(\frac{\log N}{\sqrt{N}}\right)$
(27)

holds for case C3.

Finally, Corollary 3.3 states the iteration complexity of Algorithms 1 and 2.

######  Corollary 3.3 (Iteration Complexity).

Algorithms 1 and 2 have the following worst case iteration complexity:

  1. (i)

Let $(\boldsymbol{\theta}_{n})_{n\geq 0}$ be an output of Algorithm 1. Assume
case C1 or C2. Then
$N_{\varepsilon}(\boldsymbol{\theta}_{0})=O(\varepsilon^{-2})$.

  2. (ii)

Let $(\boldsymbol{\theta}_{n})_{n\geq 0}$ be an output of Algorithm 2. Assume
case C3. Then
$N_{\varepsilon}(\boldsymbol{\theta}_{0})=O(\varepsilon^{-2}(\log\varepsilon^{-1})^{2})$.

As summarized in Table 1, the iteration complexity of our method matches the
iteration complexity of variance-reduced algorithms such as SAG [53] and its
extension to Markovian sampling [15] exactly for Algorithm 1 and up to
logarithmic factors for Algorithm 2. This may be compared with
$\tilde{O}(\varepsilon^{-4})$ for (P)SGD and its variants in the dependent
setting. However, our method requires the trade-off, shared with SAG, that one
must store information for each of the surrogate functions $g_{n}^{v}$ which
generally means storing past gradients. Thus, the memory requirement of our
algorithm scales linearly with $|\mathcal{V}|$. Our algorithm is also only
applicable to expected loss minimization for problems with finite sample
space.

|  Iteration complexity | Memory | Sampling | Sampling dependence
---|---|---|---|---
AdaGrad [2] | $\tilde{O}(\varepsilon^{-4})$ | $O(1)$ | Markovian | $O(\sqrt{t_{\text{mix}}})$
SGD [55, 15] | $\tilde{O}(\varepsilon^{-4})$ | $O(1)$ | Markovian | $O(\sqrt{t_{\text{mix}}})$
SGD [38, 29] | $O(\varepsilon^{-3})$ | $O(1)$ | Reshuffling | $O(\sqrt{|\mathcal{V}|})$
SAG [15] | $O(\varepsilon^{-2})$ | $O(|\mathcal{V}|)$ | Markovian | $O(\sqrt{t_{\text{hit}}})$
RMISO Case C1 | $O(\varepsilon^{-2})$ | $O(|\mathcal{V}|)$ | Recurrent | $O(\sqrt{t_{\odot}})$
RMISO Case C2 | $O(\varepsilon^{-2})$ | $O(|\mathcal{V}|)$ | Recurrent | $O(\sqrt{t_{\odot}+t_{\text{hit}}})$
RMISO Case C3 | $\tilde{O}(\varepsilon^{-2})$ | $O(|\mathcal{V}|)$ | Recurrent | $O(t_{\odot})$
Table 1. Comparison of iteration complexity for non-convex optimization with
non-i.i.d data. The notation $\tilde{O}(\cdot)$ omits logarithmic factors.
Here $t_{\text{mix}}$ represents dependence on the mixing time of the Markov
chain [31, Sec. 4.5].

Our last main result, Theorem 3.4, concerns the asymptotic behavior of
Algorithms 1 and 2. Though Algorithm 1 with constant proximal regularization
is the simplest of our proposed methods and has the best dependence on the
constants $L$ and $t_{\odot}$, it appears that stronger regularization schemes
as in cases C2 and C3 are needed for obtaining asymptotic convergence
guarantees. We discuss this in more detail in subsection 3.4 below as well as
in remark 8.3 in Section 8.

######  Theorem 3.4 (Global Convergence).

Algorithms 1 and 2 have the following asymptotic convergence properties:

  1. (i)

For case C2

$\displaystyle\mathop{\rm
lim\,sup}\limits_{n\to\infty}\,\,\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}-\nabla
f(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\right]\leq
0.$ (28)

  2. (ii)

For case C3 almost surely every limit point of
$(\boldsymbol{\theta}_{n})_{n\geq 1}$ is stationary for $f$ over
$\boldsymbol{\Theta}$.

Theorem 3.4 shows that although RMISO with diminishing radius requires
computing a projection at each step and has higher-order dependence on
$t_{\odot}$, it enjoys the strongest asymptotic guarantees. RMISO with dynamic
proximal regularization is somewhere in the middle. It has lower order
dependence on $t_{\odot}$ than the diminishing radius version but also depends
on $t_{\text{hit}}$ and has slightly higher dependence on $|\mathcal{V}|$ than
the constant proximal regularization version. However, we are able to show
that the optimality gap converges in the $L^{1}$ sense for this method. In
particular, notice that in the familiar case that
$\boldsymbol{\Theta}=\mathbb{R}^{p}$ and $f$ is differentiable (28) implies
that $\lim_{n\to\infty}\mathbb{E}[\lVert\nabla
f(\boldsymbol{\theta}_{n})\rVert]=0$.

###  3.4. Further remarks on main results

######  Remark 3.5 (Comparison with i.i.d sampling).

If the sequence $(v_{n})$ is formed by sampling vertices uniformly at random
from $\mathcal{V}$ then, as previously mentioned, the return times
$\tau_{n,v}$ are i.i.d geometric random variables with parameter
$\frac{1}{|\mathcal{V}|}$. Then
$\mathbb{E}[\tau_{n,v}|\mathcal{F}_{n}]=|\mathcal{V}|$ for each $n$ and $v$ so
$t_{\odot}=|\mathcal{V}|$. Substituting $|\mathcal{V}|$ for $t_{\odot}$ in the
optimal bound given for case C1 in Theorem 3.1 we recover the result given for
MISO in the i.i.d setting in [26] up to a factor of two. This is in spite of
the fact that our analytical approach is necessarily different to handle
general recurrent sampling and shows that our results are tight.

######  Remark 3.6 (Iterate stability and regularization).

Here we give some remarks on the use of diminishing radius and proximal
regularization in Algorithms 1 and 2.

The diminishing radius restriction in Algorithm 2 is a "hard" regularization
technique. It bakes necessary iterate stability directly into the problem by
enforcing the stronger condition
$\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert\leq r_{n}$. In
comparison with proximal regularization, diminishing radius bounds the one-
step iterate difference by a deterministic quantity. Moreover, as is argued in
the proof of Theorem 3.4 (ii) and the preceding propositions, sufficiently
often the iterate $\boldsymbol{\theta}_{n}$ obtained by minimizing the
$\bar{g}_{n}$ over the trust region in fact minimizes $\bar{g}_{n}$ over the
entire feasible set $\boldsymbol{\Theta}$. This allows us to prove that limit
points of the iterates produced by Algorithm 2 are stationary even for general
recurrent sampling schemes. However, diminishing radius introduces the
drawback of needing to compute a projection at each step of the optimization
process.

Compared to diminishing radius, proximal regularization is a form of "soft"
regularization. It is less restrictive than diminishing radius, but only
allows us to derive a weaker form of iterate stability: we only have
$\sum_{n=1}^{\infty}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}<\infty$
instead of the stronger
$\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert\leq r_{n}$
which makes some aspects of the analysis slightly more challenging.

As mentioned in Section 3.2, the use of dynamic proximal regularization in
case C2 adapts to the sampling process and increasingly penalizes large values
of $\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert$ as the
amount of time any vertex is left un-visited increases. The drawback is that
we are unable to prove almost sure asymptotic convergence in Theorem 3.4 for
case C2 as we are for case C3. If we could show
$\rho_{n}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert\to 0$,
then asymptotic convergence would follow. However, $\rho_{n}$ is not bounded
and can take arbitrarily large values, albeit with low probability. But as we
show in Lemma 8.1, $\mathbb{E}[\rho_{n}]$ is uniformly bounded, which allows
us to deduce
$\mathbb{E}[\rho_{n}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert]\to
0$ and leads to the $L^{1}$-convergence result in Theorem 3.4.

For case C1 it is relatively straightforward to show
$\rho_{n}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert\to 0$
since $\rho_{n}$ is constant. In this case, the difficulty lies in showing
$\lVert\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n})\rVert\to 0$. We _are_ able
to prove this for case C2 however because the use of dynamic proximal
regularization allows us to show that the sequence
$\max_{v\in\mathcal{V}}(n-k^{v}(n))\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert^{2}$
is summable. More detail is given in the proofs and discussion in subsection
8.1.

##  4\. Overview of convergence analysis

In this section, we give an overview of the standard pipeline for analyzing
stochastic optimization algorithms and highlight our new approach to handle
the general recurrent sampling schemes.

###  4.1. Convergence analysis using mixing times

The analysis of first-order methods such as SGD generally relies on
conditionally unbiased gradient estimates. In the context of solving problem
(1), this is to require

$\displaystyle\mathbb{E}[\nabla
f^{v_{n}}(\boldsymbol{\theta}_{n-1})|\mathcal{F}_{n-1}]=\nabla
f(\boldsymbol{\theta}_{n-1})$ (29)

where $\mathcal{F}_{n}$ is the filtration of information up to time $n$.
However, in the dependent data setting (29) does not hold which complicates
the analysis significantly. Previous works (e.g [55, 7, 45, 30, 33, 34, 2])
use a "conditioning on the distant past" argument. Specifically, for SGD, one
assumes that $(v_{n})$ is a Markov chain on state space $\mathcal{V}$ which
mixes exponentially fast to its stationary distribution $\pi$ with parameter
$\lambda$. One then considers the quantity

$\displaystyle\mathbb{E}[\nabla
f^{v_{n}}(\boldsymbol{\theta}_{n-a_{n}})|\mathcal{F}_{n-a_{n}}]$ (30)

where $a_{n}$ is a slowly growing sequence satisfying $\sum_{n\geq
1}\lambda^{a_{n}}<\infty$. By further assuming either uniformly bounded
gradients as in [55] or that the conditional expectations
$\mathbb{E}[\lVert\nabla
f^{v_{n+1}}(\boldsymbol{\theta})\rVert|\mathcal{F}_{n}]$ are uniformly bounded
as in [2], one can show that

$\displaystyle\left\lVert\nabla
f(\boldsymbol{\theta}_{n-a_{n}})-\mathbb{E}[\nabla
f^{v_{n}}(\boldsymbol{\theta}_{n-a_{n}})|\mathcal{F}_{n-a_{n}}]\right\rVert=O(\lambda^{a_{n}}).$
(31)

Using Lipschitz continuity of gradients, it is then established that

$\displaystyle\lVert\nabla f(\boldsymbol{\theta}_{n})-\mathbb{E}[\nabla
f^{v_{n}}(\boldsymbol{\theta}_{n})|\mathcal{F}_{n-a_{n}}]\rVert=O(\lambda^{a_{n}})+O(\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-a_{n}}\rVert),$
(32)

allowing one to control the bias in the stochastic gradient estimate.
Conventional analysis may then be used to prove

$\displaystyle\sum_{n=1}^{\infty}\gamma_{n}\mathbb{E}\left[\langle\nabla
f(\boldsymbol{\theta}_{n}),\nabla
f^{v_{n}}(\boldsymbol{\theta}_{n})\rangle\right]<\infty,$ (33)

where $\gamma_{n}$ denotes the stepsize at iteration $n$. Combining this with
(32), it can finally be deduced that $\sum_{n\geq
1}\gamma_{n}\mathbb{E}[\lVert\nabla
f(\boldsymbol{\theta}_{n})\rVert^{2}]<\infty$ for an appropriate stepsize
$\gamma_{n}$ which gives the desired convergence rate.

This technique was further developed to analyze convergence rates of
_stochastic majorization minimization_ (SMM) [35] algorithms under Markovian
sampling in [34]. In the context of problem (1), SMM proceeds by minimizing a
recursively defined majorizing surrogate of the empirical loss function. The
algorithm is stated concisely as follows:

$\displaystyle(\textbf{SMM}):\begin{cases}&\text{Sample }v_{n}\text{ from the
conditional distribution }\pi(\cdot|\mathcal{F}_{n-1})\\\
&g_{n}\leftarrow\text{Strongly convex majorizing surrogate of $f^{v_{n}}$}\\\
&\boldsymbol{\theta}_{n}\in\operatorname*{arg\,min}_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left(\bar{g}_{n}(\boldsymbol{\theta}):=(1-w_{n})\bar{g}_{n-1}(\boldsymbol{\theta})+w_{n}g_{n}(\boldsymbol{\theta})\right)\end{cases}$
(34)

where $(w_{n})_{n\geq 1}$ is a non-increasing sequence of weights. A crucial
step in the analysis of SMM is to show that

$\displaystyle\sum_{n=1}^{\infty}w_{n}\mathbb{E}[|\bar{g}_{n}(\boldsymbol{\theta}_{n})-\bar{f}_{n}(\boldsymbol{\theta}_{n})|]<\infty$
(35)

where $\bar{f}_{n}$ is the empirical loss function satisfying the recursion
$\bar{f}_{n}(\boldsymbol{\theta}):=(1-w_{n})\bar{f}_{n-1}(\boldsymbol{\theta})+w_{n}f^{v_{n}}(\boldsymbol{\theta})$.
To do so, the problem is reduced to showing that

$\displaystyle\mathbb{E}\left[\mathbb{E}[f^{v_{n}}(\boldsymbol{\theta}_{n-a_{n}})-\bar{f}_{n}(\boldsymbol{\theta}_{n-a_{n}})|\mathcal{F}_{n-a_{n}}]^{+}\right]=O(w_{n-a_{n}})+O(\lambda^{a_{n}})$
(36)

(see [34, Proposition 8.1]) which is similar to (31). By additionally assuming
Lipschitz continuity of the individual loss functions $f^{v}$, it is then
shown that

$\displaystyle\mathbb{E}\left[\mathbb{E}[f^{v_{n}}(\boldsymbol{\theta}_{n})-\bar{f}_{n}(\boldsymbol{\theta}_{n})|\mathcal{F}_{n-a_{n}}]\right]=O(w_{n-a_{n}})+O(\lambda^{a_{n}})+O(\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-a_{n}}\rVert)$
(37)

which may be compared with (32). Finally, (35) is proved by showing

$\displaystyle\sum_{n=1}^{\infty}w_{n}\mathbb{E}[|\bar{g}_{n}(\boldsymbol{\theta}_{n})-\bar{f}_{n}(\boldsymbol{\theta}_{n})|]\leq\bar{g}_{1}(\boldsymbol{\theta}_{1})+\sum_{n=1}^{\infty}w_{n}\mathbb{E}[f^{v_{n}}(\boldsymbol{\theta}_{n})-\bar{f}_{n}(\boldsymbol{\theta}_{n})]$
(38)

and using the bound (37) to conclude that the latter sum is finite.

MISO and its extensions proposed in this paper are similar to SMM in their use
of the majorization-minimization principle but contain a few key differences.
Put shortly, the original implementation of MISO is

$\displaystyle\textbf{(MISO)}:\begin{cases}\text{Sample }v_{n}\text{ from the
conditional distribution }\pi(\cdot|\mathcal{F}_{n-1})\\\
g_{n}^{v_{n}}\leftarrow\text{Convex surrogate of $f^{v_{n}}$ at
$\boldsymbol{\theta}_{n-1}$; $g_{n}^{v}=g_{n-1}^{v}$ for $v\neq v_{n}$}\\\
\boldsymbol{\theta}_{n}\in\operatorname*{arg\,min}_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left(\bar{g}_{n}(\boldsymbol{\theta}):=\sum_{v\in\mathcal{V}}g_{n}^{v}\pi(v)\right)\end{cases}$
(39)

In contrast with SMM, each surrogate defining $\bar{g}_{n}$ is given a
constant weight. Consequently, the additional control provided by the
decreasing weights $(w_{n})$ in SMM is not present, which makes the adaptation
of the techniques used in the analysis of SMM non-trivial.

A key step in the original analysis of MISO given in [36] is to show

$\displaystyle\sum_{n=1}^{\infty}\mathbb{E}[\bar{h}_{n}(\boldsymbol{\theta}_{n})]<\infty.$
(40)

This is similar to (35) and shows that the averaged surrogate is an
asymptotically accurate approximation of the true objective at all
$\boldsymbol{\theta}_{n}$s. To prove this, one needs to relate the averaged
error $\bar{h}_{n}(\boldsymbol{\theta}_{n})$ to another quantity proven to be
summable through other means. This is, in analogy, the role that mixing rate
analysis plays for SGD and SMM. Using techniques from [36] one can prove

$\displaystyle\sum_{n=1}^{\infty}\mathbb{E}[h_{n}^{v_{n+1}}(\boldsymbol{\theta}_{n})]<\infty.$
(41)

If the $(v_{n})$ are drawn i.i.d from $\pi$, or more generally if the
probability of transitioning between any two vertices is uniformly bounded
below by a positive quantity, then $h_{n}^{v_{n+1}}(\boldsymbol{\theta}_{n})$
is a conditionally unbiased estimate of $\bar{h}_{n}(\boldsymbol{\theta}_{n})$
up to a constant. Then (40) follows by conditioning on the most recent
information $\mathcal{F}_{n}$.

To adapt the analysis to a more general setting, one may attempt to use Markov
chain mixing to show

$\displaystyle|\bar{h}_{n}(\boldsymbol{\theta}_{n-a_{n}})-\mathbb{E}[h_{n}^{v_{n+1}}(\boldsymbol{\theta}_{n-a_{n}})|\mathcal{F}_{n-a_{n}}]|=O(\lambda^{a_{n}})$
(42)

similar to (31) and (36). However, an additional complication arises because
the function $h_{n}^{v_{n+1}}$ conditional on $\mathcal{F}_{n-a_{n}}$ depends
on both the history of data samples $v_{k}$ and the estimated parameters
$\boldsymbol{\theta}_{k}$ for $n-a_{n}\leq k\leq n$. In contrast, $f^{v_{n}}$
for SGD (in (31)) depends only on the last data sample $v_{n}$ and
$\bar{f}_{n}$ for SMM (in (36)) depends only on the histrory of data samples
$v_{n-a_{n}},\dots,v_{n}$. So, one cannot use the Markov property to isolate
the randomness in $h_{n}^{v_{n+1}}(\boldsymbol{\theta}_{n-a_{n}})$ due to the
Markov chain transition over the interval $[n-a_{n},n+1]$. To alleviate this
problem, one may attempt to control the difference

$\displaystyle|h_{n}^{v_{n+1}}(\boldsymbol{\theta}_{n-a_{n}})-h_{n-a_{n}}^{v_{n+1}}(\boldsymbol{\theta}_{n-a_{n}})|$
(43)

but doing so is not straightforward without access to something resembling the
weights $(w_{n})$ in SMM.

###  4.2. Sketch of our approach

In this section, we provide a short sketch of our analysis in order to convey
the main ideas. Instead of showing that (40) holds we focus directly on
$\lVert\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n})\rVert$ and prove

$\displaystyle\sum_{n=1}^{N}c_{n}\mathbb{E}[\lVert\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n})\rVert]=O\left(\left(\sum_{n=1}^{N}c_{n}^{2}\right)^{1/2}\right),$
(44)

where $c_{n}$ is any non-increasing sequence. For simplicity, assume that the
surrogate functions $g_{n}^{v}$ and the objective functions $f^{v}$ are
differentiable and we are in the unconstrained setting, i.e.
$\boldsymbol{\Theta}=\mathbb{R}^{p}$. Then if each iterate is produced
according to (39), $\boldsymbol{\theta}_{n}$ is a minimizer of $\bar{g}_{n}$
so $\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})=0$. Therefore
$\lVert\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n})\rVert=\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})-\nabla
f(\boldsymbol{\theta}_{n})\rVert=\lVert\nabla
f(\boldsymbol{\theta}_{n})\rVert$ and (44) implies

$\displaystyle\sum_{n=1}^{N}c_{n}\mathbb{E}[\lVert\nabla
f(\boldsymbol{\theta}_{n})\rVert]=O\left(\left(\sum_{n=1}^{N}c_{n}^{2}\right)^{1/2}\right).$
(45)

If we take $c_{n}=1$, we can then conclude that $\min_{1\leq n\leq
N}\mathbb{E}[\lVert\nabla f(\boldsymbol{\theta}_{n})\rVert]=O(N^{-1/2})$.

In order to prove (44) in the presence of recurrent data samples, we use the
regularization schemes introduced earlier to obtain the following iterate
stability

$\displaystyle\sum_{n=1}^{\infty}\mathbb{E}[\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}]<\infty,$
(46)

an idea also used in [28, 33, 34]. The addition of regularization introduces
an added complication because we are no longer directly minimizing
$\bar{g}_{n}$ on the entire feasible set and so do not have
$\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})=0$. However, as we argue in
Section 7, the added regularization is not too strong asymptotically.

Returning to $\eqref{eq: sketch_gap_sum}$, the main idea is to focus on the
individual gradient errors $\nabla h_{n}^{v}(\boldsymbol{\theta}_{n})$. The
advantage here is that the error at each individual surrogate has the
additional property of being zero at the most recent point of relaxation. That
is, if $g_{n}^{v}$ is a surrogate function of $f^{v}$ at
$\boldsymbol{\theta}_{k}$ for $k\leq n-1$, then $\nabla
h_{n}^{v}(\boldsymbol{\theta}_{k})=0$. Because $h_{n}^{v}$ has $L$-Lipschitz
continuous gradients, we have the inequality $\lVert\nabla
h_{n}^{v}(\boldsymbol{\theta}_{n})\rVert\leq
L\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{k}\rVert$. Since
$k^{v}(n)$ is the most recent time $v$ was visited by the sampling process, it
is also the last time the surrogate at $v$ has been updated. Then $g_{n}^{v}$
is a surrogate of $f^{v}$ at $\boldsymbol{\theta}_{k^{v}(n)-1}$. We then have
$\lVert\nabla h_{n}^{v}(\boldsymbol{\theta}_{n})\rVert\leq
L\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{k^{v}(n)-1}\rVert$, so to
show (44) we only need to prove

$\displaystyle\sum_{n=1}^{N}c_{n}\mathbb{E}[\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{k^{v}(n)-1}\rVert]=O\left(\left(\sum_{n=1}^{N}c_{n}^{2}\right)^{1/2}\right)$
(47)

By the triangle inequality and monotonicity of $(c_{n})$,

$\displaystyle
c_{n}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{k^{v}(n)-1}\rVert\leq
c_{n}\sum_{i=k^{v}(n)}^{n}\lVert\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{i-1}\rVert\leq\sum_{i=k^{v}(n)}^{n}c_{i}\lVert\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{i-1}\rVert,$
(48)

so we can relate the error $\lVert\nabla
h_{n}^{v}(\boldsymbol{\theta}_{n})\rVert$ to the sequence
$(\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2})_{n\geq
1}$ which is known to be summable.

Under the recurrence assumption A1 one can expect $\mathbb{E}[n-k^{v}(n)]\leq
M$ for some $M$. One can then intuitively view the expectation of the right
most quantity in (48) similarly to

$\displaystyle\sum_{i=n-M}^{n}c_{i}\mathbb{E}[\lVert\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{i-1}\rVert].$
(49)

So we have

$\displaystyle\sum_{n=1}^{N}\sum_{i=n-M}^{n}c_{i}\mathbb{E}[\lVert\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{i-1}\rVert]\approx
M\sum_{n=1}^{N}c_{n}\mathbb{E}[\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert]=O\left(\left(\sum_{n=1}^{N}c_{n}^{2}\right)^{1/2}\right),$
(50)

where the last equality is due to Cauchy-Schwartz and the iterate stability
(46). Full details of our analysis are given in Sections 6, 7, and 8.

##  5\. Applications and Experiments

###  5.1. Applications

In this subsection, we give some applications of our general framework. These
include applications to matrix factorization as well as a double-averaging
version of RMISO derived by using prox-linear surrogates.

####  5.1.1. Distributed Matrix Factorization

Before beginning this section we need to define some additional notation. For
a finite collection of matrices $\\{A_{i}\\}_{i\in I}$ each having the same
number of rows, we let $[A_{i};i\in I]$ denote the matrix formed by
concatenating the $A_{i}$ along the horizontal axis. For example, if $I=[n]$
and each $A_{i}$ is in $\mathbb{R}^{p\times d}$ then $[A_{i};i\in
I]=[A_{1},\ldots,A_{n}]\in\mathbb{R}^{p\times nd}$. For a set
$\boldsymbol{\Theta}\subseteq\mathbb{R}^{p\times d}$, we define
$\boldsymbol{\Theta}^{I}:=\\{[A_{i};i\in I]:A_{i}\in\boldsymbol{\Theta}\text{
for all }i\\}$. We denote $\boldsymbol{\Theta}^{I}=\boldsymbol{\Theta}^{n}$ if
$I=[n]$.

We consider the matrix factorization loss $f(W,H)=\frac{1}{2}\lVert
X-WH\rVert_{F}^{2}+\alpha\lVert H\rVert_{1}$ where $X\in\mathbb{R}^{p\times
d}$ is a given data matrix to be factored into the product of dictionary
$W\in\boldsymbol{\Theta}_{W}\subseteq\mathbb{R}^{p\times r}$ and code
$H\in\boldsymbol{\Theta}_{H}\subseteq\mathbb{R}^{r\times d}$ with $\alpha\geq
0$ being the $L_{1}$-regularization parameter for $H$. Here
$\boldsymbol{\Theta}_{W}$ and $\boldsymbol{\Theta}_{H}$ are convex constraint
sets.

Suppose we have a connected graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$
where each vertex stores a collection of samples $\\{X^{i,v}\\}_{i=1}^{k_{v}}$
of data matrices. For simplicity, further suppose the vertices all store an
equal amount of samples, i.e. $k_{v}\equiv k$ for all $v\in\mathcal{V}$. Let
$X^{v}=[X^{1,v},\dots,X^{k,v}]\in\mathbb{R}^{p\times kd}$. For each
$v\in\mathcal{V}$ define the loss function

$\displaystyle
f^{v}(W):=\sum_{i=1}^{k}\inf_{H\in\boldsymbol{\Theta}_{H}}\left(\frac{1}{2}\lVert
X^{i,v}-WH\rVert_{F}^{2}+\alpha\lVert
H\rVert_{1}\right)=\inf_{H\in\boldsymbol{\Theta}_{H}^{k}}\frac{1}{2}\lVert
X^{v}-WH\rVert_{F}^{2}+\alpha\lVert H\rVert_{1},$ (51)

which is the minimum reconstruction error for factorizing $X_{v}$ using the
dictionary $W$. In this context, the empirical loss minimization problem
becomes

$\displaystyle
W^{*}\in\operatorname*{arg\,min}_{W\in\boldsymbol{\Theta}_{W}}\frac{1}{|\mathcal{V}|}\sum_{v\in\mathcal{V}}f^{v}(W),$
(52)

where $\boldsymbol{\Theta}_{W}\subseteq\mathbb{R}^{p\times r}$ is a convex
constraint set. Note that this problem is not convex. Indeed, letting
$m=k|\mathcal{V}|$ and $X=[X_{v};v\in V]$ it may be written as

$\displaystyle(W^{*},H^{*})\in\operatorname*{arg\,min}_{W\in\boldsymbol{\Theta}_{W},H\in\boldsymbol{\Theta}_{H}^{m}}\frac{1}{2}\lVert
X-WH\rVert_{F}^{2}+\alpha\lVert H\rVert_{1},$ (53)

which is a constrained non-convex optimization problem with a bi-convex loss
function.

In order to apply RMISO suppose $W^{\prime}\in\boldsymbol{\Theta}_{W}$ is
fixed and denote
$H^{\prime}\in\operatorname*{arg\,min}_{H\in\boldsymbol{\Theta}_{H}^{k}}\lVert
X_{v}-W^{\prime}H\rVert+\alpha\lVert H\rVert_{1}$. Then the function
$g^{v}(W):=\frac{1}{2}\lVert X_{v}-WH^{\prime}\rVert_{F}^{2}+\alpha\lVert
H^{\prime}\rVert_{1}$ is a majorizing surrogate of $f^{v}$ at $W^{\prime}$ and
belongs to $\mathcal{S}_{L^{\prime}}(f^{v},W^{\prime})$ for some
$L^{\prime}>0$ (see Ex. B.5). Then Algorithms 1 and 2 reduce to computing

$\displaystyle
W_{n}\in\operatorname*{arg\,min}_{W\in\boldsymbol{\Theta}_{W}}\bar{g}_{n}(W)+\frac{\rho_{n}}{2}\lVert
W-W_{n-1}\rVert_{F}^{2}\quad\text{ and }\quad
W_{n}\in\operatorname*{arg\,min}_{W\in\boldsymbol{\Theta}_{W}\cap
B_{r_{n}}(W_{n-1})}\bar{g}_{n}(W),$ (54)

respectively. Here

$\displaystyle\bar{g}_{n}(W)=\frac{1}{|\mathcal{V}|}\sum_{v\in\mathcal{V}}g_{n}^{v}(W),\quad
g_{n}^{v}(W)=\frac{1}{2}\lVert X_{v}-WH_{n}^{v}\rVert_{F}^{2}+\alpha\lVert
H_{n}^{v}\rVert_{1}$ (55)

and

$\displaystyle\begin{cases}H_{n}^{v}\in\operatorname*{arg\,min}_{H\in\boldsymbol{\Theta}_{H}^{k}}\frac{1}{2}\lVert
X_{v}-W_{n-1}H\rVert_{F}^{2}+\alpha\lVert H\rVert_{1}&v_{n}=v\\\
H_{n}^{v}=H_{n-1}^{v}&\text{otherwise}.\end{cases}$ (56)

One can show that minimizing the averaged surrogate is equivalent to
minimizing

$\displaystyle\text{tr}(WA_{n}W^{T})-2\text{tr}(WB_{n}),$ (57)

with $A_{n}$ and $B_{n}$, defined recursively as

$\displaystyle A_{n}$
$\displaystyle:=A_{n-1}+\frac{1}{|\mathcal{V}|}\Big{[}H_{n}^{v_{n}}(H_{n}^{v_{n}})^{T}-H_{n-1}^{v_{n}}(H_{n-1}^{v_{n}})^{T}\Big{]}$
(58) $\displaystyle B_{n}$
$\displaystyle:=B_{n-1}+\frac{1}{|\mathcal{V}|}\Big{[}H_{n}^{v_{n}}X_{v}^{T}-H_{n-1}^{v_{n}}X_{v}^{T}\Big{]}$
(59)

The algorithms are re-stated formally in Appendix C. The above algorithm may
be compared to the celebrated online matrix factorization algorithm in [37],
which is a special instance of SMM (34).

####  5.1.2. Prox-Linear Surrogates

Suppose each $f^{v}$ is differentiable and has $L$-Lipschitz continuous
gradients. Then the functions

$\displaystyle
g_{n}^{v}(\boldsymbol{\theta})=f^{v}(\boldsymbol{\theta}_{k^{v}(n)-1})+\langle\nabla
f^{v}(\boldsymbol{\theta}_{k^{v}(n)-1}),\boldsymbol{\theta}-\boldsymbol{\theta}_{k^{v}(n)-1}\rangle+\frac{L}{2}\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{k^{v}(n)-1}\rVert^{2}.$
(60)

are in $\mathcal{S}_{2L}(f^{v},\boldsymbol{\theta}_{k^{v}(n)-1})$ (see Example
B.2) . Further, suppose that $\pi$ is the uniform distribution. Then the
update according to Algorithm 1 is

$\begin{cases}\bar{\boldsymbol{\theta}}_{n-1}&\leftarrow\frac{1}{|\mathcal{V}|}\sum_{v\in\mathcal{V}}\boldsymbol{\theta}_{k^{v}(n)-1}\\\
\bar{\nabla}_{n-1}&\leftarrow\frac{1}{|\mathcal{V}|}\sum_{v\in\mathcal{V}}\nabla
f^{v}(\boldsymbol{\theta}_{k^{v}(n)-1})\\\
\tilde{\boldsymbol{\theta}}_{n-1}&\leftarrow\frac{\rho_{n}}{L+\rho_{n}}\boldsymbol{\theta}_{n-1}+\frac{L}{L+\rho_{n}}\bar{\boldsymbol{\theta}}_{n-1}\\\
\boldsymbol{\theta}_{n}&\leftarrow\text{Proj}_{\boldsymbol{\Theta}}\left(\tilde{\boldsymbol{\theta}}_{n-1}-\frac{1}{L+\rho_{n}}\bar{\nabla}_{n-1}\right).\end{cases}$
(61)

By setting $\rho_{n}\equiv 0$ for all $n$ we get the original MISO iteration
with prox-linear surrogates [36]:

$\displaystyle\begin{cases}\bar{\boldsymbol{\theta}}_{n-1}\leftarrow\frac{1}{|\mathcal{V}|}\sum_{v\in\mathcal{V}}\boldsymbol{\theta}_{k^{v}(n)-1}\\\
\bar{\nabla}_{n-1}\leftarrow\frac{1}{|\mathcal{V}|}\sum_{v\in\mathcal{V}}\nabla
f^{v}(\boldsymbol{\theta}_{k^{v}(n)-1})\\\
\boldsymbol{\theta}_{n}\leftarrow\text{Proj}_{\boldsymbol{\Theta}}\left(\bar{\boldsymbol{\theta}}_{n-1}-\frac{1}{L}\bar{\nabla}_{n-1}\right).\end{cases}$
(62)

Comparing MISO (62) with RMISO (61), we see that the additional proximal
regularization has the effect of further averaging the iterates, putting
additional weight of $\frac{\rho_{n}}{L+\rho_{n}}$ on the most recent
parameter $\boldsymbol{\theta}_{n-1}$. Our algorithm can also be compared to
SAG [53]

$\begin{cases}\bar{\nabla}_{n-1}&\leftarrow\frac{1}{|\mathcal{V}|}\sum_{v\in\mathcal{V}}\nabla
f^{v}(\boldsymbol{\theta}_{k^{v}(n)-1})\\\
\boldsymbol{\theta}_{n}&\leftarrow\text{Proj}_{\boldsymbol{\Theta}}\left(\boldsymbol{\theta}_{n-1}-\alpha_{n}\bar{\nabla}_{n-1}\right)\end{cases}$
(63)

where the most recent iterate $\boldsymbol{\theta}_{n-1}$ is used rather than
the averaged iterate $\tilde{\boldsymbol{\theta}}_{n-1}$. Though not supported
by our theory, it is interesting to point out that RMISO can be seen as an
interpolation between MISO and SAG. In fact if we set $L=0$ in (61), we get
(63) with stepsize $\frac{1}{\rho_{n}}$.

###  5.2. Experiments

In this section, we provide experimental results for Algorithms 1 and 2. We
consider a distributed non-negative matrix factorization problem for the MNIST
data set [12] and a binary classification problem for the a9a data set [5].

####  5.2.1. Distributed Nonnegative Matrix Factorization

In this section, we compare the performance of the distributed matrix
factorization version of RMISO from Section 5.1.1 against other well-known
methods for solving the problem (53). We consider a randomly drawn collection
of $n=5000$ samples where each sample $X_{i}\in\mathbb{R}^{28\times 28}_{+}$
represents one image in the MNIST data set [12]. In all experiments, we set
the $L_{1}$ regularization parameter $\alpha=\frac{1}{28}$ and the number of
components $p=15$. The feasible set is

$\displaystyle\boldsymbol{\Theta}_{W}=\\{W\in\mathbb{R}_{+}^{p\times r}:\lVert
w_{j}\rVert\leq 1\text{ for }1\leq j\leq p\\}$ (64)

meaning we require $W$ to have non-negative entries and each row to have norm
and most one.

The set of vertices $\mathcal{V}$ is arranged in a cycle graph with
$|\mathcal{V}|=55$. Each vertex is restricted to only contain data samples
with the same label. In order to demonstrate RMISO’s dependence on $t_{\odot}$
and $t_{\text{hit}}$ we consider two different sampling algorithms: the first
is the standard random walk on the graph and the second is cyclic sampling
where the vertices are sampled in order without replacement. For the random
walk both $t_{\odot}$ and $t_{\text{hit}}$ are $O(|\mathcal{V}|^{2})$ whereas
for cyclic sampling both are $O(|\mathcal{V}|)$. Our theory suggests we should
expect better performance for cyclic sampling versus the random walk.

Random Walk |  Cyclic
---|---
|
Figure 2. Plot of reconstruction error against both times for NMF using two
sampling algorithms. Results display the performance of algorithms RMISO,
MISO, ONMF, and AdaGrad in factorizing a collection of data matrices from the
MNIST dataset [12].

We compare RMISO with MISO [36] which is Algorithm 1 with $\rho=0$, as well as
the online nonnegative matrix factorization (ONMF) algorithm in [37] and
AdaGrad [13] with constant step size $\eta=0.5$. Our experiments include all
three versions of RMISO: constant proximal regularization (RMISO-CPR), dynamic
proximal regularization (RMISO-DPR), and diminishing radius (RMISO-DR). For
constant proximal regularization, we set $\rho=2500$ for the random walk and
$\rho=50$ for cyclic sampling. For the diminishing radius version of RMISO we
set $r_{n}=\frac{1}{\sqrt{n}\log(n+1)}$.

We plot the reconstruction error versus compute time in Figure 2. We see that
RMISO outperforms ONMF and shows competitive performance against AdaGrad in
both sampling schemes. As expected, there is a dramatic improvement in
performance under cyclic sampling versus random walk sampling.

####  5.2.2. Binary classification with nonconvex loss

We consider binary classification with the non-convex loss function

$\displaystyle\ell_{\boldsymbol{\theta}}(x,y)=\left(1-\frac{1}{1+\exp(-y(x^{T}\boldsymbol{\theta}))}\right)^{2}$
(65)

where $y\in\\{-1,1\\}$ is the class label, $x\in\mathbb{R}^{d}$ is a feature
vector, and $\boldsymbol{\theta}\in\mathbb{R}^{d}$ is the parameter to be
optimized. We use the a9a dataset [5]. Here we consider two separate graph
topologies: the complete graph and the "lonely" graph as in Figure 1. Both
graphs have $|\mathcal{V}|=50$ and each vertex only stores data with the same
label. We sampled the vertices according to a simple random walk on the graph.

Complete Graph |  Lonely Graph
---|---
|
Figure 3. Plot of objective loss against the test dataset for a9a for two
graph topologies and various optimization algorithms- RMISO, AdaGrad, MCSAG,
SGD, Adam, and SGD-HB.

We compared five different optimization algorithms: (1) The prox-linear
version of Algorithm 1 (61); (2) AdaGrad [13]; (3) Markov Chain SAG (MCSAG)
[15], which is the SAG update (63) with an adaptive step size policy where
$\alpha_{n}=\frac{1}{L(t_{\text{hit}}+\max_{v\in\mathcal{V}}(n-k^{v}(n))}$;
(4) SGD with decaying step size $\alpha_{n}=\frac{\alpha}{n^{\gamma}}$; and
(5)SGD-HB (SGD with momentum). AdaGrad and SGD-HB were run with stepsize
$\alpha=0.05$ while SGD with a decaying stepsize was run with a constant step
size parameter $\alpha=0.1$ and decay parameter $\gamma=0.5$. For MCSAG and
RMISO we took $L=1$. The random walk on the complete graph has
$t_{\text{hit}}=O(|\mathcal{V}|)$ while $t_{\text{hit}}=O(|\mathcal{V}|^{2})$
for the lonely graph but $t_{\odot}=O(|\mathcal{V}|)$ for both. Accordingly,
for MCSAG we set the hitting time parameter in the step size
$t_{\text{hit}}=50$ for the complete graph and $t_{\text{hit}}=2500$ for the
lonely graph. For RMISO we set $\rho=50$ for both the constant proximal
regularization version and the dynamic proximal regularization version.

The results are plotted in Figure 3 showing the loss against the test dataset
for both experiments. In the complete graph setting, AdaGrad and SGD-HB
converge quickly but suffer from unstable trajectories while MCSAG and RMISO
are much more stable. Overall though all algorithms display comparable
performance with the constant proximal regularization version of RMISO only
maintaining a slight advantage over the dynamic proximal regularization
version and MCSAG. In the lonely graph scenario however, RMISO-DPR and MCSAG
display poorer performance as $t_{\text{hit}}$ increases from
$O(|\mathcal{V}|)$ to $O(|\mathcal{V}|^{2})$. The performance of RMISO-CPR is
unchanged since it only depends on $t_{\odot}$.

##  6\. Preliminary Lemmas

In this section, we state and prove some preliminary lemmas which will be used
to prove in the proofs of both Theorem 3.1 and 3.4.

We first introduce some additional notation. Throughout this section as well
as the remainder of the paper we let

$\displaystyle t_{\text{cov}}:=\sup_{n\geq
1}\left\lVert\mathbb{E}\left[\max_{v\in\mathcal{V}}\tau_{n,v}\Big{|}\mathcal{F}_{n}\right]\right\rVert_{\infty}.$
(66)

We recall that the $L_{\infty}$ norm here is taken for the conditional
expectation viewed as a random variable. This quantity is a generalization of
the worst-case expected cover time from Markov chain theory [31, Ch. 11] and
will be important in the analysis of Algorithm 1 with dynamic proximal
regularization.

Our first result, Proposition 6.1, states that under A1 the return times have
finite moments of all orders and give an upper bound on $t_{\text{cov}}$ in
terms of $t_{\text{hit}}$. The first item will be used in Section 8 to prove
asymptotic results while the second is used in the proof of iteration
complexity for case C2.

######  Proposition 6.1 (Recurrence implies finite exponential moments of
return time).

Assume A1 and let $t_{\text{hit}}$ and $t_{\text{cov}}$ be as in (13) and (66)
respectively.

  1. (i)

There exists $s_{0}>0$ so that for all $0<s<s_{0}$ there is a constant
$C_{s}>0$ with

$\displaystyle\max_{v\in\mathcal{V}}\,\sup_{n\geq
1}\,\,\mathbb{E}\left[e^{s\tau_{n,v}}\Big{|}\mathcal{F}_{n}\right]\leq
C_{s}<\infty.$ (67)

Consequently, for each $p\geq 1$, $\max_{v\in\mathcal{V}}\sup_{n\geq
1}\mathbb{E}[\tau_{n,v}^{p}|\mathcal{F}_{n}]\leq C<\infty$ for some $C>0$.

  2. (ii)

We have the following bound on $t_{\text{cov}}$:

$\displaystyle t_{\text{cov}}\leq(2t_{\text{hit}}+1)\log_{2}(4|\mathcal{V}|).$
(68)

###### Proof.

Let $m$ be the smallest integer satisfying $m\geq 2t_{\text{hit}}$. For any
$n\geq 1$ and $v\in\mathcal{V}$ notice that if $\tau_{n,v}\geq km$ then we
must have $\tau_{n+jm,v}\geq m$ for each $0\leq j\leq k-1$. Then

$\displaystyle\mathbb{P}(\tau_{n,v}\geq
km|\mathcal{F}_{n})\leq\mathbb{P}\left(\bigcap_{j=0}^{k-1}\\{\tau_{n+jm,v}\geq
m\\}\Big{|}\mathcal{F}_{n}\right)=\mathbb{E}\left[\prod_{j=0}^{k-1}\mathbbm{1}(\tau_{n+jm,v}\geq
m)\Big{|}\mathcal{F}_{n}\right].$ (69)

We have

$\displaystyle\mathbb{E}\left[\prod_{j=0}^{k-1}\mathbbm{1}(\tau_{n+jm,v}\geq
m)\Big{|}\mathcal{F}_{n}\right]=\mathbb{E}\left[\mathbb{E}\big{[}\mathbbm{1}(\tau_{n+(k-1)m,v}\geq
m)|\mathcal{F}_{n+(k-1)m}\big{]}\prod_{j=0}^{k-2}\mathbbm{1}(\tau_{n+jm,v}\geq
m)\Big{|}\mathcal{F}_{n}\right]$ (70)

where we have used that $\\{\tau_{n+jm,v}\geq m\\}$ is measurable with respect
to $\mathcal{F}_{n+(k-1)m}$ for each $j\leq k-2$. By Markov’s inequality and
A1

$\displaystyle\mathbb{E}\big{[}\mathbbm{1}(\tau_{n+(k-1)m,v}\geq
m)|\mathcal{F}_{n+(k-1)m}\big{]}=\mathbb{P}\left(\tau_{n+(k-1)m,v}\geq
m|\mathcal{F}_{n+(k-1)m}\right)\leq\frac{t_{\text{hit}}}{m}.$

So

$\displaystyle\mathbb{E}\left[\prod_{j=0}^{k-1}\mathbbm{1}(\tau_{n+jm,v}\geq
m)\Big{|}\mathcal{F}_{n}\right]\leq\frac{t_{\text{hit}}}{m}\mathbb{E}\left[\prod_{j=0}^{k-2}\mathbbm{1}(\tau_{n+jm,v}\geq
m)\Big{|}\mathcal{F}_{n}\right].$

Proceeding by induction it follows that

$\displaystyle\mathbb{P}(\tau_{n,v}\geq
km|\mathcal{F}_{n})\leq\left(\frac{t_{\text{hit}}}{m}\right)^{k}\leq 2^{-k}$

with the second inequality using our choice of $m$. Now,

$\displaystyle\mathbb{E}\left[e^{s\tau_{n,v}}|\mathcal{F}_{n}\right]=\sum_{\ell=1}^{\infty}e^{s\ell}\mathbb{P}(\tau_{n,v}=\ell|\mathcal{F}_{n})\leq\sum_{k=0}^{\infty}e^{s(k+1)m}\mathbb{P}(\tau_{n,v}\geq
km|\mathcal{F}_{n})\leq\sum_{k=0}^{\infty}e^{s(k+1)m}2^{-k}.$ (71)

The latter sum is finite if $s<\frac{\log 2}{2m}$ and does not depend on $n$
or $v$ which shows (i).

With $n$ still fixed let $\tau_{\text{cov}}=\max_{v\in\mathcal{V}}\tau_{n,v}$.
We have

$\displaystyle\mathbb{E}[\tau_{\text{cov}}|\mathcal{F}_{n}]=\sum_{\ell=1}^{\infty}\mathbb{P}(\tau_{\text{cov}}\geq\ell|\mathcal{F}_{n})\leq\sum_{k=0}^{\infty}m\mathbb{P}(\tau_{\text{cov}}\geq
km|\mathcal{F}_{n})\leq(2t_{\text{hit}}+1)\sum_{k=0}^{\infty}\mathbb{P}(\tau_{\text{cov}}\geq
km|\mathcal{F}_{n})$ (72)

since $\mathbb{P}(\tau_{\text{cov}}\geq\ell|\mathcal{F}_{n})$ is a decreasing
function of $\ell$ and $m\leq 2t_{\text{hit}}+1$. By a union bound

$\displaystyle\mathbb{P}(\tau_{\text{cov}}\geq km|\mathcal{F}_{n})\leq
1\wedge\sum_{v\in\mathcal{V}}\mathbb{P}(\tau_{n,v}\geq km|\mathcal{F}_{n})\leq
1\wedge|\mathcal{V}|2^{-k}.$ (73)

Summing a geometric series we get

$\displaystyle\sum_{k=0}^{\infty}\mathbb{P}(\tau_{\text{cov}}\geq
km|\mathcal{F}_{n})\leq\log_{2}|\mathcal{V}|+|\mathcal{V}|\sum_{k>\log_{2}\mathcal{|V|}}2^{-k}\leq\log_{2}|\mathcal{V}|+2.$
(74)

Combining (72) and (74) shows (ii). ∎

The next proposition states some general properties of first-order surrogate
functions.

######  Proposition 6.2 (Properties of Surrogates).

Fix $\bar{\boldsymbol{\theta}}\in\boldsymbol{\Theta}$ and
$f:\boldsymbol{\Theta}\to\mathbb{R}$. Let
$g\in\mathcal{S}_{L}(f,\bar{\boldsymbol{\theta}})$ and let
$\boldsymbol{\theta}^{\prime}$ be a minimizer of $g$ over
$\boldsymbol{\Theta}$. Then for all
$\boldsymbol{\theta}\in\boldsymbol{\Theta}$

$\displaystyle|h(\boldsymbol{\theta})|\leq\frac{L}{2}\lVert\boldsymbol{\theta}-\bar{\boldsymbol{\theta}}\rVert^{2}$
(75)

###### Proof.

This follows from using the classical upper bound for $L$-smooth functions

$\displaystyle h(\boldsymbol{\theta})\leq
h(\bar{\boldsymbol{\theta}})+\langle\nabla
h(\bar{\boldsymbol{\theta}}),\boldsymbol{\theta}-\bar{\boldsymbol{\theta}}\rangle+\frac{L}{2}\lVert\boldsymbol{\theta}-\bar{\boldsymbol{\theta}}\rVert^{2}$
(76)

(see Lemma A.1) and noting that $h(\bar{\boldsymbol{\theta}})$ and $\nabla
h(\bar{\boldsymbol{\theta}})$ are both equal to zero according to Definition
2.1.

∎

Next, we show that the surrogate objective value
$\bar{g}_{n}(\boldsymbol{\theta}_{n})$ evaluated at $\boldsymbol{\theta}_{n}$
is non-increasing.

######  Lemma 6.3 (Surrogate Monotonicity).

Let $(\boldsymbol{\theta}_{n})_{n\geq 0}$ be an output of either Algorithm 1
or 2 . Then
$\bar{g}_{n}(\boldsymbol{\theta}_{n-1})\leq\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1})$
for all $n\geq 1$. Moreover, the sequence
$(\bar{g}_{n}(\boldsymbol{\theta}_{n}))_{n\geq 0}$ is non-increasing. As a
consequence, by Assumption A2 and Definition 2.1,
$(\bar{g}_{n}(\boldsymbol{\theta}_{n}))_{n\geq 0}$ is bounded below with
probability one and therefore
$\lim_{n\to\infty}\bar{g}_{n}(\boldsymbol{\theta}_{n})$ exists almost surely.

###### Proof.

Since $g_{n}^{v_{n}}\in\mathcal{S}_{L}(f^{v_{n}},\boldsymbol{\theta}_{n-1})$,
Definition 2.1 implies that
$g_{n}^{v_{n}}(\boldsymbol{\theta}_{n-1})=f^{v_{n}}(\boldsymbol{\theta}_{n-1})$.
Then

$\displaystyle\bar{g}_{n}(\boldsymbol{\theta}_{n-1})$
$\displaystyle=\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1})+\Big{[}g_{n}^{v_{n}}(\boldsymbol{\theta}_{n-1})-g_{n-1}^{v_{n}}(\boldsymbol{\theta}_{n-1})\Big{]}\pi(v_{n})$
(77)
$\displaystyle=\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1})+\Big{[}f^{v_{n}}(\boldsymbol{\theta}_{n-1})-g_{n-1}^{v_{n}}(\boldsymbol{\theta}_{n-1})\Big{]}\pi(v_{n})$
(78) $\displaystyle\leq\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1})$ (79)

where the last inequality used $g_{n-1}^{v_{n}}$ is a majorizing surrogate of
$f^{v_{n}}$.

Suppose now that $(\boldsymbol{\theta}_{n})_{n\geq 0}$ is an output of
Algorithm 1. Then by definition of $\boldsymbol{\theta}_{n}$,

$\displaystyle\bar{g}_{n}(\boldsymbol{\theta}_{n})$
$\displaystyle\leq\bar{g}_{n}(\boldsymbol{\theta}_{n})+\frac{\rho_{n}}{2}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}$
$\displaystyle\leq\bar{g}_{n}(\boldsymbol{\theta}_{n-1})+\frac{\rho_{n}}{2}\lVert\boldsymbol{\theta}_{n-1}-\boldsymbol{\theta}_{n-1}\rVert^{2}$
$\displaystyle=\bar{g}_{n}(\boldsymbol{\theta}_{n-1})$
$\displaystyle\leq\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1}).$

If instead $(\boldsymbol{\theta}_{n})_{n\geq 1}$ is an output of Algorithm 2,
then we can directly conclude
$\bar{g}_{n}(\boldsymbol{\theta}_{n})\leq\bar{g}_{n}(\boldsymbol{\theta}_{n-1})$
by definition of $\boldsymbol{\theta}_{n}$. The remainder of the proof is
identical to the above. ∎

The next lemma establishes the summability of the sequence
$h_{n}^{v_{n+1}}(\boldsymbol{\theta}_{n})$. This was used in [36] to prove
asymptotic convergence of MISO under i.i.d. sampling. We use it primarily in
the analysis of Algorithm 2.

######  Lemma 6.4.

Let $(\boldsymbol{\theta}_{n})_{n\geq 0}$ be an output of either Algorithm 1
or 2. Then almost surely

$\displaystyle\sum_{n=1}^{\infty}h_{n}^{v_{n+1}}(\boldsymbol{\theta}_{n})\leq\frac{1}{\pi_{\text{min}}}\Delta_{0}.$
(80)

###### Proof.

By Definition 2.1, for each $n$ the quantity
$h_{n}^{v_{n+1}}(\boldsymbol{\theta}_{n})$ is non-negative. Therefore, it
suffices to show that the sequence of partial sums,
$\sum_{n=1}^{N}h_{n}^{v_{n+1}}(\boldsymbol{\theta}_{n})$ is uniformly bounded.

Recall that

$\displaystyle\bar{g}_{n+1}(\boldsymbol{\theta}_{n+1})\leq\bar{g}_{n+1}(\boldsymbol{\theta}_{n})$
$\displaystyle=\bar{g}_{n}(\boldsymbol{\theta}_{n})+(g_{n+1}^{v_{n+1}}(\boldsymbol{\theta}_{n})-g_{n}^{v_{n+1}}(\boldsymbol{\theta}_{n}))\pi(v_{n+1})$
(81)
$\displaystyle=\bar{g}_{n}(\boldsymbol{\theta}_{n})+(f^{v_{n+1}}(\boldsymbol{\theta}_{n})-g_{n}^{v_{n+1}}(\boldsymbol{\theta}_{n}))\pi(v_{n+1}).$
(82)

We then have

$\displaystyle\sum_{n=1}^{N}h_{n}^{v_{n+1}}(\boldsymbol{\theta}_{n})$
$\displaystyle=\sum_{n=1}^{N}g_{n}^{v_{n+1}}(\boldsymbol{\theta}_{n})-f^{v_{n+1}}(\boldsymbol{\theta}_{n})$
$\displaystyle\leq\sum_{n=1}^{N}\frac{1}{\pi(v_{n+1})}\left(\bar{g}_{n}(\boldsymbol{\theta}_{n})-\bar{g}_{n+1}(\boldsymbol{\theta}_{n+1})\right)$
$\displaystyle\leq\frac{1}{\pi_{\text{min}}}\sum_{n=1}^{N}\bar{g}_{n}(\boldsymbol{\theta}_{n})-\bar{g}_{n+1}(\boldsymbol{\theta}_{n+1})$
$\displaystyle\leq\frac{1}{\pi_{\text{min}}}\Delta_{0}$

which is what we needed to show. ∎

######  Proposition 6.5.

Suppose $(\boldsymbol{\theta}_{n})_{n\geq 0}$ is an output of Algorithm 2. The
for all $n\geq 1$,
$\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert\leq r_{n}$.

###### Proof.

This follows directly from the definition of $\boldsymbol{\theta}_{n}$ in
Algorithm 2. ∎

Lemma 6.6 establishes the iterate stability. These results are crucially used
to control the surrogate error gradient
$\lVert\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n})\rVert$ in Lemma 7.1 as well
as in the asymptotic analysis of RMISO in Section 8.

######  Lemma 6.6 (Finite variation of iterate differences).

The following hold almost surely:

  1. (i)

For case C3,

$\displaystyle\sum_{n=1}^{\infty}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}\leq\sum_{n=1}^{\infty}r_{n}^{2}<\infty.$
(83)

  2. (ii)

In either of the cases C1 or C2,

$\displaystyle\sum_{n=1}^{\infty}\frac{\rho_{n}}{2}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}\leq\Delta_{0}<\infty.$
(84)

###### Proof.

The proof of (i) can be deduced from Proposition 6.5 and Assumption A3.

Now assume either of the cases C1 or C2. By definition of
$\boldsymbol{\theta}_{n}$,

$\displaystyle\bar{g}_{n}(\boldsymbol{\theta}_{n})+\frac{\rho_{n}}{2}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}\leq\bar{g}_{n}(\boldsymbol{\theta}_{n-1})+\frac{\rho_{n}}{2}\lVert\boldsymbol{\theta}_{n-1}-\boldsymbol{\theta}_{n-1}\rVert^{2}=\bar{g}_{n}(\boldsymbol{\theta}_{n-1})\leq\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1})$
(85)

where the last inequality is due to Lemma 6.3. So

$\displaystyle\frac{\rho_{n}}{2}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}\leq\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1})-\bar{g}_{n}(\boldsymbol{\theta}_{n}).$
(86)

Hence,

$\displaystyle\sum_{n=1}^{N}\frac{\rho_{n}}{2}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}\leq\sum_{n=1}^{N}\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1})-\bar{g}_{n}(\boldsymbol{\theta}_{n})$
$\displaystyle=\bar{g}_{0}(\boldsymbol{\theta}_{0})-\bar{g}_{N}(\boldsymbol{\theta}_{N})\leq\Delta_{0}.$
(87)

Letting $N\to\infty$ shows that

$\displaystyle\sum_{n=1}^{\infty}\frac{\rho_{n}}{2}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}\leq\Delta_{0}$
(88)

as desired. This shows (ii). ∎

The remaining results in this section concern Algorithm 2 and are used both in
the convergence rate analysis in Section 7 as well as the asymptotic analysis
of Section 8. Recall that in this case, we are assuming that $\nabla f^{v}$ is
$L$-Lipschitz continuous for each $v\in\mathcal{V}$. Proposition 6.7 states
that this assumption implies $\nabla\bar{g}_{n}$ is differentiable and $2L$
Lipschitz for each $n$.

######  Proposition 6.7.

Let $\\{\boldsymbol{\theta}_{v}\\}_{v\in\mathcal{V}}$ be a collection of
$|\mathcal{V}|$ points in $\boldsymbol{\Theta}$. Suppose the optional
condition A2’ and that
$g_{n}^{v}\in\mathcal{S}_{L}(f^{v},\boldsymbol{\theta}_{v})$ for each $v$.
Then

  1. (i)

The gradient of the objective function $\nabla f=\sum_{v\in\mathcal{V}}\nabla
f^{v}\pi(v)$ is $L$-Lipschitz over $\boldsymbol{\Theta}$.

  2. (ii)

For each $v$, $\nabla g_{n}^{v}$ is $2L$-Lipchitz over $\boldsymbol{\Theta}$.
In addition, $\nabla\bar{g}_{n}$ is $2L$-Lipchitz.

###### Proof.

Since $\pi$ is a probability distribution, (i) follows easily from the
triangle inequality.

For (ii) note that $\nabla(g_{n}^{v}-f^{v})=\nabla h_{n}^{v}$ is $L$-Lipshitz
by Definition 2.1. Then since $\nabla g_{n}^{v}=\nabla h_{n}^{v}+\nabla f^{v}$
it follows from the triangle inequality that $\nabla g_{n}^{v}$ is
$2L$-Lipschitz. Then recalling that
$\nabla\bar{g}_{n}=\sum_{v\in\mathcal{V}}\nabla g_{n}^{v}\pi(v)$ another
application of the triangle inequality shows that $\nabla\bar{g}_{n}$ is
$2L$-Lipschitz continuous. ∎

######  Proposition 6.8.

Assume case C3 and let $(\boldsymbol{\theta}_{n})_{n\geq 0}$ be an output of
Algorithm 2. Then

$\displaystyle\sum_{n=1}^{N}|\langle\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rangle|\leq\Delta_{0}+L\sum_{n=1}^{N}r_{n}^{2}$
(89)

almost surely.

###### Proof.

Since $\nabla\bar{g}_{n}$ is $2L$-Lipschitz continuous by Proposition 6.7, by
Lemma A.1

$\displaystyle|\bar{g}_{n}(\boldsymbol{\theta}_{n})-\bar{g}_{n}(\boldsymbol{\theta}_{n-1})-\langle\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rangle|\leq
L\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}.$ (90)

Under Algorithm 2 we have
$\bar{g}_{n}(\boldsymbol{\theta}_{n})\leq\bar{g}_{n}(\boldsymbol{\theta}_{n-1})$
by the definition of $\boldsymbol{\theta}_{n}$ and
$\bar{g}_{n}(\boldsymbol{\theta}_{n-1})\leq\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1})$
by Lemma 6.3. These observations together with (90) imply

$\displaystyle\left|\langle\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rangle|\right|$
$\displaystyle\leq|\bar{g}_{n}(\boldsymbol{\theta}_{n-1})-\bar{g}_{n}(\boldsymbol{\theta}_{n})|+L\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}$
(91)
$\displaystyle=\bar{g}_{n}(\boldsymbol{\theta}_{n-1})-\bar{g}_{n}(\boldsymbol{\theta}_{n})+L\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}$
(92)
$\displaystyle\leq\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1})-\bar{g}_{n}(\boldsymbol{\theta}_{n})+L\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}.$
(93)

We have

$\displaystyle\sum_{n=1}^{N}\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1})-\bar{g}_{n}(\boldsymbol{\theta}_{n})\leq\Delta_{0}$
(94)

almost surely. Therefore

$\displaystyle\sum_{n=1}^{N}\left|\langle\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rangle\right|$
$\displaystyle\leq\sum_{n=1}^{N}\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1})-\bar{g}_{n}(\boldsymbol{\theta}_{n})+L\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}$
(95)
$\displaystyle\leq\Delta_{0}+L\sum_{n=1}^{N}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}$
(96) $\displaystyle\leq\Delta_{0}+L\sum_{n=1}^{N}r_{n}^{2},$ (97)

where the last line uses
$\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert\leq r_{n}$. ∎

The next lemma is a key to establishing the iteration complexity of Algorithm
2. A similar lemma was used to analyze block majorization-minimization with
diminishing radius in [28].

######  Lemma 6.9 (Approximate first order optimality).

Let $(\boldsymbol{\theta}_{n})_{n\geq 0}$ be an output of Algorithm 2 and let
$b_{n}=\min\\{1,r_{n}\\}$. Then

$\displaystyle
b_{n}\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert\leq
1}\langle-\nabla\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rangle\leq\langle-\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rangle+r_{n}\lVert\nabla
h_{n-1}^{v_{n}}(\boldsymbol{\theta}_{n-1})\rVert+2Lr_{n}^{2}.$ (98)

###### Proof.

Fix $\boldsymbol{\theta}\in\boldsymbol{\Theta}$ with
$\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert\leq b_{n}$. By
definition of $\boldsymbol{\theta}_{n}$ we have
$\bar{g}_{n}(\boldsymbol{\theta}_{n})\leq\bar{g}_{n}(\boldsymbol{\theta})$.
Subtracting $\bar{g}_{n}(\boldsymbol{\theta}_{n-1})$ from both sides and using
Proposition 6.7 and Lemma A.1,

$\displaystyle\langle\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rangle-L\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}$
$\displaystyle\leq\bar{g}_{n}(\boldsymbol{\theta}_{n})-\bar{g}_{n}(\boldsymbol{\theta}_{n-1})$
(99)
$\displaystyle\leq\bar{g}_{n}(\boldsymbol{\theta})-\bar{g}_{n}(\boldsymbol{\theta}_{n-1})$
(100)
$\displaystyle\leq\langle\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rangle+L\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert^{2}.$
(101)

Notice that

$\displaystyle\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n-1})$
$\displaystyle=\nabla\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1})+\left[\nabla
g_{n}^{v_{n}}(\boldsymbol{\theta}_{n-1})-\nabla
g_{n-1}^{v_{n}}(\boldsymbol{\theta}_{n-1})\right]\pi(v_{n})$ (102)
$\displaystyle=\nabla\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1})+\left[\nabla
f^{v_{n}}(\boldsymbol{\theta}_{n-1})-\nabla
g_{n-1}^{v_{n}}(\boldsymbol{\theta}_{n-1}\right]\pi(v_{n})$ (103)
$\displaystyle=\nabla\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1})-\nabla
h_{n-1}^{v_{n}}(\boldsymbol{\theta}_{n-1})\pi(v_{n}).$ (104)

where the second line used
$g_{n}^{v_{n}}\in\mathcal{S}_{L}(f^{v_{n}},\boldsymbol{\theta}_{n-1})$ and
item (ii) of Definition 2.1, and the third line used the definition of
$h_{n-1}^{v_{n}}$. Therefore, adding and subtracting
$\langle\nabla\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rangle$
from the right hand side of (101) we get

$\displaystyle\langle\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rangle$
$\displaystyle\leq\langle\nabla\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rangle-\pi(v_{n})\langle\nabla
h_{n-1}^{v_{n}}(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rangle+L\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert^{2}$
(105)
$\displaystyle\qquad\qquad\qquad\qquad+L\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}$
(106)
$\displaystyle\leq\langle\nabla\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rangle+\lVert\nabla
h_{n-1}^{v_{n}}(\boldsymbol{\theta}_{n-1})\rVert\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert+L\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert^{2}$
(107)
$\displaystyle\qquad\qquad\qquad\qquad+L\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}$
(108)
$\displaystyle\leq\langle\nabla\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rangle+r_{n}\lVert\nabla
h_{n-1}^{v_{n}}(\boldsymbol{\theta}_{n-1})\rVert+2Lr_{n}^{2}$ (109)

where the last line used
$\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert\leq r_{n}$ and
$\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert\leq b_{n}\leq
r_{n}$. Since the above holds for all
$\boldsymbol{\theta}\in\boldsymbol{\Theta}$ with
$\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert\leq b_{n}$ we
obtain

$\displaystyle\langle\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rangle\leq\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert\leq
b_{n}}\langle\nabla\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rangle+r_{n}\lVert\nabla
h_{n-1}^{v_{n}}(\boldsymbol{\theta}_{n-1})\rVert+2Lr_{n}^{2}.$ (110)

Finally notice that since $b_{n}\leq 1$, the convexity of
$\boldsymbol{\Theta}$ implies that
$\boldsymbol{\theta}_{n-1}+b_{n}(\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1})\in\boldsymbol{\Theta}$
for any $\boldsymbol{\theta}\in\boldsymbol{\Theta}$. Thus, if there exists
$\boldsymbol{\theta}\in\boldsymbol{\Theta}$ with
$\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert\leq 1$ there is
$\boldsymbol{\theta}^{\prime}\in\boldsymbol{\Theta}$ with
$\lVert\boldsymbol{\theta}^{\prime}-\boldsymbol{\theta}_{n-1}\rVert\leq b_{n}$
such that the direction of
$\boldsymbol{\theta}^{\prime}-\boldsymbol{\theta}_{n-1}$ agrees with that of
$\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}$. Therefore

$\displaystyle
b_{n}\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert\leq
1}\langle-\nabla\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rangle$
$\displaystyle=\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert\leq
1}\langle-\nabla\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1}),b_{n}(\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1})\rangle$
(111)
$\displaystyle\leq\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert\leq
b_{n}}\langle-\nabla\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rangle.$
(112)

combining this with (110) we complete the proof. ∎

##  7\. Convergence Rate Analysis

In this subsection, we prove the convergence rate guarantees of Theorem 3.1.

###  7.1. The key lemma

First, we state and prove Lemma 7.1 which lies at the heart of our analysis.
It allows us to relate the surrogate error gradient
$\lVert\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n})\rVert$ to the sequence of
parameter differences
$(\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2})$ which is
known to be summable by Lemma 6.6. It is important to note that we only use
the recurrence of data sampling A1 and the structure of the algorithm in the
proof.

######  Lemma 7.1 (Key lemma).

Let $(c_{n})_{n\geq 1}$ be a non-increasing sequence of positive numbers. For
any of the cases C1-C3 and any $v\in V$,

$\displaystyle\mathbb{E}\left[\sum_{n=1}^{N}c_{n}\lVert\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n})\rVert\right]\leq
Lt_{\odot}\left(\sum_{n=1}^{N}c_{n}^{2}\right)^{1/2}\mathbb{E}\left[\left(\sum_{n=1}^{N}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}\right)^{1/2}\right].$
(113)

###### Proof.

Fix some $v\in V$. We recall that $k^{v}(n)$ is the last time before $n$ that
the sampling process visited $v$ and therefore the last time the surrogate
$g_{n}^{v}$ was updated. We then have
$g_{n}^{v}\in\mathcal{S}_{L}(f^{v},\boldsymbol{\theta}_{k^{v}(n)-1})$ so by
the definition of first-order surrogates (Definition 2.1) $\nabla
h_{n}^{v}(\boldsymbol{\theta}_{k^{v}(n)-1})=0$. Combining this with the
Lipschitz continuity of $\nabla h_{n}^{v}$ we get

$\displaystyle c_{n}\lVert\nabla
h_{n}^{v}(\boldsymbol{\theta}_{n})\rVert=c_{n}\lVert\nabla
h_{n}^{v}(\boldsymbol{\theta}_{n})-\nabla
h_{n}^{v}(\boldsymbol{\theta}_{k^{v}(n)-1})\rVert\leq
Lc_{n}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{k^{v}(n)-1}\rVert\leq
L\sum_{i=k^{v}(n)}^{n}c_{n}\lVert\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{i-1}\rVert.$
(114)

Therefore by the triangle inequality,

$\displaystyle c_{n}\lVert\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n})\rVert\leq
c_{n}\sum_{v\in\mathcal{V}}\lVert\nabla
h_{n}^{v}(\boldsymbol{\theta}_{n})\rVert\pi(v)\leq
L\sum_{v\in\mathcal{V}}\left(\sum_{i=k^{v}(n)}^{n}c_{n}\lVert\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{i-1}\rVert\right)\pi(v).$
(115)

For an integer $n$, let $p^{v}(n)=\inf\\{j>n:v_{j}=v\\}$ be the first time
strictly after $n$ that the sampling algorithm visits $v$. Denote $a\land
b:=\min(a,b)$. We have

$\displaystyle\sum_{n=1}^{N}\sum_{v\in\mathcal{V}}\left(\sum_{i=k^{v}(n)}^{n}c_{n}\lVert\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{i-1}\rVert\right)\pi(v)$
$\displaystyle=\sum_{v\in\mathcal{V}}\left(\sum_{n=1}^{N}\sum_{i=k^{v}(n)}^{n}c_{n}\lVert\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{i-1}\rVert\right)\pi(v)$
(116)
$\displaystyle=\sum_{v\in\mathcal{V}}\left(\sum_{i=1}^{N}\sum_{n=i}^{N\wedge
p^{v}(i)-1}c_{n}\lVert\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{i-1}\rVert\right)\pi(v)$
(117)
$\displaystyle\leq\sum_{v\in\mathcal{V}}\left(\sum_{i=1}^{N}c_{i}\lVert\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{i-1}\rVert(p^{v}(i)-i)\right)\pi(v)$
(118)
$\displaystyle=\sum_{v\in\mathcal{V}}\left(\sum_{i=1}^{N}c_{i}\lVert\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{i-1}\rVert\tau_{i,v}\right)\pi(v),$
(119)

where the third line used that $(c_{n})$ is non-increasing. So we get

$\displaystyle\mathbb{E}\left[\sum_{n=1}^{N}c_{n}\lVert\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n})\rVert\right]$
$\displaystyle\leq
L\mathbb{E}\left[\sum_{v\in\mathcal{V}}\left(\sum_{n=1}^{N}c_{n}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert\tau_{n,v}\right)\pi(v)\right]$
(120)
$\displaystyle=L\mathbb{E}\left[\sum_{n=1}^{N}c_{n}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert\left(\sum_{v\in
V}\tau_{n,v}\pi(v)\right)\right]$ (121)
$\displaystyle=L\mathbb{E}\left[\sum_{n=1}^{N}c_{n}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert\left(\sum_{v\in\mathcal{V}}\mathbb{E}[\tau_{n,v}|\mathcal{F}_{n}]\pi(v)\right)\right]$
(122) $\displaystyle\leq
Lt_{\odot}\mathbb{E}\left[\sum_{n=1}^{N}c_{n}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert\right]$
(123) $\displaystyle\leq
Lt_{\odot}\left(\sum_{n=1}^{N}c_{n}^{2}\right)^{1/2}\mathbb{E}\left[\left(\sum_{n=1}^{N}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}\right)^{1/2}\right]$
(124)

with the second to last line using A1 and the last using Cauchy-Schwartz. ∎

###  7.2. The constant proximal regularization case C1

In this section, we prove Theorem 3.1 for case C1.

###### Proof of Theorem 3.1 for case C1.

We first use the linearity of the limit and the differentiability of the
average surrogate error $\bar{h}_{n}$ to get

$\displaystyle\left|\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})-\nabla
f(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\right|=\left|\langle\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n}),\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rangle\right|\leq\lVert\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n})\rVert\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq\lVert\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n})\rVert$
(125)

for all $\boldsymbol{\theta}\in\boldsymbol{\Theta}$ with
$\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq 1$. It follows
from the triangle inequality, taking supremum, and then expectations that

$\displaystyle\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}-\nabla
f(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\right]\leq\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}-\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\right]+\mathbb{E}\left[\lVert\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n})\rVert\right].$
(126)

Our goal will be to control the sum of the right hand side.

We first address the first term on the right-hand side of (126). For any
$\boldsymbol{\theta}\in\boldsymbol{\Theta}$

$\displaystyle\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})+\rho\langle\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\geq
0$ (127)

since $\boldsymbol{\theta}_{n}$ is a minimizer of
$\bar{g}_{n}(\boldsymbol{\theta})+\frac{\rho}{2}\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert^{2}$
over $\boldsymbol{\Theta}$. Then

$\displaystyle-\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\leq\rho\langle\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1},\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rangle\leq\rho\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq\rho\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert$
(128)

for any $\boldsymbol{\theta}\in\boldsymbol{\Theta}$ with
$\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq 1$. Therefore,

$\displaystyle\mathbb{E}\left[\sum_{n=1}^{N}\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}-\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\right]$
$\displaystyle\leq\rho\mathbb{E}\left[\sum_{n=1}^{N}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert\right]$
(129)
$\displaystyle\leq\rho\sqrt{N}\mathbb{E}\left[\left(\sum_{n=1}^{N}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}\right)^{1/2}\right]$
(130)

where we used the Cauchy-Schwartz inequality in the last line. By Lemma 6.6

$\displaystyle\left(\sum_{n=1}^{N}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}\right)^{1/2}\leq\sqrt{\frac{2\Delta_{0}}{\rho}}$
(131)

almost surely. Thus,

$\displaystyle\mathbb{E}\left[\sum_{n=1}^{N}\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}-\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\right]\leq\sqrt{2N\rho\Delta_{0}}.$
(132)

We now turn to the second term on the right-hand side of (126). By Lemma 7.1
with $c_{n}=1$ and Lemma 6.6

$\displaystyle\mathbb{E}\left[\sum_{n=1}^{N}\lVert\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n})\rVert\right]\leq\sqrt{N}Lt_{\odot}\mathbb{E}\left[\left(\sum_{n=1}^{N}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}\right)^{1/2}\right]\leq\sqrt{\frac{2N\Delta_{0}}{\rho}}Lt_{\odot}.$
(133)

Now, summing both sides of (126) and using (132) and (133),

$\displaystyle\sum_{n=1}^{N}\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}-\nabla
f(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\right]$
$\displaystyle=\mathbb{E}\left[\sum_{n=1}^{N}\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}-\nabla
f(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\right]$
(134)
$\displaystyle\leq\sqrt{2N\Delta_{0}}\left(\sqrt{\rho}+\frac{Lt_{\odot}}{\sqrt{\rho}}\right)$
(135)

so

$\displaystyle\min_{1\leq n\leq
N}\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}-\nabla
f(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\right]\leq\frac{\sqrt{2\Delta_{0}}\left(\sqrt{\rho}+\frac{Lt_{\odot}}{\sqrt{\rho}}\right)}{\sqrt{N}}.$
(136)

The proof is complete.

∎

###  7.3. The dynamic proximal regularization case C2

In this section, we prove Theorem 3.1 for case C2. Recall the definition of
$t_{\text{cov}}$ from (66). Before proving the theorem we introduce a Lemma
adapted from [15, Lemma 4]. This is used to bound the expected sum of the
first $N$ dynamic regularization parameters $\rho_{n}$ in terms of
$t_{\text{cov}}$.

######  Lemma 7.2.

Let $a_{n}=\max_{v\in V}(n-k^{v}(n))$. Then

$\displaystyle\sum_{n=1}^{N}\mathbb{E}[a_{n}]\leq Nt_{\text{cov}}.$ (137)

###### Proof.

Since $a_{n}\leq n-1$ we have

$\displaystyle\sum_{n=1}^{N}\mathbb{E}[a_{n}]=\sum_{n=1}^{N}\sum_{i=1}^{n-1}\mathbb{P}(a_{n}\geq
i).$ (138)

Let $b_{n}=\max_{v\in\mathcal{V}}\tau_{n,v}$. We note that if $a_{n}\geq i$
then there is $v\in V$ with $v_{j}\neq v$ for all $n-i<j\leq n$ and so
$\tau_{n-i,v}\geq i$. So we have the inclusion $\\{a_{n}\geq
i\\}\subseteq\\{b_{n-i}\geq i\\}$. Therefore

$\displaystyle\sum_{n=1}^{N}\sum_{i=1}^{n-1}\mathbb{P}(a_{n}\geq i)$
$\displaystyle=\sum_{i=1}^{N}\sum_{n=i+1}^{N}\mathbb{P}(a_{n}\geq i)$ (139)
$\displaystyle\leq\sum_{i=1}^{N}\sum_{n=i+1}^{N}\mathbb{P}(b_{n-i}\geq i)$
(140) $\displaystyle=\sum_{s=1}^{N-1}\sum_{t=1}^{N-s}\mathbb{P}(b_{s}\geq t)$
(141)
$\displaystyle\leq\sum_{s=1}^{N-1}\sum_{t=1}^{\infty}\mathbb{P}(b_{s}\geq t)$
(142) $\displaystyle\leq N\sup_{n\geq 1}\mathbb{E}[b_{n}].$ (143)

We have

$\displaystyle\mathbb{E}[b_{n}]=\mathbb{E}\left[\mathbb{E}\left[\max_{v\in\mathcal{V}}\tau_{n,v}\Big{|}\mathcal{F}_{n}\right]\right]\leq
t_{\text{cov}}$ (144)

so we are done. ∎

###### Proof of Theorem 3.1 for case C2.

The proof in this case follows the same strategy as case C1, but is slightly
more complicated due to the randomness of the dynamic proximal regularization
parameter.

Define
$G_{n}:=\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1})-\bar{g}_{n}(\boldsymbol{\theta}_{n})$.
By optimality of $\boldsymbol{\theta}_{n}$,

$\displaystyle\bar{g}_{n}(\boldsymbol{\theta}_{n})+\frac{\rho_{n}}{2}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}\leq\bar{g}_{n-1}(\boldsymbol{\theta}_{n-1})$
(145)

so
$\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert\leq\sqrt{2\rho_{n}^{-1}G_{n}}$.
Using similar reasoning as in the proof for case C1

$\displaystyle-\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\leq\rho_{n}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert\leq\sqrt{2\rho_{n}G_{n}}.$
(146)

We have used Cauchy-Schwartz twice,

$\displaystyle\sum_{n=1}^{N}\mathbb{E}[\sqrt{\rho_{n}G_{n}}]$
$\displaystyle\leq\sum_{n=1}^{N}(\mathbb{E}[\rho_{n}])^{1/2}(\mathbb{E}[G_{n}])^{1/2}$
(147)
$\displaystyle\leq\left(\sum_{n=1}^{N}\mathbb{E}[\rho_{n}]\right)^{1/2}\left(\sum_{n=1}^{N}\mathbb{E}[G_{n}]\right)^{1/2}$
(148) $\displaystyle\leq\sqrt{N(\rho+t_{\text{cov}})\Delta_{0}}.$ (149)

The last inequality here uses $\rho_{n}=\rho+\max_{v\in V}(n-k^{v}(n))$ and
Lemma 7.2 as well as $\sum_{n=1}^{N}G_{n}\leq\Delta_{0}$ a.s. It follows that

$\displaystyle\sum_{n=1}^{N}\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert\leq
1}-\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\right]\leq\sqrt{2N(\rho+t_{\text{cov}})\Delta_{0}}.$
(150)

To handle the gradient error
$\lVert\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n})\rVert$ we first use Lemma
6.6 and $\rho_{n}\geq\rho$ to conclude

$\displaystyle\sum_{n=1}^{N}\frac{\rho}{2}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}\leq\sum_{n=1}^{N}\frac{\rho_{n}}{2}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}\leq\Delta_{0}$
(151)

almost surely. It then follows from Lemma 7.1 that

$\displaystyle\sum_{n=1}^{N}\mathbb{E}\left[\lVert\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n})\rVert\right]\leq\sqrt{N}Lt_{\odot}\mathbb{E}\left[\left(\sum_{n=1}^{N}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert\right)^{1/2}\right]\leq\sqrt{\frac{2N\Delta_{0}}{\rho}}Lt_{\odot}.$
(152)

Finally, combining (150) and (152) we get

$\displaystyle\sum_{n=1}^{N}\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}-\nabla
f(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\right]\leq\sqrt{2N\Delta_{0}}\left(\sqrt{\rho+t_{\text{cov}}}+\frac{Lt_{\odot}}{\sqrt{\rho}}\right)$
(153)

and so we deduce

$\displaystyle\min_{1\leq n\leq
N}\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}-\nabla
f(\boldsymbol{\theta}_{n},\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\right]\leq\frac{\sqrt{2\Delta_{0}}\left(\sqrt{\rho+t_{\text{cov}}}+\frac{Lt_{\odot}}{\sqrt{\rho}}\right)}{\sqrt{N}}.$
(154)

We complete the proof by substituting the bound for $t_{\text{cov}}$ in
Proposition 6.1. ∎

###  7.4. The diminishing radius case C3

Here we prove Theorem 3.1 for case C3.

###### Proof of Theorem 3.1 for case C3.

Similar to the proof in case C1 we have

$\displaystyle\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}-\langle\nabla
f(\boldsymbol{\theta}_{n}),\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rangle\right]\leq\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}-\langle\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n}),\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rangle\right]+\mathbb{E}\left[\lVert\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n})\rVert\right]$
(155)

Let $b_{n}=r_{n}\wedge 1$. Then by Lemma 6.9

$\displaystyle\sum_{n=1}^{N}b_{n+1}\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}-\langle\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n}),\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rangle\right]$
(156)
$\displaystyle\leq\sum_{n=1}^{N}\mathbb{E}\left[\langle-\nabla\bar{g}_{n+1}(\boldsymbol{\theta}_{n}),\boldsymbol{\theta}_{n+1}-\boldsymbol{\theta}_{n}\rangle\right]+\sum_{n=1}^{N}r_{n+1}\mathbb{E}\left[\lVert\nabla
h_{n}^{v_{n+1}}(\boldsymbol{\theta}_{n})\rVert\right]+\sum_{n=1}^{N}2Lr_{n+1}^{2}.$
(157)

Because $h_{n}^{v_{n+1}}$ is non-negative and has $L$-Lipschitz continuous
gradients $\lVert\nabla
h_{n}^{v_{n+1}}(\boldsymbol{\theta}_{n})\rVert\leq\sqrt{2Lh_{n}^{v_{n+1}}(\boldsymbol{\theta}_{n})}$
(see Lemma A.2). Then

$\displaystyle\sum_{n=1}^{N}r_{n+1}\mathbb{E}\left[\lVert\nabla
h_{n}^{v_{n+1}}(\boldsymbol{\theta}_{n})\rVert\right]$
$\displaystyle\leq\sum_{n=1}^{N}r_{n+1}\mathbb{E}\left[\sqrt{2Lh_{n}^{v_{n+1}}(\boldsymbol{\theta}_{n})}\right]$
(158)
$\displaystyle\leq\left(\sum_{n=1}^{N}r_{n+1}^{2}\right)^{1/2}\left(\sum_{n=1}^{N}\mathbb{E}\left[2Lh_{n}^{v_{n+1}}(\boldsymbol{\theta}_{n})\right]\right)^{1/2}$
(159)
$\displaystyle\leq\sqrt{\frac{2L\Delta_{0}}{\pi_{\text{min}}}\sum_{n=1}^{N}r_{n+1}^{2}}.$
(160)

Here the second line used Cauchy-Schwartz and then Jensen’s inequality to move
the square inside the expectation and the last line used Lemma 6.4. From
Proposition 6.8,

$\displaystyle\sum_{n=1}^{N}\mathbb{E}\left[\langle-\nabla\bar{g}_{n+1}(\boldsymbol{\theta}_{n}),\boldsymbol{\theta}_{n+1}-\boldsymbol{\theta}_{n}\rangle\right]\leq\Delta_{0}+L\sum_{n=1}^{N}r_{n+1}^{2}$
(161)

so from (157)

$\displaystyle\sum_{n=1}^{N}b_{n+1}\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}-\langle\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n}),\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rangle\right]\leq\Delta_{0}+\sqrt{\frac{2L\Delta_{0}}{\pi_{\text{min}}}\sum_{n=1}^{N}r_{n+1}^{2}}+3L\sum_{n=1}^{N}r_{n+1}^{2}.$
(162)

From Lemma 7.1 with $c_{n}=r_{n+1}$

$\displaystyle\sum_{n=1}^{N}r_{n+1}\mathbb{E}\left[\lVert\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n})\rVert\right]\leq\sum_{v\in
V}\mathbb{E}\left[\sum_{n=1}^{N}r_{n+1}\lVert\nabla
h_{n}^{v}(\boldsymbol{\theta}_{n})\rVert\right]\pi(v)\leq
Lt_{\odot}\left(\sum_{n=1}^{N}r_{n+1}^{2}\right)^{1/2}\mathbb{E}\left[\left(\sum_{n=1}^{N}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}\right)^{1/2}\right].$
(163)

Since $\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert\leq
r_{n}$ and $(r_{n})$ is non-increasing, this bound reduces to

$\displaystyle\sum_{n=1}^{N}r_{n+1}\mathbb{E}\left[\lVert\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n})\rVert\right]\leq
Lt_{\odot}\sum_{n=1}^{N}r_{n}^{2}.$ (164)

Combining (162) and (164) with (155) we have

$\displaystyle\sum_{n=1}^{N}b_{n}\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}-\langle\nabla
f(\boldsymbol{\theta}_{n}),\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rangle\right]\leq\Delta_{0}+\sqrt{\frac{2L\Delta_{0}}{\pi_{\text{min}}}\sum_{n=1}^{N}r_{n}^{2}}+\Big{(}3+t_{\odot}\Big{)}L\sum_{n=1}^{N}r_{n}^{2}$
(165)

so

$\displaystyle\min_{1\leq n\leq
N}\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert\leq
1}-\langle\nabla
f(\boldsymbol{\theta}_{n}),\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rangle\right]\leq\frac{\Delta_{0}+\sqrt{\frac{2L\Delta_{0}}{\pi_{\text{min}}}\sum_{n=1}^{N}r_{n}^{2}}+\Big{(}3+t_{\odot}\Big{)}L\sum_{n=1}^{N}r_{n}^{2}}{\sum_{n=1}^{N}b_{n}}.$
(166)

∎

###  7.5. Proofs of Corollaries

In this section, we prove Corollaries 3.2 and 3.3.

###### Proof of corollary 3.2.

Fix $N$ and let $k$ be the integer recognizing the minimum in (18). If
$\boldsymbol{\Theta}=\mathbb{R}^{p}$ then we may choose
$\boldsymbol{\theta}^{*}$ so that
$\boldsymbol{\theta}^{*}-\boldsymbol{\theta}_{k}=-\frac{\nabla
f(\boldsymbol{\theta}_{k})}{\lVert\nabla f(\boldsymbol{\theta}_{k})\rVert}$.
Thus,

$\displaystyle\min_{1\leq n\leq N}\mathbb{E}[\lVert\nabla
f(\boldsymbol{\theta}_{n})\rVert]\leq\mathbb{E}[\lVert\nabla
f(\boldsymbol{\theta}_{k})\rVert]$
$\displaystyle=\mathbb{E}\left[\langle-\nabla
f(\boldsymbol{\theta}_{k}),\boldsymbol{\theta}^{*}-\boldsymbol{\theta}_{k}\rangle\right]$
(167)
$\displaystyle\leq\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{k}\rVert\leq
1}\langle-\nabla
f(\boldsymbol{\theta}_{k}),\boldsymbol{\theta}-\boldsymbol{\theta}_{k}\rangle\right]=O\left(N^{-1/2}\right).$
(168)

If instead $\boldsymbol{\Theta}\neq\mathbb{R}^{p}$ but the second condition
$\text{dist}(\boldsymbol{\theta}_{k},\partial\boldsymbol{\Theta})\geq c$
holds, we can take $\boldsymbol{\theta}-\boldsymbol{\theta}_{k}=-c\frac{\nabla
f(\boldsymbol{\theta}_{k})}{\lVert\nabla f(\boldsymbol{\theta}_{k})\rVert}$.
In doing so we obtain

$\displaystyle\min_{1\leq n\leq N}\mathbb{E}[\lVert\nabla
f(\boldsymbol{\theta}_{n})\rVert]\leq\mathbb{E}[\lVert\nabla
f(\boldsymbol{\theta}_{k})\rVert]$
$\displaystyle=\frac{1}{c}\mathbb{E}\left[\langle-\nabla
f(\boldsymbol{\theta}_{k}),\boldsymbol{\theta}^{*}-\boldsymbol{\theta}_{k}\rangle\right]$
(169)
$\displaystyle\leq\frac{1}{c}\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{k}\rVert\leq
1}\langle-\nabla
f(\boldsymbol{\theta}_{k}),\boldsymbol{\theta}-\boldsymbol{\theta}_{k}\rangle\right]=O\left(N^{-1/2}\right).$
(170)

This shows (25). The proof of (26) is similar. ∎

###### Proof of corollary 3.3.

The convergence rates of Theorem 3.1 are of order $O(N^{-1/2})$ for Algorithm
1. Then we can prove (i) by choosing $N$ large enough so that
$N^{-1/2}\leq\varepsilon$.

If we take $r_{n}=\frac{1}{\sqrt{n}\log n}$ in Algorithm 2 then its
corresponding convergence rate in 3.1 is of order $O(n^{-1/2}\log n)$. Then
(ii) follows by using the fact that
$n\geq\varepsilon^{-2}(3\log\varepsilon^{-1})^{2}$ implies $n^{-1/2}\log
n\leq\varepsilon$ for sufficiently large $\varepsilon$. Indeed since
$\frac{\log x}{\sqrt{x}}$ is decreasing for sufficiently large $x$ we have

$\displaystyle n^{-1/2}\log
n\leq\frac{\varepsilon}{3\log\varepsilon^{-1}}\left(2\log\varepsilon^{-1}+2\log(3\log\varepsilon^{-1})\right)\leq\varepsilon$
(171)

for $\varepsilon$ sufficiently small. ∎

##  8\. Asymptotic Analysis

We use this section to prove Theorem 3.4. Recall that by Proposition 6.1,
$\sup_{n\geq 1}\mathbb{E}[\tau_{n,v}^{2}|\mathcal{F}_{n}]\leq C<\infty$ for
each $v\in\mathcal{V}$. Accordingly, we let $\sigma^{2}=\max_{v\in
V}\sup_{n\geq
1}\left\lVert\mathbb{E}[\tau_{n,v}^{2}|\mathcal{F}_{n}]\right\rVert_{\infty}$.

The first Lemma of this section states that the expected value of the dynamic
regularization parameter $\rho_{n}$ in Algorithm 1 is uniformly bounded. While
$\rho_{n}$ only appears in Algorithm 1, the random variable $\max_{v\in
V}(n-k^{v}(n))$ is also present in the analysis of Algorithm 2. Therefore,
this Lemma is used in the analysis of both algorithms in this section.

######  Lemma 8.1.

Assume A1. Then there is a constant $C>0$ such that

$\displaystyle\sup_{n\geq 1}\,\,\mathbb{E}[\rho_{n}]\leq C.$ (172)

###### Proof.

Fix $v\in\mathcal{V}$. For a positive integer $j$ we have

$\displaystyle\\{n-k^{v}(n)\geq j\\}=\\{k^{v}(n)\leq
n-j\\}\subseteq\\{\tau_{n-j,v}\geq j\\}.$ (173)

Therefore, we get

$\displaystyle\mathbb{E}[(n-k^{v}(n))]$
$\displaystyle=\sum_{j=1}^{\infty}\mathbb{P}(n-k^{v}(n)\geq j)$ (174)
$\displaystyle\leq\sum_{j=1}^{\infty}\mathbb{P}(\tau_{n-j,v}\geq j)$ (175)
$\displaystyle\leq\sum_{j=1}^{\infty}\frac{\mathbb{E}[\tau_{n-j,v}^{2}]}{j^{2}}$
(176) $\displaystyle\leq\sigma^{2}\sum_{j=1}^{\infty}j^{-2}$ (177)

since $\mathbb{E}[\tau_{n-j,v}^{2}]\leq\sigma^{2}$. Finally,

$\displaystyle\mathbb{E}[\rho_{n}]=\rho+\mathbb{E}\left[\max_{v\in\mathcal{V}}(n-k^{v}(n))\right]\leq\rho+\sum_{v\in
V}\mathbb{E}[n-k^{v}(n)]\leq\rho+|\mathcal{V}|\sigma^{2}\sum_{j=1}^{\infty}j^{-2}.$
(178)

∎

###  8.1. The dynamic proximal regularization case C2

Here we prove Theorem 3.4 (i). Our first lemma, Lemma 8.2, is similar to 7.1
and key to showing that
$\lVert\nabla\bar{h}_{n}(\boldsymbol{\theta}_{n})\rVert\to 0$. The difference
is that we must deal with
$\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{k^{v}(n)-1}\rVert^{2}$
instead of
$\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{k^{v}(n)-1}\rVert$. In
order to relate this to the sequence of one-step iterate differences
$(\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2})$ we need
to use the Cauchy-Schwartz inequality which introduces a dependence on
$\sigma^{2}$ as well as $t_{\text{hit}}$.

######  Lemma 8.2.

Let $(\boldsymbol{\theta}_{n})_{n\geq 0}$ be an output of Algorithm 1. Assume
case C2. Then

$\displaystyle\sum_{n=1}^{\infty}\mathbb{E}[\bar{h}_{n}(\boldsymbol{\theta}_{n})]<\infty.$
(179)

###### Proof.

Since $\bar{h}_{n}(\boldsymbol{\theta}_{n})=\sum_{v\in
V}h_{n}^{v}(\boldsymbol{\theta}_{n})\pi(v)$ and $\mathcal{V}$ is finite, it is
sufficient to show
$\sum_{n=1}^{\infty}\mathbb{E}[h_{n}^{v}(\boldsymbol{\theta}_{n})]<\infty$ for
each $v\in V$. Before starting recall that by Lemma 6.6, and
$\rho_{n}\geq\rho>0$

$\displaystyle\sum_{n=1}^{\infty}\frac{\rho_{n}}{2}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}\leq\Delta_{0}\quad\text{
and
}\quad\sum_{n=1}^{\infty}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}\leq\frac{2}{\rho}\Delta_{0}$
(180)

almost surely. This implies

$\displaystyle\mathbb{E}\left[\sum_{n=1}^{\infty}\rho_{n}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}\right]+\mathbb{E}\left[\sum_{n=1}^{\infty}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}\right]<\infty.$
(181)

Fix $v\in V$. For each $n$ we have
$g_{n}^{v}\in\mathcal{S}_{L}(f^{v},\boldsymbol{\theta}_{k^{v}(n)-1})$. Then
using Proposition 6.2, the triangle inequality and Cauchy Schwartz

$\displaystyle|h_{n}^{v}(\boldsymbol{\theta}_{n})|\leq\frac{L}{2}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{k^{v}(n)-1}\rVert^{2}\leq\frac{L}{2}(n-k^{v}(n)+1)\sum_{i=k^{v}(n)}^{n}\lVert\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{i-1}\rVert^{2}.$
(182)

Let
$B_{n}=(n-k^{v}(n)+1)\sum_{i=k^{v}(n)}^{n}\lVert\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{i-1}\rVert^{2}$.
We claim that $\mathbb{E}\left[\sum_{n=1}^{\infty}B_{n}\right]<\infty$. As in
Lemma 7.1 let $p^{v}(n)=\inf\\{j>n:v_{j}=v\\}$ be the next time strictly after
time $n$ the sampling algorithm visits $v$. Exchanging the order of summation
we have

$\displaystyle\mathbb{E}\left[\sum_{n=1}^{\infty}(n-k^{v}(n)+1)\sum_{i=k^{v}(n)}^{n}\lVert\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{i-1}\rVert^{2}\right]$
$\displaystyle=\mathbb{E}\left[\sum_{i=1}^{\infty}\lVert\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{i-1}\rVert^{2}\sum_{n=i}^{p^{v}(i)-1}(n-k^{v}(n)+1)\right]$
(183)
$\displaystyle=\mathbb{E}\left[\sum_{i=1}^{\infty}\lVert\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{i-1}\rVert^{2}\sum_{n=i}^{\infty}(n-k^{v}(n)+1)\mathbbm{1}(p^{v}(i)>n)\right].$
(184)

The equality $\\{p^{v}(i)>n\\}=\\{\tau_{i,v}\geq n-i+1\\}$ holds as both are
equal to the event $\\{v_{j}\neq v:i<j\leq n+1\\}$. Moreover,
$k^{v}(n)=k^{v}(i)$ on $\\{p^{v}(i)>n\\}$ since there is no visit to $v$
between times $i$ and $n$. Therefore,