# Convergence and Complexity of
Stochastic Block Majorization-Minimization

Hanbaek Lyu  Hanbaek Lyu, Department of Mathematics, University of Wisconsin -
Madison, WI 53709, USA<EMAIL_ADDRESS>
###### Abstract.

Stochastic majorization-minimization (SMM) is an online extension of the
classical principle of majorization-minimization, which consists of sampling
i.i.d. data points from a fixed data distribution and minimizing a recursively
defined majorizing surrogate of an objective function. In this paper, we
introduce stochastic block majorization-minimization, where the surrogates can
now be only block multi-convex and a single block is optimized at a time
within a diminishing radius. Relaxing the standard strong convexity
requirements for surrogates in SMM, our framework gives wider applicability
including online CANDECOMP/PARAFAC (CP) dictionary learning and yields greater
computational efficiency especially when the problem dimension is large. We
provide an extensive convergence analysis on the proposed algorithm, which we
derive under possibly dependent data streams, relaxing the standard i.i.d.
assumption on data samples. We show that the proposed algorithm converges
almost surely to the set of stationary points of a nonconvex objective under
constraints at a rate $O((\log n)^{1+\varepsilon}/n^{1/2})$ for the empirical
loss function and $O((\log n)^{1+\varepsilon}/n^{1/4})$ for the expected loss
function, where $n$ denotes the number of data samples processed. Under some
additional assumption, the latter convergence rate can be improved to $O((\log
n)^{1+\varepsilon}/n^{1/2})$. Our results provide first convergence rate
bounds for various online matrix and tensor decomposition algorithms under a
general Markovian data setting.

##  1\. Introduction

Empirical loss minimization is a classical problem setting regarding parameter
estimation with a growing number of observations, where one seeks to minimize
a recursively defined empirical loss function as new data arrives. Some of its
well-known applications include maximum likelihood estimation, or more
generally, $M$-estimation [14, 18, 44], as well as the online dictionary
learning literature [37, 35, 38, 30]. On the other hand, the expected loss
minimization seeks to estimate a parameter by minimizing the loss function
with respect to random data. It provides a general framework for stochastic
optimization literature [45, 36, 1, 43]. Optimization algorithms for empirical
or expected loss minimization are in nature ‘online’, meaning that sampling
new data points and adjusting the current estimation occurs recursively. Such
onilne algorithms have proven to be particularly efficient in large-scale
problems in statistics, optimization, and machine learning [5, 12, 15, 24].

First-order methods for expected loss minimization includes (projected or
proximal) Stochastic Gradient Descent (SGD) have been extensively studied for
many decades, which usually consist of stochastically estimating the full
gradient of the expected loss function and performing a gradient descent step
followed by a projection onto the parameter space. For a general constrained
and nonconvex (weakly convex) setting, the global convergence to stationary
(first-order optimal) points of such methods was recently established in [11]
with rate of convergence $O(\log n/\sqrt{n})$, where $n$ denotes the number of
iterations (processed samples). This result assumes arbitrary initialization,
although faster convergence with special initialization is known in a matrix
factorization setting [51]. Throughout this paper, we are concerned with the
convergence of online algorithms with arbitrary initialization, which is often
referred to as a ‘global convergence’ in the literature.

On the other hand, Stochastic Majorization-Minimization (SMM) is one of the
most popular approaches that directly solve empirical loss minimization [35]
by sampling data points from a target data distribution and minimizing a
recursively defined majorizing surrogate of the empirical loss function. This
method is an online extension of the classical majorization-minimization [27]
principle and encompasses Expectation-Minimization in statistics [42, 10], and
also generalizes the celebrated online matrix factorization algorithm in [37].
While it was observed empirically that SMM shows competitive performance while
requiring less parameter tuning than SGD for online dictionary learning
problems [37], theoretical convergence guarantee of SMM in the literature only
ensures asymptotic convergence to stationary points and lacks any convergence
rate bounds [37, 35, 34, 38].

Most theoretical analyses on online optimization algorithms assume the ability
to sample i.i.d. data points from the target distribution. This is a classical
and convenient assumption for the analysis, but it is violated frequently in
practice, especially when the samples are accessed through Markov chain-based
methods. A common practice is to first sample a single Markov chain trajectory
of dependent data samples (after some burn-out period), which is then thinned
(subsampled) to reduce the dependence. However, such practice only weakens the
dependence between data samples and does make it truly independent in general.
One may re-initialize independent Markov chains for every sample to attain
true independence, but this approach suffers from a huge computational burden.
Optimization algorithms with Markovian data samples were studied in [21, 22]
in the context of distributed optimization in networks. More recently, it was
shown in [46] that arbitrarily initialized SGD almost surely converges to
critical points of unconstrained nonconvex objectives at rate $O((\log
n)^{2}/\sqrt{n})$, even when the data samples have a Markovian dependence.

On the other hand, in [30], the online matrix factorization algorithm in [37]
based on SMM is extended and shown to converge to stationary points for
constrained matrix factorization problems in the Markovian setting. Based on
the result and combining with an MCMC network sampling algorithm in [29], an
online dictionary learning algorithm for learning latent motifs in networks is
proposed in [28]. More recently, an online algorithm for nonnegative tensor
dictionary learning utilizing CANDECOMP/PARAFAC (CP) decomposition is
developed in [32], where convergence to stationary points under Markovian
setting is established. The algorithm is based on SMM but uses new components
such as block coordinate descent with diminishing radius [33]. Although the
works [30, 32] are the first to establish global convergence of SMM-type
methods to stationary points in the Markovian setting, there has not been a
convergence rate bounds, which is still the case in the i.i.d. setting. A rate
of convergence result is important also in practice since it provides bounds
on the number of iterations and data samples sufficient to guarantee a to
obtain an approximate solution up to the desired precision.

###  1.1. Contribution

In this paper, we intend to provide the missing convergence rate analysis of
SMM for general constrained nonconvex objectives, both in the context of
empirical and expected loss minimization. Our analysis is for a generalized
version of SMM, where the surrogates do not have to be convex but only convex
restricted on block coordinates. We also allow using approximate surrogate
functions for robustness under inexact computations for solving convex sub-
problems. For a parameter update, averaged surrogate functions are only block-
minimized either with cyclic or randomized scheduling for more efficient
computation. We call our algorithm stochastic block majorization-minimization
(SBMM), which generalizes the original SMM in [34], online matrix
factorization in [37, 30], stochastic approximate majorization-minimization
and subsampled online matrix factorization in [38], and the online CP-
dictionary learning algorithm in [32].

We provide an extensive convergence analysis on the proposed algorithm, which
we derive under possibly dependent data streams, relaxing the standard i.i.d.
assumption on data samples. We obtain global convergence to stationary points
of rate $O((\log n)^{1+\varepsilon}/n^{1/2})$, matching the optimal
convergence rates for SGD based methods [46, 11]. Interestingly, our analysis
shows that SBMM (and hence SMM) is more adapted to solve empirical loss
minimization than expected loss minimization, in the sense that the
aforementioned rate of convergence holds for the empirical loss functions
almost surely and in expectation for the expected loss function; an almost
sure convergence for the empirical loss function is obtained at a slower rate
of $O((\log n)^{1+\varepsilon}/n^{1/4})$. This is the opposite of SGD based
methods, which converges almost surely for the expected loss at rate $O((\log
n)^{2}/n^{1/2})$ and rate $O((\log n)^{1+\varepsilon}/n^{1/4})$ for the
empirical loss. However, we also show that the optimal almost sure convergence
rate $O((\log n)^{1+\varepsilon}/n^{1/2})$ is obtained when the stationary
points of the surrogate functions are in the interior of the constraint set.
As an application, our result yields the first convergence rate bounds for
online matrix factorization algorithms [37, 30, 38] as well as for online CP-
dictionary learning algorithm [32] under Markovian data setting, which has not
been known even under the i.i.d. data setting.

###  1.2. Related work

In [35], convergence analysis of SMM with i.i.d. data samples is given. There,
it is shown that when the objective function $f$ is convex, SMM converges to
the global minimum with rate of convergence in expectation of order $O(\log
n/\sqrt{n})$ (see [35, Prop 3.1]) and of order $O(1/n)$ when $f$ is strongly
convex (see [35, Prop 3.2]). Also, when $f$ is nonconvex, SMM is shown to
converge almost surely to the set of stationary points of $f$ over a convex
constraint set (see [35, Prop 3.4]) but no result on the rate of convergence
is given. The latter nonconvex result were later extended to SMM with
approximate surrogate functions in [38], which was applied to developing
subsampled online matrix factorization algorithms.

In [30], the classical online matrix factorization algorithm in [37] is shown
to converge when the input data matrices are given by a function of some
underlying Markov chain. A similar almost sure global convergence in Markovian
setting for online tensor factorization is obtained in [32]. In both works,
bounds on the rate of convergence are not given. In [33], block coordinate
descent with diminishing radius is proposed for deterministic nonconvex
problems and shown to converge to the stationary points of the objective
function. Also, a rate of convergence of order $O(\log n/\sqrt{n})$ is
obtained.

Stochastic Gradient Descent (SGD) is another popular method for various
optimization problems. In [46], a convergence of SGD under Markovian data
assumption is obtained. For the convex case, [46, Thm. 1] shows convergence of
SGD to a global minimum with rate of convergence of order $O((\log
n)^{1+\varepsilon}/\sqrt{n})$ for each fixed $\varepsilon>0$ . Moreover, [46,
Thm. 2] shows that SGD for unconstrained nonconvex optimization converge
almost surely to the stationary points with rate $O((\log
n)^{1+\varepsilon}/\sqrt{n})$ for each fixed $\varepsilon>0$. A similar rate
of convergence for nonconvex and constrained projected SGD is also known in
[11]. In [54], an online NMF algorithm based on projected SGD with general
divergence in place of the squared $\ell_{2}$-loss is proposed, and
convergence to stationary points to the expected loss function for i.i.d. data
samples is shown.

###  1.3. Organization

In Section 2, we state the problem settings of empirical and expected loss
minimization and then give a breif introduction to stochastic (block)
majorization-minimization. In Section 3, we define multi-convex approximate
surrogate functions and state our main algorithm in Algorithm 1. Next in
Section 4, we state our main results (Theorems 4.1, 4.2, and 4.3) and their
corollaries (Corollaries 4.4 and 4.5). In Section 5, we list various examples
of multi-convex surrogates and discuss applications of our results in the
contex of online matrix factorization and CP-dictionary learning problems. In
Section 6 we establish some preliminary lemmas. Following in Section 7, we
state five key lemmas (Lemmas 7.1, 7.2, 7.3, 7.4, and 7.5) and derive all main
results stated in Section 4 from them. The first three key lemmas will be
proved in Section 8, whereas the other two key lemmas will be established in
Section 9. Lastly, we give some backgrounds on Markov chains and MCMC sampling
in Appendix A and provide some auxiliary lemmas in Appendix B.

##  2\. Problem statement

###  2.1. Empirical Loss Minimization

Suppose we have a loss function
$\ell:\mathfrak{X}\times\boldsymbol{\Theta}\rightarrow[0,\infty)$ that
measures the fitness of a parameter
$\boldsymbol{\theta}\in\boldsymbol{\Theta}\subseteq\mathbb{R}^{p}$ with
respect to an observed data $\mathbf{x}\in\mathfrak{X}$. Consider a sequence
of newly observed data $(\mathbf{x}_{n})_{n\geq 1}$ in $\mathfrak{X}$. We
would like to estimate a sequence of parameters
$(\boldsymbol{\theta}_{n})_{n\geq 1}$ such that $\boldsymbol{\theta}_{n}$ is
in some sense the best fit to all data $\mathbf{x}_{1},\cdots,\mathbf{x}_{n}$
up to time $n$. In the empirical loss minimization framework, one seeks to
minimize a recursively defined empirical loss function as a new data defined
empirical loss function as a new data arrives:

$\displaystyle\textup{Upon arrival of
$\mathbf{x}_{n}$:}\quad\boldsymbol{\theta}_{n}\in\operatorname*{arg\,min}_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\big{(}\bar{f}_{n}(\boldsymbol{\theta}):=(1-w_{n})\bar{f}_{n-1}(\boldsymbol{\theta})+w_{n}\,\ell(\mathbf{x}_{n},\boldsymbol{\theta})\big{)},$
(1)

where $(w_{n})_{n\geq 1}$ is a sequence of weights in $(0,1]$ and the function
$\bar{f}_{n}$ is the empirical loss function recursively defined by the
weighted average as in (1) with $f_{0}\equiv 0$. More explicitely, we can
write

$\displaystyle\bar{f}_{n}(\boldsymbol{\theta})=\sum_{k=1}^{n}\ell(\mathbf{x}_{k},\boldsymbol{\theta})w^{n}_{k},\qquad
w^{n}_{k}:=w_{k}\prod_{i=k+1}^{n}(1-w_{i}).$ (2)

The weight $w_{n}$ in (1) controls how much we want our new estimate
$\boldsymbol{\theta}_{n}$ deviate from minimizing the previous empirical loss
$\bar{f}_{n-1}$ to adapting to the newly observed tensor data
$\mathbf{x}_{n}$. In one extreme case of $w_{n}\equiv 1$,
$\boldsymbol{\theta}_{n}$ is a minimizer of the time-$n$ loss
$\ell(\mathbf{x}_{n},\cdot)$ and ignores the past $\bar{f}_{n-1}$. If
$w_{n}\equiv\alpha\in(0,1)$ then the history is forgotten exponentially fast,
that is,
$\bar{f}_{n}(\cdot)=\sum_{k=1}^{n}\alpha(1-\alpha)^{n-k}\,\ell(\mathbf{x}_{k},\cdot)$.
On the other hand, the ‘balanced weight’ $w_{n}=1/n$ makes the empirical loss
to be the the arithmetic mean:
$\bar{f}_{n}(\cdot)=\frac{1}{n}\sum_{k=1}^{n}\ell(\mathbf{x}_{s},\cdot)$,
which is the canonical choice in the literature including maximum likelihood
estimation and online NMF problem in [37]. Hence, one may choose the sequence
of weights $(w_{n})_{n\geq 1}$ in (1) to decay fast for learning average
features and decay slow (or keep it constant) for learning trending features).

For instance, consider the linear regression problem where
$\ell(\mathbf{x},\boldsymbol{\theta})=\lVert\mathbf{x}-\mathbf{D}\boldsymbol{\theta}\rVert^{2}$,
where $\mathbf{x}\in\mathbb{R}^{q}$ and $\mathbf{D}\in\mathbb{R}^{q\times p}$
is a fixed matrix of basis features. Then by writing
$\ell(\mathbf{x},\boldsymbol{\theta})=\boldsymbol{\theta}^{T}\mathbf{D}^{T}\mathbf{D}\boldsymbol{\theta}-2\boldsymbol{\theta}^{T}\mathbf{D}^{T}\mathbf{x}+\mathbf{x}^{T}\mathbf{x}$,
it is easy to see that

$\displaystyle\bar{f}_{n}(\boldsymbol{\theta})=\lVert\bar{\mathbf{x}}_{n}-\mathbf{D}\boldsymbol{\theta}\rVert_{F}^{2}+C_{n},$
(3)

where $\bar{\mathbf{x}}_{n}:=\sum_{k=1}^{n}\mathbf{x}_{k}w^{n}_{k}$ and
$C_{n}$ depends only on $\mathbf{x}_{1},\dots,\mathbf{x}_{n}$. Hence
$\boldsymbol{\theta}_{n}$ in (1) in this case is the vector of best linear
regression coefficients that fits the basis $\mathbf{D}$ to the averaged data
point $\bar{\mathbf{x}}_{n}$.

###  2.2. Expected Loss Minimization

Instead of fitting the model underlying the loss function $\ell$ to a sequence
of data points $(\mathbf{x}_{n})_{n\geq 1}$, consider fitting it to a single
but random data point $\mathbf{x}$ following some probability distribution
$\pi$. Then we would be seeking a single estimate $\boldsymbol{\theta}^{*}$
that minimizes the expected loss function $f$ as below:

$\displaystyle\boldsymbol{\theta}^{*}\in\operatorname*{arg\,min}_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left(f(\boldsymbol{\theta}):=\mathbb{E}_{\mathbf{x}\sim\pi}\left[\ell(\mathbf{x},\boldsymbol{\theta})\right]\right).$
(4)

We call the above problem setting expected loss minimization. This is a
popular setting in the optimization literature for stochastic programs. A
popular optimization algorithm for solving (4) is projected stochastic
gradient descent, which proceeds by first drawing a sample
$\mathbf{x}_{n}\sim\pi$, estimating the full gradient $\nabla
f(\boldsymbol{\theta}_{n})$ by the stochastic gradient
$\nabla\ell(\mathbf{x}_{n},\boldsymbol{\theta}_{n-1})$, and then updating
$\boldsymbol{\theta}_{n}\leftarrow\textup{Proj}_{\boldsymbol{\Theta}}(\boldsymbol{\theta}_{n-1}-\alpha_{n}\nabla\ell(\mathbf{x}_{n},\boldsymbol{\theta}_{n-1}))$,
where $\alpha_{n}$ is a chosen stepsize and
$\textup{Proj}_{\boldsymbol{\Theta}}$ is the projection operator onto the
parameter space $\boldsymbol{\Theta}$.

In the context of linear regression we discussed before, the expected loss
function $f$ becomes

$\displaystyle
f(\boldsymbol{\theta})=\lVert\mathbb{E}_{\mathbf{x}\sim\pi}[\mathbf{x}]-\mathbf{D}\boldsymbol{\theta}\rVert_{F}^{2}+C_{n}^{\prime},$
(5)

where $C_{n}^{\prime}$ does not depend on $\boldsymbol{\theta}$ and
$\mathbf{D}$. By comparing the empiricial (3) and the expected (5) loss
functions for linear regression, one can see that the two problem settings are
asymptotically equivalent when
$\bar{\mathbf{x}}_{n}\rightarrow\mathbb{E}_{\mathbf{x}\sim\pi}[\mathbf{x}]$ as
the sample size $n$ for the empricial loss minimization tends to infinity.
This is certainly true when $\mathbf{x}_{k}$’s are i.i.d. from $\pi$ and
$w_{k}=1/k$ by the strong law of large numbers. This holds true in a more
general setting where $\mathbf{x}_{n}$ form a Markov chain with stationary
distribution $\pi$ and the weights $w_{n}$ are non-increasing and decays
sufficiently faster than $1/\sqrt{n}$ (see Lemma B.4).

###  2.3. Stochastic (Block) Majorization-Minimization

Majorization-minimization (MM) is a large class classical approaches for
solving nonconvex optimization problems (see [27] for a recent review) that
encludes classical gradient descent algorithms as well as expectation-
maximizaiton (EM) for solving maximu likelihood estimate (MLE) problems in
statistics [10, 42]. A key idea is local convex relaxation. For an
illustration, suppose we would like to minimize a differentiable function $f$
with $L$-Lipschitz gradient (i.e., $L$-smooth) and the current estimate is
$\boldsymbol{\theta}_{n-1}$. In order to compute the new estimate
$\boldsymbol{\theta}_{n}$, instead of directly minimizing $f$ near
$\boldsymbol{\theta}_{n-1}$, we may minimize its following quadratic expansion

$\displaystyle g_{n}(\boldsymbol{\theta}):\boldsymbol{\theta}\mapsto
f(\boldsymbol{\theta}_{n-1})+\langle\nabla
f(\boldsymbol{\theta}_{n-1}),\,\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rangle+\frac{1}{2\eta}\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert^{2},$
(6)

where $\eta>0$ is a parameter. It is easy to see that $g_{n}$ is minimized at
$\boldsymbol{\theta}_{n}:=\boldsymbol{\theta}_{n-1}-\eta\nabla
f(\boldsymbol{\theta}_{n-1})$, yielding the classical gradient descent update
with stepsize $\eta$. In fact, if one choses $\eta\leq L$, then $g_{n}$
majorizes $f$ (i.e., $g_{n}\geq f$, see Lemma B.1) and the resulting gradient
descent algorithm converges to the global minimum at rate $O(1/k)$ when $f$ is
convex (see, e.g., [6]). In general, MM algorithms proceed by iteratively
minimizing a majorizing surrogate (such as $g_{n}$ above) of the objective
function $f$.

Next, consider using a similar idea for solving the ELM problem (1). There,
instead of a single objective function $f$, we have a sequence of emprical
loss functions $\bar{f}_{n}$ to be minimized online. It would be natural to
first find a majorizing surrogate $g_{n}$ of the one-point loss function
$f_{n}(\cdot):=\ell(\mathbf{x}_{n},\cdot)$ at the current estimate
$\boldsymbol{\theta}_{n}$ and then undergo the same recursive averaging step
as in (1) to compute an averaged surrogate $\bar{g}_{n}$ so that $\bar{g}_{n}$
is a majorizing surrogate of the current objective $\bar{f}_{n}$. We can then
minimize $\bar{g}_{n}$ to compute the next estimate $\boldsymbol{\theta}_{n}$.
This method is called the stochastic marjorization-minimization (SMM) [35],
which we state concisely as follow:

$\displaystyle(\textbf{{SMM}})\quad\textup{Upon arrival of
$\mathbf{x}_{n}$:}\quad\begin{cases}\textup{$g_{n}\leftarrow$ Strongly convex
majorizing surrogate of $f_{n}(\cdot)=\ell(\mathbf{x}_{n},\cdot)$;}\\\
\boldsymbol{\theta}_{n}\in\operatorname*{arg\,min}_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\big{(}\bar{g}_{n}(\boldsymbol{\theta}):=(1-w_{n})\bar{g}_{n-1}(\boldsymbol{\theta})+w_{n}\,g_{n}(\boldsymbol{\theta})\big{)}.\end{cases}$
(7)

SMM has been successful in a variety of settings including online versions of
dictionary learning, matrix factorization, proximal gradient, and DC
programming (see [37, 34, 38]).

A premise of SMM is that $\bar{g}_{t}$ is strongly convex so that it is easy
to find a unique minimizer, which is the case for online matrix factorization
problems [37, 34]. When the surrogate $\bar{g}_{n}$ is only block multi-
convex, that is,
$\bar{g}_{n}:\boldsymbol{\Theta}=\Theta^{(1)}\times\dots\times\Theta^{(m)}\rightarrow[0,\infty)$
and it is convex in each of the $m$ blocks, then finding a global minimum of
$\bar{g}_{n}$ is not easy, if not impossible. Moreover, even when
$\bar{g}_{n}$ is convex, one can still exploit its block multi-convex
structure and use a more efficient coordiante descent method to minimize
$\bar{g}_{n}$. After all, it may be enough to minimize $\bar{g}_{n}$ only
approximately at each $n$ using a cheap coordiante descent method since we are
interested in solving an online problem.

To this end, we propose the following generalization of SMM, where the
surrogate $g_{n}$ functions can only be block-convex (not necessarily jointly
convex) and the averaged surrogate $\bar{g}_{n}$ can only be approximately
minimized (e.g., by using a single round of block coordinate descent) within a
trust region. We call our algorithm Stochastic Block Marjorization-
Minimizaiton (SBMM):

$\displaystyle(\textbf{{SBMM}})\quad\textup{Upon arrival of
$\mathbf{x}_{n}$:}\quad\begin{cases}\textup{$g_{n}\leftarrow$ Block multi-
convex majorizing surrogate of $f_{n}(\cdot)=\ell(\mathbf{x}_{n},\cdot)$;}\\\
\boldsymbol{\theta}_{n}\approx\operatorname*{arg\,min}_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\big{(}\bar{g}_{n}(\boldsymbol{\theta}):=(1-w_{n})\bar{g}_{n-1}(\boldsymbol{\theta})+w_{n}\,g_{n}(\boldsymbol{\theta})\big{)}\\\
\qquad\hskip 4.26773pt\textup{subject to
$\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert\leq
c^{\prime}w_{n}$},\end{cases}$ (8)

where $c^{\prime}\in(0,\infty]$ is a constant. For the definition of block-
convex surrogates, see Definition 3.1.

A motivating example for SBMM is the recent work on online nonnegative
CP111CANDECOMP/PARAFAC for tensor decomposition-dictionary learning in [32],
which encludes online nonnegative tensor CP-decomposition as a special case.
There, the variational surrogate is only convex in each of the $n$ loading
matrices. The method developed in [32] in order to handle a similar issue for
a tensor factorization setting is that, at each round, we only approximately
minimize the surrogate $\bar{g}_{n}$ by a single round of cyclic block
coordinate descent (BCD) in the $n$ loading matrices. This additional layer of
relaxation causes a number of technical difficulties in convergence analysis.
One of the key innovations used to ensure convergence of the algorithm is the
use of “search radius restriction” developed in [33], which can be regarded as
a trust region method [53, 9] with diminishing radius. In this paper, we
generalize and significantly improve this approach ana analysis in [32]. Most
importantly, we will establish a rate of convergence of SBMM in the general
nonconvex, constrained, and Markovian data setting. As a corollary, we obtain
a rate of convergence results for the classical SMM in the Markovian data
case, which has not been available even in the i.i.d. data case.

##  3\. Statement of the algorithm

In order to state the main algorithm, we first define block multi-convex
first-order surrogates. For each $J\subseteq\\{1,\dots,p\\}$, we denote
$\mathbb{R}^{J}$ by the $|J|$-dimensional subspace of $\mathbb{R}^{p}$
generated by the coordinates in $J$. We identify
$\mathbb{R}^{\\{1,\dots,p\\}}$ and $\mathbb{R}^{p}$. We call a subset
$\mathbb{J}\subseteq 2^{\\{1,\dots,p\\}}\setminus\\{\emptyset\\}$ a set of
coordinate blocks if $\bigcup\mathbb{J}=\\{1,\dots,p\\}$. In this case, each
element $J\in\mathbb{J}$ is called a coordinate block (or block coordinate).
Note that two distinct coordinate blocks do not need to be disjoint. For
$\boldsymbol{\Theta}\subseteq\mathbb{R}^{p}$,
$\boldsymbol{\theta}\in\mathbb{R}^{p}$, and a block coordinate
$J\subseteq\\{1,\dots,p\\}$, denote
$\boldsymbol{\Theta}^{J}:=\textup{Proj}_{\mathbb{R}^{J}}(\boldsymbol{\Theta})$
and
$\boldsymbol{\theta}^{J}:=\textup{Proj}_{\mathbb{R}^{J}}(\boldsymbol{\theta})$.
Next, we define the class of surrogate functions we use in this paper:

######  Definition 3.1 (First-order Multi-convex Surrogates).

Fix a set $\mathbb{J}$ of coordinate blocks and parameters $L>0$, $\rho\geq
0$, and $\varepsilon\geq 0$. A function
$g:\mathbb{R}^{p}\rightarrow\mathbb{R}$ is a (first-order) multi-convex
$\varepsilon$-approximate surrogate of $f:\mathbb{R}^{p}\rightarrow\mathbb{R}$
at $\boldsymbol{\theta}^{*}\in\mathbb{R}^{p}$ on block coordinates in
$\mathbb{J}$ if the following hold:

(i)

($\varepsilon$-majorization) $g(\boldsymbol{\theta})+\varepsilon\geq
f(\boldsymbol{\theta})$ for all $\boldsymbol{\theta}\in\mathbb{R}^{p}$;

(ii)

(Smoothness of error) The approximation error $h:=g-f$ is differentiable and
$\nabla h$ is $L$-Lipschitz continuous. Furthermore,
$h(\boldsymbol{\theta}^{*})\leq\varepsilon$ and $\lVert\nabla
h^{*}(\boldsymbol{\theta})\lVert\leq\varepsilon$.

(iii)

(Block multi-convexity) For each $J\in\mathbb{J}$, $g$ is $\rho$-strongly
convex on $\mathbb{R}^{J}$.

We denote
$\operatorname{\textup{{Srg}}}_{L,\rho}^{\mathbb{J}}(f,\boldsymbol{\theta}^{*},\varepsilon)$
for the set of all multi-convex $\varepsilon$-approximate surrogates of $f$ at
$\boldsymbol{\theta}^{*}$ with parameters $L$ and $\mu$ and coordinate blocks
$\mathbb{J}$. If $\\{1,\dots,p\\}\in\mathbb{J}$, then the surrogates are
convex functions and we write
$\operatorname{\textup{{Srg}}}_{L,\rho}(f,\boldsymbol{\theta}^{*},\varepsilon):=\operatorname{\textup{{Srg}}}_{L,\rho}^{\mathbb{J}}(f,\boldsymbol{\theta}^{*},\varepsilon)$,
which consists of $\rho$-strongly convex surrogates of $f$ at
$\boldsymbol{\theta}^{*}$. Note that we allow $\rho=0$, in which case $g_{n}$
is convex in each block. We denote
$\operatorname{\textup{{Srg}}}_{L,\rho}^{\mathbb{J}}(f,\boldsymbol{\theta}^{*}):=\operatorname{\textup{{Srg}}}_{L,\rho}^{\mathbb{J}}(f,\boldsymbol{\theta}^{*},0)$
and call the elements of
$\operatorname{\textup{{Srg}}}_{L,\rho}^{\mathbb{J}}(f,\boldsymbol{\theta}^{*})$
multi-convex surrogates of $f$ at $\boldsymbol{\theta}^{\star}$.

Now we state our main algorithm, Algorithm 1.

Algorithm 1 Stochastic Majorization-Minimization (SMM) with Diminishing Radius

1:Input:
$\boldsymbol{\theta}_{0}\in\boldsymbol{\Theta}\subseteq\mathbb{R}^{p}$
(initial estimate); $N$ (number of iterations); $(w_{n})_{n\geq 1}$, (non-
increasing weights in $(0,1]$); $(\varepsilon_{n})_{n\geq 1}$, (non-increasing
surrogate tolerance in $[0,1]$); $c^{\prime}>0$ (seaerch radius constant);
$\mathbb{J}\in 2^{\\{1,\dots,p\\}}$ (set of coordinate blocks); $L>0$
(surrogate smoothness parameter); $\mu\geq 0$ (surrogate convexity parameter);

2: Initialize the approximate surrogate
$\bar{g}_{0}:\boldsymbol{\theta}\mapsto\frac{\rho}{2}\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{0}\rVert^{2}$;
$\bar{\boldsymbol{\theta}}_{0}=\boldsymbol{\theta}_{0}$;
$\hat{\boldsymbol{\theta}}=\boldsymbol{\theta}_{0}$;

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

4: sample a training point $\mathbf{x}_{n}\in\mathfrak{X}$; define
$f_{n}:\boldsymbol{\theta}\mapsto\ell(\mathbf{x}_{n},\boldsymbol{\theta})$;

5: choose a surrogate function
$g_{n}\in\operatorname{\textup{{Srg}}}_{L,\mu}^{\mathbb{J}}(f_{n},\boldsymbol{\theta}_{n-1},\varepsilon_{n})$;

6: update the average surrogate:
$\bar{g}_{n}=(1-w_{n})\bar{g}_{n-1}+w_{n}g_{n}$;

7: use Algorithm 2 to update the current estimate $\boldsymbol{\theta}_{n-1}$
to $\boldsymbol{\theta}_{n}$ such that
$\displaystyle\qquad\qquad\boldsymbol{\theta}_{n}\approx\operatorname*{arg\,min}_{\boldsymbol{\theta}\in\boldsymbol{\Theta},\,\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert\leq
c^{\prime}w_{n}}\bar{g}_{n}(\boldsymbol{\theta})\qquad\left(\triangleright\,\textup{Inexact
surrogate minimization within trust region}\right)$ (9)

8: end for

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

The algorithm iteratively execute the following: 1) Sample a new data point
$\mathbf{x}_{n}\in\mathfrak{X}$ using an a priori sampling algorithm (e.g.,
MCMC); 2) Choose a new $m$-block multi-convex first-order surrogate $g_{n}$ of
the loss function
$\boldsymbol{\theta}\mapsto\ell(\mathbf{x}_{n},\boldsymbol{\theta})$; 3)
Update and aggregate surrogate $\bar{g}_{n}$ by taking an weighted average of
$\bar{g}_{n-1}$ and $g_{n}$; 4) Perform a single round of block coordinate
descent on $\bar{g}_{n}$ over the $m$ blocks with search radius restricted by
$r_{n}$, where $(r_{n})_{n\geq 1}$ is a sequence of prescribed radii (e.g.,
$r_{n}=1/\sqrt{n}$). There are three options in the algorithm, according to
which the final estimate is appropriately averaged. The averaging scheme is
used when $f$ is convex or strongly convex. Note that Algorithm 1 reduces to
SMM in [35] when $m=1$, $c^{\prime}=\infty$, and the surrogates $g_{n}$ are
$\rho$-strongly convex for some $\rho>0$.

Algorithm 2 Block Surrogate Minimization

1:Input: all input of Algorithm 1;
$\boldsymbol{\theta}_{n-1}\in\boldsymbol{\Theta}$ (current estimate);
$m\in\mathbb{N}$ (number of sub-iterations);

2: For $i=1,\cdots,m$ do:

3: choose a block coordinate $J=J_{i}(n)\in\mathbb{J}$ independently of
everything else;
$\displaystyle\qquad\boldsymbol{\Theta}^{J}_{n}\leftarrow\left\\{\boldsymbol{\theta}\in\boldsymbol{\Theta}\,\bigg{|}\,\boldsymbol{\theta}^{J^{c}}=\boldsymbol{\theta}_{n-1}^{J^{c}},\,\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n-1}\rVert\leq
c^{\prime}w_{n}/m\right\\}\quad\left(\triangleright\begin{matrix}\textup{Freeze
coordinates in $J^{c}$ and }\\\ \textup{restrict the search radius by
$c^{\prime}w_{n}/m$}\end{matrix}\right);$ (10)
$\displaystyle\qquad\boldsymbol{\theta}_{n-1}^{(i)}\in\operatorname*{arg\,min}_{\boldsymbol{\theta}\in\boldsymbol{\Theta}^{J}_{n}}\,\bar{g}_{n}\left(\boldsymbol{\theta}\right)\qquad\left(\triangleright\begin{matrix}\textup{Update
$\boldsymbol{\theta}_{n-1}$ on coordinates in $J$ }\\\ \textup{while holding
the rest of coordinates in $J^{c}$}\end{matrix}\right);$ (11)

4: $\boldsymbol{\theta}_{n}\leftarrow\boldsymbol{\theta}_{n-1}^{(m)}$;

5:output: $\boldsymbol{\theta}_{n}\in\boldsymbol{\Theta}$

Algorithm 1 combines two key elements: stochastic majorization-minimization
(SMM) [35] and block coordinate descent with diminishing radius (BCD-DR) [33].
Since when the surrogates $g_{n}$ (and hence their weighted averages
$\bar{g}_{n}$) are block multi-convex, it is infeasible to minimize
$\bar{g}_{n}$ over $\boldsymbol{\Theta}$ as in the standard SMM [35]. In order
to circumvent this issue, one could try to perform BCD on $\bar{g}_{n}$.
However, this results in sub-optimal estimates in each iteration, causing a
number of difficulties in convergence analysis. Moreover, global convergence
of BCD to stationary points is not guaranteed in general even for the
deterministic cases (see, e.g., [23]), and such a guarantee is known only with
some additional regularity conditions [16, 17, 3]. There are other popular
strategies of using proximal [17] or prox-linear [52] modifications of BCD to
improve convergence properties. While these methods only ensure square-
summability
$\sum_{t=1}^{\infty}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert^{2}<\infty$
of changes, we find it critical for our stochastic analysis that we are able
to control the change
$\lVert\boldsymbol{\theta}_{t}-\boldsymbol{\theta}_{n-1}\rVert$ we make to the
loading matrices in each iteration. This motivates to use BCD with diminishing
radius in [33].

Each of the block minimization problems in (11) is a convex optimization
problem over convex constriant $\boldsymbol{\Theta}_{n}^{(i)}$, so it can be
easily solved by a number of known algorithms (e.g., LARS [13], LASSO [47],
and feature-sign search [26]).

##  4\. Main results

###  4.1. Optimality conditions and convergence measures

In this subsection, we introduce some notions on optimality conditions and
related quantities. Here we denote $f$ to be a generic objective function
$\boldsymbol{\Theta}\rightarrow\mathbb{R}$, but elsewhere $f$ will denote the
expected loss function in (4) unless otherwise menteiond.

Recall that we say $\boldsymbol{\theta}^{*}\in\boldsymbol{\Theta}$ is a
stationary point of $f$ over $\boldsymbol{\Theta}$ if

$\displaystyle\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\,\langle\nabla
f(\boldsymbol{\theta}^{*}),\,\boldsymbol{\theta}-\boldsymbol{\theta}^{*}\rangle\geq
0,$ (12)

where $\langle\cdot,\,\cdot\rangle$ denotes the dot project on
$\mathbb{R}^{p}\supseteq\boldsymbol{\Theta}$. This is equivalent to saying
that $-\nabla f(\boldsymbol{\theta}^{*})$ is in the normal cone of
$\boldsymbol{\Theta}$ at $\boldsymbol{\theta}^{*}$. If
$\boldsymbol{\theta}^{*}$ is in the interior of $\boldsymbol{\Theta}$, then it
implies $\lVert\nabla f(\boldsymbol{\theta}^{*})\rVert=0$. For iterative
algorithms, such a first-order optimality condition may hardly be satisfied
exactly in a finite number of iterations, so it is more important to know how
the worst-case number of iterations required to achieve an
$\varepsilon$-approximate solution scales with the desired precision
$\varepsilon$. More precisely, we say
$\boldsymbol{\theta}^{*}\in\boldsymbol{\Theta}$ is an $\varepsilon$-approxiate
stationary point of $f$ over $\boldsymbol{\Theta}$ if

$\displaystyle-\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\,\left\langle\nabla
f(\boldsymbol{\theta}^{*}),\,\frac{(\boldsymbol{\theta}-\boldsymbol{\theta}^{*})}{\lVert\boldsymbol{\theta}-\boldsymbol{\theta}^{*}\rVert}\right\rangle\leq\sqrt{\varepsilon}.$
(13)

This notion of $\varepsilon$-approximate solution is consistant with the
corresponding notion for unconstrained problems. Indeed, if
$\boldsymbol{\theta}^{*}$ is an interior point of $\boldsymbol{\Theta}$, then
(13) reduces to $\lVert\nabla
f(\boldsymbol{\theta}^{*})\rVert^{2}\leq\varepsilon$. See also [41, Def. 1].

Next, for each $\varepsilon>0$ we define the iteration complexity
$N_{\varepsilon}$ of an algorithm for minimizing $f$ over
$\boldsymbol{\Theta}$ with initialization
$\boldsymbol{\theta}_{0}\in\boldsymbol{\Theta}$ as

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

where $(\boldsymbol{\theta}_{n})_{n\geq 0}$ is a sequence of estimates
produced by the algorithm with initial estimate $\boldsymbol{\theta}_{0}$. An
upper bound on $N_{\varepsilon}$ can be regarded as the worst-case bound on
the number of iterations for an algorithm to achieve an
$\varepsilon$-approximate solution.

###  4.2. Assumptions

In this subsection, we state all assumptions we use for establishing the main
results.

######  (A1).

The loss function
$\ell:\mathfrak{X}\times\boldsymbol{\Theta}\rightarrow[0,\infty)$ is bounded,
$R$-Lipscthiz continuous for some $R>0$, and differentiable with Lipscthiz
gradients.

######  (A2).

The observed data points $\mathbf{x}_{n}\in\mathfrak{X}$ are given by
$\mathbf{x}_{n}=\varphi(Y_{n})$, where $Y_{n}$ is Markov chain on state space
$\Omega$ and $\varphi:\Omega\rightarrow\mathbb{R}^{p}$ is a function with
compact image. Furthermore, $Y_{n}$ has a unique stationary distribution $\pi$
and satisfies exponential mixing with rate $\lambda\in[0,1)$:

$\displaystyle\sup_{\mathbf{y}\in\Omega}\,\lVert
P^{m}(\mathbf{y},\cdot)-\pi\rVert_{TV}\leq\lambda^{m}.$ (15)

######  (A3).

The constraint set $\boldsymbol{\Theta}$ is a compact and convex subset of
$\mathbb{R}^{p}$. Let $\mathbb{J}\subseteq
2^{\\{1,\dots,p\\}}\setminus\\{\emptyset\\}$ denote the set of coordinate
blocks used in Algorithm 1. Then for each block coodinate $J\in\mathbb{J}$,
$\boldsymbol{\Theta}^{J}$ contains an open ball in $\mathbb{R}^{J}$.

######  (A4).

The weights $w_{n}\in(0,1]$ are non-increasing and satisfy
$w_{n}^{-1}-w_{n-1}^{-1}\leq 1$ for all sufficiently large $n\geq 1$,
$\sum_{n=1}^{\infty}w_{n}=\infty$, and $w_{n}\sqrt{n}=O(1/(\log
n)^{1+\varepsilon})$ for some $\varepsilon>0$. Also, there exists a constant
$C>0$ such that for $a_{n}=\lfloor C\log n\rfloor$, the following hold:

$\displaystyle\sum_{n=1}^{\infty}w_{n}^{2}\sum_{k=1}^{n+1}\mathbb{E}[\varepsilon_{k}]<\infty,\qquad\sum_{n=1}^{\infty}a_{n}w_{n}w_{n-a_{n}}<\infty,\qquad\sum_{n=1}^{\infty}w_{n}\lambda^{a_{n}}<\infty,$
(16) $\displaystyle\hskip 56.9055pt\qquad\left(\textup{optionally,
}\sum_{n=1}^{\infty}w_{n}w_{n-a_{n}}\sqrt{n}<\infty\right),$ (17)

where $(\varepsilon_{n})_{n\geq 1}$ is given in Algorithm 1 and
$\lambda\in(0,1]$ is given in (A2).

######  (A5).

The joint distribution $(J_{1}(n),\dots,J_{m}(n))$ of coordinate blocks chosen
in Algorithm 2 at iteration $n$ of Algorithm 1 does not depend on $n$.
Furthermore, the expected number of each coordinate in $\\{1,\dots,p\\}$
appearing in all coordinate blocks $J_{1}(n),\dots,J_{m}(n)$ is constant.

######  (A6).

The averaged surrogates $\bar{g}_{n}$ are parameterized by some variable
$\kappa_{n}$ in some compact set $\mathcal{K}$. That is, there exists a
function $\bar{g}:\mathcal{K}\times\boldsymbol{\Theta}\rightarrow[0,\infty)$
such that
$\bar{g}_{n}(\boldsymbol{\theta})=\bar{g}(\kappa_{n},\boldsymbol{\theta},\varepsilon_{n})$
for some $\kappa_{n}\in\mathcal{K}$. Furthermore, $\bar{g}$ is Lipschitz in
the first coordinate.

Assumption (A4) states that the sequence of weights $w_{n}\in(0,1]$ we use to
recursively define the empirical loss (1) and surrogate loss (8) does not
decay too fast so that $\sum_{n=1}^{\infty}w_{n}=\infty$ but decay fast enough
so that $\sum_{n=1}^{\infty}w_{n}^{2}<\infty$. This is analogous to
requirements for stepsizes in stochastic gradient descent algorithms, where
the stepsizes are usually required to be non-summable but square-summable
(see, e.g., [46]). Note that our general results do not require the stronger
assumption $\sum_{n=1}^{\infty}w_{n}^{2}\sqrt{n}<\infty$, which is standard in
the literature [37, 35, 38, 30, 32]. Also, the condition
$w_{n}^{-1}-w_{n-1}^{-1}\leq 1$ for all suficiently large $n$ is equivalent to
saying the recursively defined weights $w^{n}_{k}$ in (2) are non-decreasing
in $k$ for all sufficiently large $k$, which is required to use Lemma B.4. We
also remark that (A4) is implied by the following simpler condition:

######  (A4’).

The sequence of non-increasing weights $w_{n}\in(0,1]$ in Algorithm 1 satisfy
either $w_{n}=n^{-1}$ for $n\geq 1$ or $w_{n}=\Theta(n^{-\beta}(\log
n)^{-\delta})$ for some $\delta>1$ and $\beta\in[1/2,1)$ (Optionally,
$\beta\in[3/4,1)$). Also, the sequence $(\varepsilon_{n})_{n\geq 1}$ in
Algorithm 1 satisfies
$\sum_{n=1}^{\infty}w_{n+1}^{2}\sum_{k=1}^{n+1}\mathbb{E}[\varepsilon_{k+1}]<\infty$.

If we consider the weight $w_{n}=n^{-\beta}/(\log n)^{\delta}$ for some
$\beta\in[1/2,1)$ and $\delta>1$ as in (A4’), we have

$\displaystyle\left(\sum_{k=1}^{n}w_{n}\right)^{-1}=O(n^{\beta-1}(\log
n)^{\delta}).$ (18)

Hence the above bound is optimized when $\beta=1/2$ for each fixed $\delta>1$.

###  4.3. Statement of main results

Now we state our first main result, Theorem 4.1, which states that Algorithm 1
converges globally (w.r.t. initialization) to the set of stationary points of
both the empirical loss $\bar{f}_{n}$ and the expected loss $f$. Moreover, it
also states that the averaged surrogate $\bar{g}_{n}$ is asymptotically an
accurate approximation of the empirical and the expected loss functions at
$\boldsymbol{\theta}_{n}$ both in the function values and gradients.

######  Theorem 4.1 (Global Convergence).

Let $(\boldsymbol{\theta}_{n})_{n\geq 1}$ be an output of Algorithm 1 of one
of the following instances:

C1.

(Strongly convex surrogates without radius restriction) $c^{\prime}=\infty$
and
$g_{n}\in\operatorname{\textup{{Srg}}}_{L,\rho}(f_{n},\boldsymbol{\theta}_{n-1})$
for $n\geq 1$ for some $L,\rho>0$;

C2.

(Block multi-convex surrogates with radius restriction) $c^{\prime}<\infty$
and
$g_{n}\in\operatorname{\textup{{Srg}}}_{L,0}^{\mathbb{L}}(f_{n},\boldsymbol{\theta}_{n-1})$
for $n\geq 1$ for some set $\mathbb{J}$ of coordinate blocks and for some
$L>0$. In this case, assume (A6).

Assume (A1)-(A5) hold. Then the following hold:

(i)

(Empirical Loss Minimization) Both
$|\bar{g}_{n}(\boldsymbol{\theta}_{n})-\bar{f}_{n}(\boldsymbol{\theta}_{n})|$
and
$\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})-\nabla\bar{f}_{n}(\boldsymbol{\theta}_{n})\rVert$
converge to zero almost surely. Furthermore, $\boldsymbol{\theta}_{n}$ is
asymptotically a stationary point of $\bar{f}_{n}$ over $\boldsymbol{\Theta}$
almost surely.

(ii)

(Expected Loss Minimization) All of
$|\bar{g}_{n}(\boldsymbol{\theta}_{n})-f(\boldsymbol{\theta}_{n})|$,
$\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})-\nabla
f(\boldsymbol{\theta}_{n})\rVert$ and
$\lVert\mathbb{E}[\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})]-\nabla
f(\boldsymbol{\theta}_{n})\rVert$ converge to zero alsmot surely. Furthermore,
$\boldsymbol{\theta}_{n}$ converges to the set of a stationary points of $f$
over $\boldsymbol{\Theta}$ almost surely.

We note that standard convergence results in the literature of SMM [37, 34,
30] asserts that the SMM algorithm converges globally to the set of the
stationary points of the expected loss function $f$ almost surely, which is
recovered by Theorem 4.1 (ii) for the first case of strongly convex surrogates
without radius restriction. Note that we establish Theorem 4.1 without the
optional condition in (A4), which essentially states that
$\sum_{n=1}^{\infty}w_{n}^{2}\sqrt{n}<\infty$. Such condition was standard in
the literature (see, e.g., [35, Assumption (E)]). The same statement for the
second case of block multi-convex surrogates with radius restriction was
recently obtained in [32] for the contex of online CP-dictionary learning.

The convergence of gradient norms
$\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})-\nabla\bar{f}_{n}(\boldsymbol{\theta}_{n})\rVert$
and $\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})-\nabla
f(\boldsymbol{\theta}_{n})\rVert$ are new and the asymptotic stationarity for
the empirical loss functions in Theorem 4.1 (i) has not been elaborated very
much in the aforementioned literature, since the main focus of using SMM was
to solve the expected loss minimizaiton. It is worth noting that the
hypothesis for the expected loss minimization stated in Theorem 4.1 (ii) is a
bit stronger than that for the empirical loss minimization stated in Theorem
4.1 (i). This is an indication that SBMM (and hence SMM) is generically more
suited to solve the empirical loss minimization than the expected loss
minimization, which we will elaborate on the forthcoming results.

In the following two results, Theorems 4.2 and 4.3, we establish bounds on the
rate of convergence of surrogate gaps as well as approximate optimality.
Below, we denote by $\mathbb{E}[\cdot]$ the coditional expectation
$\mathbb{E}[\cdot\,|\,\mathcal{F}_{0}]$ with respect to the time-0 information
that contains the initial estimate $\boldsymbol{\theta}_{0}$.

######  Theorem 4.2 (Rate of Convergence of Surrogate Gaps and Variation).

Let $(\boldsymbol{\theta}_{n})_{n\geq 1}$ be an output of Algorithm 1. Make
the same assumption as in Theorem 4.1. Then the following hold:

(i)

(Empirical Loss Minimization) Asymtotically almost surely,

$\displaystyle\min_{1\leq k\leq
n}\begin{pmatrix}\left|\bar{g}_{k}(\boldsymbol{\theta}_{k})-\bar{f}_{k}(\boldsymbol{\theta}_{k})\right|+\lVert\nabla\bar{g}_{k}(\boldsymbol{\theta}_{k})-\nabla\bar{f}_{n}(\boldsymbol{\theta}_{k})\rVert^{2}\end{pmatrix}$
$\displaystyle=O\left(\left(\sum_{k=1}^{n}w_{k}\right)^{-1}\right).$ (19)

Furthermore, if the optional condition in (A4) holds, then asymptotically
almost surely,

$\displaystyle\min_{1\leq k\leq
n}\begin{pmatrix}\,\,\left|\mathbb{E}[\bar{g}_{k}(\boldsymbol{\theta}_{k})]-\bar{f}_{k}(\boldsymbol{\theta}_{k})\right|+\lVert\mathbb{E}[\nabla\bar{g}_{k}(\boldsymbol{\theta}_{k})]-\nabla\bar{f}_{k}(\boldsymbol{\theta}_{k})\rVert^{2}\end{pmatrix}$
$\displaystyle=O\left(\left(\sum_{k=1}^{n}w_{k}\right)^{-1}\right).$ (20)

(ii)

(Expected Loss Minimization) Asymptotically almost surely,

$\displaystyle\min_{1\leq k\leq
n}\begin{pmatrix}\,\,\left|\mathbb{E}[\bar{g}_{k}(\boldsymbol{\theta}_{k})]-f(\boldsymbol{\theta}_{k})\right|+\lVert\mathbb{E}[\nabla\bar{g}_{k}(\boldsymbol{\theta}_{k})]-\nabla
f(\boldsymbol{\theta}_{k})\rVert^{2}\end{pmatrix}$
$\displaystyle=O\left(\left(\sum_{k=1}^{n}w_{k}\right)^{-1}\right).$ (21)

Furthermore, if the optional condition in (A4) holds, then asymptotically
almost surely,

$\displaystyle\min_{1\leq k\leq
n}\left(\left|\bar{g}_{k}(\boldsymbol{\theta}_{k})-f(\boldsymbol{\theta}_{k})\right|+\lVert\nabla\bar{g}_{k}(\boldsymbol{\theta}_{k})-\nabla
f(\boldsymbol{\theta}_{k})\rVert^{2}\right)=O\left(\left(\sum_{k=1}^{n}w_{k}\right)^{-1}\right).$
(22)

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

Let $(\boldsymbol{\theta}_{n})_{n\geq 1}$ be an output of Algorithm 1. Make
the same assumption as in Theorem 4.1. Then the following hold:

(i)

(Surrogate and Empirical Loss Stationarity) Asymtotically almost surely,

$\displaystyle\min_{1\leq k\leq
n}\,\,\left[-\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\langle\nabla\bar{g}_{k}(\boldsymbol{\theta}_{k}),\,\frac{(\boldsymbol{\theta}-\boldsymbol{\theta}_{k})}{\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{k}\rVert}\right\rangle\right]=O\left(\left(\sum_{k=1}^{n}w_{k}\right)^{-1}\right).$
(23) $\displaystyle\min_{1\leq k\leq
n}\,\,\left[-\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\langle\nabla\bar{f}_{k}(\boldsymbol{\theta}_{k}),\,\frac{(\boldsymbol{\theta}-\boldsymbol{\theta}_{k})}{\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{k}\rVert}\right\rangle\right]=O\left(\left(\sum_{k=1}^{n}w_{k}\right)^{-1/2}\right).$
(24)

(ii)

(Expected Loss Stationarity) Further assume the optional condition in (A4)
holds. Then asymtotically almost surely,

$\displaystyle\min_{1\leq k\leq
n}\,\,\left[-\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\langle\nabla
f(\boldsymbol{\theta}_{k}),\,\frac{(\boldsymbol{\theta}-\boldsymbol{\theta}_{k})}{\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{k}\rVert}\right\rangle\right]=O\left(\left(\sum_{k=1}^{n}w_{k}\right)^{-1/2}\right).$
(25)

To our best knowledge, the rate of convergence results in Theorems 4.2 and 4.3
are entirely new for SMM-type algorithms even under the i.i.d. data setting.
Only almost sure convergence to stationary points were known before [37, 35,
54, 30, 32] even under the i.i.d. data assumption with strongly conex
surrogates. Moreover, (20) and (21) in Theorem 4.2 give bounds on the
variation of the random objective values
$\bar{f}_{n}(\boldsymbol{\theta}_{n})$ and $f(\boldsymbol{\theta}_{n})$
against the deterministic quantity
$\mathbb{E}[\bar{g}_{n}(\boldsymbol{\theta}_{n})]$ with respect to the
randomness of data samples $\mathbf{x}_{1},\dots,\mathbf{x}_{n}$ and of
Algorithm 1 (e.g., possibly randomized choice of blocks in Algorithm 2).

Next, we state a corollary of Theorems 4.2 and 4.3, which specializes these
results to a more familar setting of unconstrained nonconvex optimization.

######  Corollary 4.4.

Let $(\boldsymbol{\theta}_{n})_{n\geq 1}$ be an output of Algorithm 1. Make
the same assumption as in Theorem 4.1. Further assume that
$\boldsymbol{\theta}_{n}$ is in the interior of $\boldsymbol{\Theta}$ for all
$n\geq 1$. Then asymptotically almost surely,

$\displaystyle\min_{1\leq k\leq
n}\left\lVert\nabla\bar{g}_{k}(\boldsymbol{\theta}_{k})\right\rVert^{2}=O\left(\left(\sum_{k=1}^{n}w_{k}\right)^{-2}\right),\quad\min_{1\leq
k\leq
n}\left\lVert\nabla\bar{f}_{k}(\boldsymbol{\theta}_{k})\right\rVert^{2}=O\left(\left(\sum_{k=1}^{n}w_{k}\right)^{-1}\right),$
(26) $\displaystyle\hskip 76.82234pt\quad\min_{1\leq k\leq n}\left\lVert\nabla
f(\boldsymbol{\theta}_{k})\right\rVert^{2}=O\left(\left(\sum_{k=1}^{n}w_{k}\right)^{-1}\right).$
(27)

Below we unpack the above results and give some remarks. For a more direct
interpretation of our results, we will consider the weight
$w_{n}=n^{-\beta}/(\log n)^{1+\varepsilon}$ for some $\beta\in[1/2,1)$ and
$\varepsilon>0$ as in (A4’). Then the bound in (18) is optimized at
$\beta=1/2$ for each fixed $\varepsilon>0$ to

$\displaystyle\left(\sum_{k=1}^{n}w_{n}\right)^{-1}=O(n^{-1/2}(\log
n)^{1+\varepsilon}).$ (28)

Suppose all iterates $\boldsymbol{\theta}_{n}$ are in the interior of
$\boldsymbol{\Theta}$. This would be a reasonable assumption when all
stationary points of the averaged surrogates $\bar{g}_{n}$ are in the interior
of $\boldsymbol{\Theta}$. Then by Corollary 4.4,

$\displaystyle\min_{1\leq k\leq
n}\left\lVert\nabla\bar{g}_{k}(\boldsymbol{\theta}_{k})\right\rVert^{2}=O\left(\frac{(\log
n)^{2+2\varepsilon}}{n}\right),\quad\min_{1\leq k\leq
n}\left\lVert\nabla\bar{f}_{k}(\boldsymbol{\theta}_{k})\right\rVert^{2}=O\left(\frac{(\log
n)^{1+\varepsilon}}{\sqrt{n}}\right),$ (29) $\displaystyle\hskip
76.82234pt\quad\min_{1\leq k\leq n}\left\lVert\nabla
f(\boldsymbol{\theta}_{k})\right\rVert^{2}=O\left(\frac{(\log
n)^{1+\varepsilon}}{\sqrt{n}}\right).$ (30)

The last asymptotic expresses a well-known rate of convergence bound for
nonconvex and unconstrained SGD [46, 50] in the weaker form of $L_{1}$. A
similar rate of convergence for nonconvex and constrained projected SGD is
also known in [11], although a different measure (using Moreu envelope) of the
rate of convergence for the constrained nonconvex problems was used.

For the general case when $\bar{g}_{n}$ may have its stationary points at the
boundary of $\boldsymbol{\Theta}$, we cannot measure first-order optimality
simply by the gradient norms as above. In this case, for the empirical loss
minimization, (24) in Theorem 4.3 states that

$\displaystyle\min_{1\leq k\leq
n}\,\,\left[-\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\langle\nabla\bar{f}_{k}(\boldsymbol{\theta}_{k}),\,\frac{(\boldsymbol{\theta}-\boldsymbol{\theta}_{k})}{\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{k}\rVert}\right\rangle\right]^{2}$
$\displaystyle=O\left(\frac{(\log n)^{1+\varepsilon}}{\sqrt{n}}\right),$ (31)

which obtains the same asymptotic rate of convergence as in the interior
stationary points case above. As for the expected loss minimization, we need
to assume the strong er condition
$\sum_{n=1}^{\infty}w_{n}^{2}\sqrt{n}<\infty$ (i.e., the optional condition in
(A4)). In this case, the optimal value of $\beta$ is $3/4$, in which case we
obtain

$\displaystyle\min_{1\leq k\leq
n}\,\,\left[-\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\langle\nabla
f(\boldsymbol{\theta}_{k}),\,\frac{(\boldsymbol{\theta}-\boldsymbol{\theta}_{k})}{\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{k}\rVert}\right\rangle\right]=O\left(\frac{(\log
n)^{1+\varepsilon}}{n^{1/4}}\right).$ (32)

Hence we lose a factor of $n^{-1/4}$ from the optimal bound in (26). This
discussion gives a more quantitative indication that SBMM (and hence SMM) is
generically more suited to solve the empirical loss minimization than the
expected loss minimization for constrained nonconvex objectives.

Next, we remark on the effect of using inexact minimization of block multi-
convex surrogates on convergence rates. All results stated above entail two
cases: (C1) strongly convex surrogates are exactly minimized without search
radius restriction $(c^{\prime}=\infty)$; and (C2) block multi-convex
surrogates are cyclically minimized with search radius restriction
($c^{\prime}<\infty$). First in case C1, each $\boldsymbol{\theta}_{n}$ is the
minimizer of the strongly convex surrogates $\bar{g}_{n}$ over
$\boldsymbol{\Theta}$ so that $-\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})$ is
in the normal cone of $\boldsymbol{\Theta}$ at $\boldsymbol{\theta}_{n}$.
Hence Theorem 4.2 can be directly used to obtain (24) and (25). In particular,
$\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})\rVert$ is identically zero
when there is no boundary statinary points of $\bar{g}_{n}$.

On the other hand, in case C2, the block multi-convex surrogates $\bar{g}_{n}$
are only approximately minimized with search radius restriction (see Algorithm
1), so $-\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})$ is not necessarily in the
normal cone of $\boldsymbol{\Theta}$ at $\boldsymbol{\theta}_{n}$. However,
(23) shows that $-\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})$ should converge
to some vector the normal cone of $\boldsymbol{\Theta}$ at
$\boldsymbol{\theta}_{n}$ at rate $O((\sum_{k=1}^{n}w_{k})^{-1})$, and by
(25), the bound on the rate of convergence to stationary points of $f$ over
$\boldsymbol{\Theta}$ is slower at order $O((\sum_{k=1}^{n}w_{k})^{-1/2})$.
Hence, the asymptotic rate of convergence is unchanged when we relax strongly
convex surrogates to block multi-convex surrogates.

Lastly, we state a corollary of the earlier results on the iteration
complexity of Algorithm 1. Corollary 4.5 gives a “worst-case” rate of
convergence of Algorithm 1 until reaching an $\varepsilon$-stationary points
of the objective functions.

######  Corollary 4.5 (Iteration Complexity).

Let $(\boldsymbol{\theta}_{n})_{n\geq 1}$ be an output of Algorithm 1. Make
the same assumption as in Theorem 4.1. Then the following hold:

(i)

(Empirical Loss Minimization) Suppose $w_{n}=n^{-1/2}(\log n)^{\delta}$ for
some $\delta>1$. Then we have the following worst-case iteration complexity
for Algorithm 1 with the empricial loss objective:

$\displaystyle
N_{\varepsilon}(\bar{f}_{n})=O(\varepsilon^{-2}(\log\varepsilon^{-1})^{\delta}).$
(33)

(ii)

(Expected Loss Minimization) Suppose $w_{n}=n^{-3/4}(\log n)^{\delta}$ for
some $\delta>1$. Then we have the following worst-case iteration complexity
for Algorithm 1 with the expected loss objective:

$\displaystyle
N_{\varepsilon}(f)=O(\varepsilon^{-4}(\log\varepsilon^{-1})^{\delta}).$ (34)

Furthermore, suppose $\boldsymbol{\theta}_{n}$ is in the interior of
$\boldsymbol{\Theta}$ for $n\geq 1$. Then we may choose $w_{n}=n^{-1/2}(\log
n)^{\delta}$ for some $\delta>1$ and obtain

$\displaystyle
N_{\varepsilon}(f)=O(\varepsilon^{-2}(\log\varepsilon^{-1})^{\delta}).$ (35)

##  5\. Examples and Applications

###  5.1. Example of Block-convex Surrogate functions

In this subsection, we list some examples of block-convex surrogate functions,
which include the usual convex surrogate functions in the literature (see,
[34, 35]).

######  Example 5.1 (Lipschitz Gradient Surrogates).

If $\nabla f$ is $L$-Lipschitz, then the following function $g$ belongs to
$\operatorname{\textup{{Srg}}}_{2L,L}(f,\boldsymbol{\theta}^{\star})$:

$\displaystyle g:\boldsymbol{\theta}\mapsto
f(\boldsymbol{\theta}^{\star})+\left\langle\nabla
f(\boldsymbol{\theta}^{\star}),\,\boldsymbol{\theta}-\boldsymbol{\theta}^{\star}\right\rangle+\frac{L}{2}\lVert\boldsymbol{\theta}-\boldsymbol{\theta}^{\star}\rVert^{2}.$
(36)

If in addition $f$ is convex, then
$g\in\operatorname{\textup{{Srg}}}_{L,L}(f_{n},\boldsymbol{\theta}^{\star})$;
If $f$ is $\mu$-strongly convex, then
$g\in\operatorname{\textup{{Srg}}}_{L-\mu,L}(f_{n},\boldsymbol{\theta}^{\star})$.
$\blacktriangle$

######  Example 5.2 (Proximal Gradient Surrogates).

Suppose $f=f_{1}+f_{2}$ where $f_{1}$ is differentiable with $L$-Lipschitz
gradient and $f_{2}$ is convex over $\boldsymbol{\Theta}$. Then the following
function $g$ belongs to
$\operatorname{\textup{{Srg}}}_{2L,L}(f,\boldsymbol{\theta}^{\star})$:

$\displaystyle g:\boldsymbol{\theta}\mapsto
f_{1}(\boldsymbol{\theta}^{\star})+\left\langle\nabla
f_{1}(\boldsymbol{\theta}^{\star}),\,\boldsymbol{\theta}-\boldsymbol{\theta}^{\star}\right\rangle+\frac{L}{2}\lVert\boldsymbol{\theta}-\boldsymbol{\theta}^{\star}\rVert^{2}+f_{2}(\boldsymbol{\theta}).$
(37)

Minimizing $g$ over $\boldsymbol{\Theta}$ amounts to performing a proximal
gradient step in [7, 41]. $\blacktriangle$

######  Example 5.3 (DC programming surrogates).

Suppose $f=f_{1}+f_{2}$ where $f_{1}$ is convex and $f_{2}$ is concave and
differentiable with $L_{2}$-Lipschitz gradient over $\boldsymbol{\Theta}$. One
can also write $f=f_{1}-(-f_{2})$ which is the difference of convex (DC)
functions $f_{1}$ and $-f_{2}$. Then the following function $g$ belongs to
$\operatorname{\textup{{Srg}}}_{2L,L}(f,\boldsymbol{\theta}^{\star})$:

$\displaystyle g:\boldsymbol{\theta}\mapsto
f_{1}(\boldsymbol{\theta})+f_{2}(\boldsymbol{\theta}^{\star})+\left\langle\nabla
f_{2}(\boldsymbol{\theta}^{\star}),\,\boldsymbol{\theta}-\boldsymbol{\theta}^{\star}\right\rangle.$
(38)

Such surrogates are important in DC-programming (see, e.g., [20]).

In fact, our method allows $f_{1}$ to be only multi-convex with respect to the
coordinate blocks in $\mathbb{J}$, that is, $f_{1}$ is convex on
$\textup{Proj}_{\mathbb{R}^{J}}(\boldsymbol{\Theta})$ for each coodinate block
$J\in\mathbb{J}$. In this case, $g$ above is a multi-convex surrogate of $f$
at $\boldsymbol{\theta}^{\star}$ and belongs to
$\operatorname{\textup{{Srg}}}_{2L,L}^{\mathbb{J}}(f,\boldsymbol{\theta}^{\star})$.
$\blacktriangle$

######  Example 5.4 (Convex Variational Surrogates).

Let $f:\mathbb{R}^{p}\times\mathbb{R}^{q}\rightarrow\mathbb{R}$ be a two-block
multi-convex function and let $\Theta_{1}\subseteq\mathbb{R}^{p}$ and
$\Theta_{2}\subseteq\mathbb{R}^{q}$ be two convex sets. Define a function
$f_{*}:\boldsymbol{\theta}\mapsto\inf_{H\in\Theta_{2}}f(\boldsymbol{\theta},H)$.
Then for each $\boldsymbol{\theta}^{\star}\in\Theta$, the following function

$\displaystyle g:\boldsymbol{\theta}\mapsto
f(\boldsymbol{\theta},H^{\star}),\quad
H^{\star}\in\operatorname*{arg\,min}_{H\in\boldsymbol{\Theta}_{2}}f(\boldsymbol{\theta}^{\star},H)$
(39)

is convex over $\Theta_{1}$ and satisfies $g\geq f$ and
$g(\boldsymbol{\theta}^{\star})=f(\boldsymbol{\theta}^{\star})$. Further
assume that

(i)

$\boldsymbol{\theta}\mapsto f(\boldsymbol{\theta},H)$ is differentiable for
all $H\in\Theta_{2}$ and
$\boldsymbol{\theta}\mapsto\nabla_{\boldsymbol{\theta}}f(\boldsymbol{\theta},H)$
is $L^{\prime}$-Lipschitz for all $H\in\Theta_{2}$ ;

(ii)

$H\mapsto\nabla_{\boldsymbol{\theta}}f(\boldsymbol{\theta},H)$ is
$L$-Lipschitz for all $\boldsymbol{\theta}\in\Theta_{1}$;

(iii)

$H\mapsto f(\boldsymbol{\theta},H)$ is $\rho$-strongly convex for all
$\boldsymbol{\theta}\in\mathbb{R}^{p}$.

Then $g$ belongs to
$\operatorname{\textup{{Srg}}}_{L^{\prime\prime},\rho}(f_{*},\boldsymbol{\theta}^{\star})$
for some $L^{\prime\prime}>0$. When $f$ is jointly convex, then $f_{*}$ is
also convex and we can choose $L^{\prime\prime}=L$. $\blacktriangle$

######  Example 5.5 (Multi-convex Variational Surrogates).

Let
$f:\mathbb{R}^{p_{1}}\times\dots\times\mathbb{R}^{p_{m}}\times\mathbb{R}^{p_{m+1}}\rightarrow\mathbb{R}$
be a $(m+1)$-block multi-convex function and let
$\Theta_{1}\subseteq\mathbb{R}^{p}$ and
$\Theta_{i}\subseteq\mathbb{R}^{p_{i}}$ for $i=1,\dots,m+1$ be convex sets.
Denote $\boldsymbol{\Theta}:=\Theta_{1}\times\dots\times\Theta_{m}$. Define a
function
$f_{*}:\boldsymbol{\theta}\mapsto\inf_{H\in\Theta_{p_{m+1}}}f(\boldsymbol{\theta},H)$.
Then for each $\boldsymbol{\theta}^{\star}\in\boldsymbol{\Theta}$, the
following function

$\displaystyle g:\boldsymbol{\theta}\mapsto
f(\boldsymbol{\theta},H^{\star}),\quad
H^{\star}\in\operatorname*{arg\,min}_{H\in\boldsymbol{\Theta}_{p_{m+1}}}f(\boldsymbol{\theta}^{\star},H)$
(40)

is $m$-block multi-convex over $\boldsymbol{\Theta}$ and satisfies $g\geq f$
and $g(\boldsymbol{\theta}^{\star})=f(\boldsymbol{\theta}^{\star})$. Further
assume that

(i)

$\boldsymbol{\theta}\mapsto f(\boldsymbol{\theta},H)$ is differentiable for
all $H\in\Theta_{p_{m+1}}$ and
$\boldsymbol{\theta}\mapsto\nabla_{\boldsymbol{\theta}}f(\boldsymbol{\theta},H)$
is $L^{\prime}$-Lipschitz for all $H\in\Theta_{p_{m+1}}$ ;

(ii)

$H\mapsto\nabla_{\boldsymbol{\theta}}f(\boldsymbol{\theta},H)$ is
$L$-Lipschitz for all $\boldsymbol{\theta}\in\boldsymbol{\Theta}$;

(iii)

$H\mapsto f(\boldsymbol{\theta},H)$ is $\mu$-strongly convex for all
$\boldsymbol{\theta}\in\mathbb{R}^{p_{m+1}}$.

Let $\mathbb{J}=\\{J_{1},\dots,J_{m}\\}$, where $J_{i}$ denotes the $i$th
coordinate block corresponding to $\mathbb{R}^{p_{i}}$. Then $g$ belongs to
$\operatorname{\textup{{Srg}}}_{L^{\prime\prime},\mu}^{\mathbb{J}}(f_{*},\boldsymbol{\theta}^{\star})$
for some $L^{\prime\prime}>0$. When $f$ is jointly convex, then $f_{*}$ is
also convex and we can choose $L^{\prime\prime}=L$. $\blacktriangle$

######  Example 5.6 (Refining block structure).

Suppose $f:\boldsymbol{\Theta}\subseteq\mathbb{R}^{p}\rightarrow\mathbb{R}$
with $\boldsymbol{\Theta}$ convex, and we have a convex surrogate
$g\in\operatorname{\textup{{Srg}}}_{L,\rho}(f,\boldsymbol{\theta}^{\star})$.
Let $\mathbb{J}$ denote an arbitrary set of coordinate bloks for
$\mathbb{R}^{p}$ (e.g., $\mathbb{J}=\\{\\{1\\},\dots,\\{p\\}\\}$). Then one
can also view the fuction $g$ as a multi-convex surrogate in
$\operatorname{\textup{{Srg}}}_{L,\rho}^{\mathbb{J}}(f_{*},\boldsymbol{\theta}^{\star})$.
Instead of minimizing $g$ over the convex set $\boldsymbol{\Theta}$ with
respect to all $p$ coordinates, one can block-minimize $g$ with respect to the
block structure $\mathbb{J}$. For instance, if
$\mathbb{J}=\\{\\{1\\},\dots,\\{p\\}\\}$ and if $J_{i}(n)$ in Algorithm 2 is
chosen uniformly at random from $\mathbb{J}$, then Algorithm 2 becomes a
random coordinate descent for $m$ iterations, which choses a random coordinate
$i\in\\{1,\dots,p\\}$ and minimize $g$ over $\boldsymbol{\Theta}$ with respect
to the $i$th coordinate. (See, e.g., [49]).

A similar remark holds for multi-convex surrogates. Namely, suppose
$g\in\operatorname{\textup{{Srg}}}_{L,\rho}^{\mathbb{J}}(f,\boldsymbol{\theta}^{\star})$,
where $\mathbb{J}=\\{J_{1},\dots,J_{m}\\}$ denotes a set of coordinate blocks
for $\mathbb{R}^{p}$. Let $\mathbb{J}^{\prime}$ denote an arbitrary set of
coordinate bloks for $\mathbb{R}^{p}$ that refines $\mathbb{J}$. Namely, each
coordinate block $J\in\mathbb{J}^{\prime}$ is a subset of some $J_{i}$,
$i=1,\dots,m$. Then one can also view the fuction $g$ as a multi-convex
surrogate in
$\operatorname{\textup{{Srg}}}_{L^{\prime\prime},\rho}^{\mathbb{J}^{\prime}}(f_{*},\boldsymbol{\theta}^{\star})$
and accordingly, one can block-minimize $g$ using the finer block structure
given by $\mathbb{J}^{\prime}$ instead of the original block structure
$\mathbb{J}$. $\blacktriangle$

See [34] for other types of surrogate functions such as quadratic, Jensen, and
saddle point surrogates.

###  5.2. Applications

In this subsection, we discuss some applications of our general results in the
setting of online matrix and tensor factorization problems.

######  Application 5.7 (Online Matrix Factorization [37, 30]).

Consider the matrix factorization loss $f(W,H):=\lVert
X-WH\rVert^{2}+\lambda\lVert H\rVert_{1}$, where $X\in\mathbb{R}^{p\times d}$
is a given data matrix to be factorized into the product of dictionary
$W\in\mathbb{R}^{p\times r}$ and the code $H\in\mathbb{R}^{r\times d}$ with
$\lambda\geq 0$ being the $L_{1}$-regularization parameter for $H$. Clearly
$f$ is two-block multi-convex and differentiable with respect to $W$ with
gradient $\nabla_{W}f(W,H)=(X-WH)H^{T}$. Fix compact and convex constraint
sets $\Theta_{1}\subseteq\mathbb{R}^{p\times r}$ and
$\Theta_{2}\subseteq\mathbb{R}^{r\times d}$. Then $W\mapsto\nabla f(W,H)$ and
$H\mapsto\nabla f(W,H)$ are both Lipschitz over the compact and convex set
$\Theta_{1}\times\Theta_{2}$.

Suppose we have a sequence of data matrices $(X_{n})_{n\geq 1}$ and a sequence
of weights $(w_{n})_{n\geq 1}$. For each $n\geq 1$, the function
$f_{n}:W\mapsto\ell(X_{n},W)=\inf_{H\in\Theta_{2}}\lVert
X_{n}-WH\rVert^{2}+\lambda\lVert H\rVert_{1}$ is the minimum reconstruction
error for factorizing $X_{n}$ using the dictionary matrix $W$. The
corresponding empirical loss minimization problem is

$\displaystyle
W_{n}\in\operatorname*{arg\,min}_{W\in\Theta_{1}}\bar{f}_{n}(W):=(1-w_{n})\bar{f}_{n-1}(W)+w_{n}\ell(X_{n},W).$
(41)

If we have a target distribution $\pi$ for the data matrices $X_{n}$, then the
corresponding expected loss minimization problem is

$\displaystyle
W\in\operatorname*{arg\,min}_{W\in\Theta_{1}}\,f(W):=\mathbb{E}_{X\sim\pi}\left[\ell(X,W)\right].$
(42)

The latter is known as the online matrix factorization problem introduced in
[37].

In order to apply Algorithm 1 in this setting, denote
$H_{n}\in\operatorname*{arg\,min}_{H\in\Theta_{2}}\lVert
X_{n}-W_{n-1}H\rVert_{F}^{2}+\lambda\lVert H\rVert_{1}$, which is an optimal
code for factorizing $X_{n}$ using the previous dictionary $W_{n-1}$. Then the
function $g_{n}(W)=\lVert X_{n}-WH_{n}\rVert_{F}^{2}+\lambda\lVert
H_{n}\rVert_{1}$ is a surrogate of $f_{n}:W\mapsto\ell(X_{n},W)$ at $W_{n-1}$
and belongs to $\operatorname{\textup{{Srg}}}_{L^{\prime\prime},0}(f_{n})$ for
some $L^{\prime\prime}>0$. A simple calculation shows that the resulting
averaged surrogate function $\bar{g}_{n}$ becomes

$\displaystyle\bar{g}_{n}(W)=\textup{tr}(WA_{n}W^{T})-2\textup{tr}(WB_{n})+C_{n},$
(43)

where the matrices $A_{n},B_{n},C_{n}$ are recursively defined as

$\displaystyle A_{n}=(1-w_{n})A_{n-1}+w_{n}H_{n}H_{n}^{T};\quad
B_{n}=(1-w_{n})B_{n-1}+w_{n}H_{n}X_{n}^{T};\quad
C_{n}=(1-w_{n})\bar{X}_{n-1}+w_{n}X_{n}X_{n}^{T}.$ (44)

Thus Algorithm 1 in this case reduces to

$\displaystyle\textit{Upon arrival of
$X_{n}$:\qquad}\begin{cases}H_{n}=\operatorname*{arg\,min}_{H\in\Theta_{2}\subseteq\mathbb{R}^{r\times
d}}\lVert X_{n}-W_{n-1}H\rVert_{F}^{2}+\lambda\lVert H\rVert_{1}\\\
A_{n}=(1-w_{n})A_{n-1}+w_{n}H_{n}H_{n}^{T}\\\
B_{n}=(1-w_{n})B_{n-1}+w_{n}H_{n}X_{n}^{T}\\\
W_{n}=\operatorname*{arg\,min}_{W\in\Theta_{1}\subseteq\mathbb{R}^{p\times
r}}\left(\textup{tr}(WA_{n}W^{T})-2\textup{tr}(WB_{n})\right),\end{cases}$
(45)

which is the online matrix factorization algorithm (OMF) proposed in [37] for
i.i.d. data matrices. Later in [30], this algorithm was analyzed in a
Markovian data setting. Theorems 4.2 and 4.3 establish rate of convergence of
the online matrix factorization algorithm (45) under Markovian data setting.
To the author’s best knowledge, such a result was not known for online matrix
factorization even under the i.i.d. data setting. $\blacktriangle$

######  Application 5.8 (Subsampled Online Matrix Factorization [38]).

While the OMF algorithm (45) is efficient in handling matrices with a large
number of columns $d$, its computational cost with respect to the data
dimension $p$ (i.e., the number of rows) is not reduced. In [38], subsampled
OMF was proposed to improve the computational efficiency with respect to $p$,
by using only a random subsample of rows. For a preliminary version, consider
the following variant of the online matrix factorization algorithm 45, which
uses a random coordinate descent on subsampled rows of $W_{n}$:

$\displaystyle\textit{Upon arrival of
$X_{n}$:\qquad}\begin{cases}H_{n}=\operatorname*{arg\,min}_{H\in\Theta_{2}\subseteq\mathbb{R}^{r\times
d}}\lVert X_{n}-W_{n-1}H\rVert_{F}^{2}+\lambda\lVert H\rVert_{1}\\\
A_{n}=(1-w_{n})A_{n-1}+w_{n}H_{n}H_{n}^{T}\\\
B_{n}=(1-w_{n})B_{n-1}+w_{n}H_{n}X_{n}^{T}\\\ J\leftarrow\textup{Random subset
of $\\{1,\dots,p\\}$}\\\
W_{n}=\operatorname*{arg\,min}_{W\in\Theta_{1}\subseteq\mathbb{R}^{p\times
r}}\left(\textup{tr}(WA_{n}W^{T})-2\textup{tr}(WB_{n})\right)\\\ \hskip
28.45274pt\textup{while freezing rows of $W_{n}$ not in $J$ to the rows of
$W_{n-1}$}.\end{cases}$ (46)

For instance, $J$ may be taken as the uniform subset of $\\{1,\dots,p\\}$ of a
fixed size $p^{\prime}<p$ or it may contain each index $i\in\\{1,\dots,p\\}$
independently with a fixed probability. All our main results in this paper
(Theorems 4.1, 4.2, and 4.3) apply to the OMF algorithm in (46).

However, note that the code computation for $H_{n}$ in (46) still involves
solving a least squares problem with $p$ dimensional data matrix $X_{n}$. In
order to fully reduce the dependency on $p$, one may also use row subsampling
to compute only $p^{\prime}$ rows of $H_{n}$. The approach used in [38] was to
replace $X_{n}$ with an averaged $p^{\prime}$-dimenisonal matrix
$\bar{X}_{n}^{(i)}$ assuming $X_{n}=X^{(i)}$, where a finite pool of data
matrices $\mathfrak{X}=\\{X^{(1)},\dots,X^{(l)}\\}$ was assumed and
$\bar{X}^{(i)}_{n}$ is a recursively defined matrix of size $p^{\prime}\times
d$ based on previous occurances of the matrix $X^{(i)}$. See [38, Alg. 3] for
more details.

An almost sure convergence to stationary points of subsampled OMF algorithm
under i.i.d. data samples was shown in [38]. The analysis there is based on a
general convergence result on stochastic approxiamate majorization-
minimization algorithm (7) using strongly convex $\varepsilon$-approximate
surrogate functions (see [38, Prop. 3]). A similar convergence result is
retained by Theorem 4.1, although the assumptions slightly different. In
addition, Theorems 4.2-4.3 also provide rate of convergence. $\blacktriangle$

######  Application 5.9 (Online CP-dictionary Learning [32]).

Suppose we have $b$ observed $n$-mode tensor-valued signals
$\mathbf{X}_{1},\dots,\mathbf{X}_{N}\in\mathbb{R}_{\geq
0}^{I_{1}\times\dots\times I_{m}}$. Consider the following tensor-valued
dictionary learning problem

$\displaystyle[\mathbf{X}_{1},\dots,\mathbf{X}_{b}]\approx[\mathbf{D}_{1},\dots,\mathbf{D}_{r}]\times_{m+1}H,$
(47)

where $\times_{m+1}$ denotes the mode-$(m+1)$ tensor-matrix product and we
impose the tensor dictionary atoms
$\mathbf{D}_{1},\dots,\mathbf{D}_{r}\in\mathbb{R}^{I_{1}\times\dots\times
I_{m}\times r}$ to be of rank-1 and $H\in\mathbb{R}^{R\times b}$ is called a
code matrix. Equivalently, we assume that there exist loading matrices
$[U^{(1)},\dots,U^{(n)}]\in\mathbb{R}_{\geq 0}^{I_{1}\times
r}\times\dots\times\mathbb{R}_{\geq 0}^{I_{n}\times r}$ such that

$\displaystyle[\mathbf{D}_{1},\dots,\mathbf{D}_{r}]$
$\displaystyle=\operatorname{\texttt{Out}}(U^{(1)},\dots,U^{(m)})$ (48)
$\displaystyle:=\left[\bigotimes_{k=1}^{m}U^{(k)}(:,1),\,\bigotimes_{k=1}^{m}U^{(k)}(:,2),\,\dots\,,\bigotimes_{k=1}^{m}U^{(k)}(:,r)\right]\in\mathbb{R}_{\geq
0}^{I_{1}\times\dots\times I_{m}\times r},$ (49)

where $U^{(k)}(:,j)$ denotes the $i^{\textup{th}}$ column of the $I_{k}\times
r$ matrix $U^{(k)}$ and $\otimes$ denotes the outer product. Since we impose a
CANDECOMP/PARAFAC (CP) [48, 19, 8] structure for the tensor-valued dictionary
$[\mathbf{D}_{1},\dots,\mathbf{D}_{r}]$, the above is called a CP-dictionary
learning problem introduced in [32]. When $m=1$, it reduces to the standard
vector-valued dictionary learning problem.

In [32], the following online CP-dictionary learning problem was proposed and
analyzed. Fix compact and convex constraint sets for code and loading matrices
$\mathcal{C}^{\textup{code}}\subseteq\mathbb{R}^{R\times b}$ and
$\mathcal{C}^{(i)}\subseteq\mathbb{R}^{I_{i}\times r}$, $i=1,\dots,n$,
respectively. Write
$\boldsymbol{\Theta}:=\mathcal{C}^{(1)}\times\cdots\times\mathcal{C}^{(n)}$.
For each $\mathcal{X}\in\mathbb{R}_{\geq 0}^{I_{1}\times\dots\times
I_{n}\times b}$,
$\mathbf{U}:=[U^{(1)},\dots,U^{(n)}]\in\mathbb{R}^{I_{1}\times
r}\times\dots\times\mathbb{R}^{I_{n}\times r}$, $H\in\mathbb{R}^{R\times b}$,
define

$\displaystyle\ell(\mathcal{X},\mathbf{U},H)$
$\displaystyle:=\lVert\mathcal{X}-\operatorname{\texttt{Out}}(\mathbf{U})\times_{n+1}H\rVert_{F}^{2}+\lambda\lVert
H\rVert_{1},$ (50) $\displaystyle\ell(\mathcal{X},\mathbf{U})$
$\displaystyle:=\inf_{H\in\mathcal{C}^{\textup{code}}}\,\,\ell(\mathcal{X},\mathbf{U},H),$
(51)

where $\lambda\geq 0$ is a regularization parameter. Fix a sequence of non-
increasing weights $(w_{t})_{t\geq 0}$ in $(0,1]$. Here $\mathcal{X}$ denotes
a minibatch of $b$ tensors in $\mathbb{R}^{I_{1}\times\dots\times I_{n}}$, so
minimizing $\ell(\mathcal{X},\mathbf{U})$ with respect to $\mathbf{U}$ amounts
to fitting the CP-dictionary $\operatorname{\texttt{Out}}(\mathbf{U})$ to the
minibatch of $b$ tensors in $\mathcal{X}$. The corresponding empirical loss
minimization problem is

$\displaystyle\mathbf{U}_{n}\in\operatorname*{arg\,min}_{\mathbf{U}\in\boldsymbol{\Theta}}\bar{f}_{n}(\mathbf{U}):=(1-w_{n})\bar{f}_{n-1}(\mathbf{U})+w_{n}\ell(\mathcal{X}_{n},\mathbf{U}).$
(52)

If we have a target distribution $\pi$ for the data tensors $\mathcal{X}_{n}$,
then the corresponding expected loss minimization problem is

$\displaystyle
W\in\operatorname*{arg\,min}_{W\in\Theta_{1}}\,f(W):=\mathbb{E}_{X\sim\pi}\left[\ell(\mathcal{X},W)\right].$
(53)

In order to solve (52) and (53), the following online CP-dictionary learning
algorithm was poposed and analyzed in [32]:

$\displaystyle\textit{Upon arrival of
$\mathcal{X}_{n}$:\quad}\begin{cases}H_{n}=\operatorname*{arg\,min}_{H\in\mathcal{C}^{\textup{code}}\subseteq\mathbb{R}^{r\times
b}}\ell(\mathcal{X}_{n},\mathbf{U}_{n-1},H)\\\
\bar{g}_{n}(\mathbf{U})=(1-w_{n})\bar{g}_{n-1}(\mathbf{U})+w_{n}\ell(\mathcal{X}_{n},\mathbf{U},H_{n})\\\
\textup{for $i=1,\dots,m$}:\\\ \quad
U_{n}^{(i)}\in\operatorname*{arg\,min}_{\begin{subarray}{c}U\in\mathcal{C}^{(i)}\\\
\lVert U-U_{n-1}^{(i)}\rVert\leq
c^{\prime}w_{n}\end{subarray}}\,\,\bar{g}_{n}(U_{n}^{(1)},\dots,U_{n}^{(i-1)},U,U_{n-1}^{(i+1)},\dots,U_{n-1}^{(m)}).\end{cases}$
(54)

An almost sure convergence to the stationary points of the above algorithm
under Markovain data setting was obtained in [32]. In fact, the algorithm (54)
is a special case of Algorithm 1 with $m$-block multi-convex variational
surrogate corresponding (see Example 5.5) to the function $\ell$ in (50) with
$\mathbb{J}=\\{J_{1},\dots,J_{m}\\}$, where $J_{i}$ denotes the coordinate
block corresponding to the $i$th loading matrix in $\mathbb{R}^{I_{i}\times
r}$ with cyclic block coordinate descent in Algorithm 2.

Our generalized algorithm and refined analysis improve the results in [32] in
multiple ways. First, Theorems 4.2 and 4.3 provide rate of convergence results
for the online CP-dictionary learning algorithm (54), which has not been known
before. Second, due to the flexibility of using approximate surrogate
functions, the convergence results of (54) hold under inexact code or factor
matrix computation, following a similar argument as in [38]. Third, one can
only optimize a small number of subsampled rows when updating the factor
matrices $U_{n}^{(i)}$ in (54) by using a refined block coordinate structure
(see Example 5.6). A similar idea of using subsampling for tensor CP-
decomposition was also used recently in [25].

Lastly, we remark that random row subsampling can be utilized for the code
computation for $H_{n}$ in (54), which essentially involves solving a least-
squares problem with the mode-$(m+1)$ unfolding of $\mathcal{X}_{n}$, which
has size $(I_{1}\times\cdots\times I_{m})\times b$. Since the number of rows
of such matrix $I_{1}\times\cdots\times I_{m}$ can be very large, the
computational gain in this approach should be significant. We believe it would
be straightforward to adapt the approach in [38] for subsampled online matrix
factorization for a subsampled online CP-dictionary learning setting. For our
theoretical results to apply, one only needs to verify that the resulting
inexact code computation for $H_{n}$ yields $\varepsilon_{n}$-approximate
surrogates $\bar{g}_{n}$ that verifies the assumption (A4). We do not proceed
with this line of research in the present paper. $\blacktriangle$

##  6\. Preliminary lemmas

In this section, we give some preliminary lemmas that we will use in our
analysis in the following sections.

######  Lemma 6.1.

Let $(\boldsymbol{\theta}_{n})_{n\geq 1}$ be an output of Algorithm 1. Then
for all $n\geq 1$, the following hold:

(i)

(Forward Monotonicity)
$\bar{g}_{n}(\boldsymbol{\theta}_{n-1})\geq\bar{g}_{n}(\boldsymbol{\theta}_{n})$;

(ii)

(Stability of Estimates I)
$\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert\leq
c^{\prime}w_{n}$;

(iii)

(Stability of Estimates II) If $c^{\prime}=\infty$ and the surrogates $g_{n}$
are $\rho$-strongly convex for some $\rho>0$, then under (A1),
$\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\rVert\leq
2Rw_{n}/\rho$;

###### Proof.

Consider the computation of $\boldsymbol{\theta}_{n+1}$ in Algorithm 2. Let
$J_{1},\dots,J_{m}\in\mathbb{J}$ denote the coordinate blocks used in order,
and let
$\boldsymbol{\theta}_{n}^{(1)},\dots,\boldsymbol{\theta}_{n}^{(m)}(=\boldsymbol{\theta}_{n+1})$
denote the outputs of (11) after each block minimization. Then
$\boldsymbol{\theta}_{n}^{(m)}=\boldsymbol{\theta}_{n}$. Denote
$\boldsymbol{\theta}_{n}^{(0)}=\boldsymbol{\theta}_{n-1}$. Note that
$\boldsymbol{\theta}_{n}^{(i)}$ is a minimizer of the convex function
$\bar{g}_{n}$ over the convex set $\boldsymbol{\Theta}_{n}^{J_{i}}$ defined in
(10). Also note that
$\boldsymbol{\theta}_{n}^{(i-1)}\in\boldsymbol{\Theta}_{n}^{J_{i}}$. Hence
$\bar{g}_{n}(\boldsymbol{\theta}_{n}^{(i)})-\bar{g}_{n}(\boldsymbol{\theta}_{n}^{(i-1)})\geq
0$ for each $i=1,\dots,m$. From this, we can easily conclude (i), since

$\displaystyle\bar{g}_{n}(\boldsymbol{\theta}_{n-1})-\bar{g}_{n}(\boldsymbol{\theta}_{n})=\sum_{i=1}^{n}\bar{g}_{n}(\boldsymbol{\theta}_{n}^{(i)})-\bar{g}_{n}(\boldsymbol{\theta}_{n}^{(i-1)})\geq
0.$ (55)

Next, (ii) is trivial by the search radius restriction in Algorithm 1. When
$c^{\prime}=\infty$, the statement is vacuous.

Now we show (iii). Our argument is inspired by the proof Mairal et al. [37].
By the second-order growth property (Lemma B.2),

$\displaystyle\frac{\rho}{2}\lVert\boldsymbol{\theta}_{n}^{(i)}-\boldsymbol{\theta}_{n}^{(i-1)}\rVert^{2}$
$\displaystyle\leq\bar{g}_{n}(\boldsymbol{\theta}_{n}^{(i-1)})-\bar{g}_{n}(\boldsymbol{\theta}_{n}^{(i)})$
(56)
$\displaystyle=w_{n}\left(g_{n+1}(\boldsymbol{\theta}_{n}^{(i-1)})-g_{n}(\boldsymbol{\theta}_{n}^{(i)})\right)+(1-w_{n})\left(\bar{g}_{n-1}(\boldsymbol{\theta}_{n}^{(i-1)})-\bar{g}_{n-1}(\boldsymbol{\theta}_{n}^{(i)})\right)$
(57) $\displaystyle\leq
w_{n}\left(g_{n}(\boldsymbol{\theta}_{n}^{(i-1)})-g_{n}(\boldsymbol{\theta}_{n}^{(i)})\right).$
(58)

Note that $g_{n}\geq f_{n}$ and
$g_{n}(\boldsymbol{\theta}_{n-1})=f_{n}(\boldsymbol{\theta}_{n-1})$ since
$g_{n}\in\operatorname{\textup{{Srg}}}_{L,\rho}^{\mathbb{J}}(f_{n},\boldsymbol{\theta}_{n-1})$.
Then since $f_{n}$ is $R$-Lipschitz continuous by Assumption (A1),

$\displaystyle
w_{n}\left(g_{n}(\boldsymbol{\theta}_{n-1})-g_{n}(\boldsymbol{\theta}_{n})\right)\leq
w_{n}\left(f_{n}(\boldsymbol{\theta}_{n-1})-f_{n}(\boldsymbol{\theta}_{n})\right)\leq
Rw_{n}\lVert\boldsymbol{\theta}_{n-1}-\boldsymbol{\theta}_{n}\rVert.$ (59)

Hence it follows that

$\displaystyle\frac{\rho}{2}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n+1}\rVert^{2}\leq
Rw_{n+1}\lVert\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n+1}\rVert,$ (60)

so the assertion follows. ∎

Next, we define a sequence $(\bar{\varepsilon}_{n})_{n\geq 1}$ recursively as

$\displaystyle\bar{\varepsilon}_{n}=(1-w_{n})\bar{\varepsilon}_{n-1}+w_{n}\varepsilon_{n},\qquad\bar{\varepsilon}_{0}=0,$
(61)

where $(\varepsilon_{n})_{n\geq 1}$ is the sequence of surrogate error
tolerance in Algorithm 1. Then

$\displaystyle\bar{\varepsilon}_{n}=\sum_{k=1}^{n}w^{n}_{k}\varepsilon_{k}\leq
w_{n}\sum_{k=1}^{n}\varepsilon_{k},$ (62)

where $w^{n}_{k}$ is defined in (2) and the inequality above uses
$w^{n}_{1}\leq w^{n}_{2}\leq\dots\leq w^{n}_{n}=w_{n}$ under (A4). Since
$g_{n}-f_{n}\geq-\varepsilon_{n}$ for all $n\geq 1$, it follows that

$\displaystyle\bar{g}_{n}-\bar{f}_{n}\geq-\bar{\varepsilon}_{n}\geq-
w_{n}\sum_{k=1}^{n}\varepsilon_{k}.$ (63)

######  Proposition 6.2.

Let $(\boldsymbol{\theta}_{n})_{n\geq 1}$ be an output of Algorithm 1. Assume
(A4). Then for each $n\geq 0$, the following hold:

(i)

$\bar{g}_{n+1}(\boldsymbol{\theta}_{n+1})-\bar{g}_{n}(\boldsymbol{\theta}_{n})\leq
w_{n+1}\left(\ell(\mathbf{x}_{n+1},\boldsymbol{\theta}_{n})-\bar{f}_{n}(\boldsymbol{\theta}_{n})\right)+w_{n}^{2}\left(\sum_{k=1}^{n+1}\varepsilon_{k}\right)$.

(ii)

$0\leq
w_{n+1}\left(\bar{g}_{n}(\boldsymbol{\theta}_{n})-\bar{f}_{n}(\boldsymbol{\theta}_{n})\right)\leq
w_{n+1}\left(\ell(\mathbf{x}_{n+1},\boldsymbol{\theta}_{n})-\bar{f}_{n}(\boldsymbol{\theta}_{n})\right)+\bar{g}_{n}(\boldsymbol{\theta}_{n})-\bar{g}_{n+1}(\boldsymbol{\theta}_{n+1})+w_{n}^{2}\left(\sum_{k=1}^{n+1}\varepsilon_{k}\right)$.

###### Proof.

We begin by observing that

$\displaystyle\bar{g}_{n+1}(\boldsymbol{\theta}_{n})$
$\displaystyle=(1-w_{n+1})\bar{g}_{n}(\boldsymbol{\theta}_{n})+w_{n+1}\ell(\mathbf{x}_{n+1},\boldsymbol{\theta}_{n})+w_{n+1}(g_{n+1}(\boldsymbol{\theta}_{n})-\ell(\mathbf{x}_{n+1},\boldsymbol{\theta}_{n}))$
(64)

for all $t\geq 0$. Hence

$\displaystyle\bar{g}_{n+1}(\boldsymbol{\theta}_{n+1})-\bar{g}_{n}(\boldsymbol{\theta}_{n})$
(65)
$\displaystyle\quad=\bar{g}_{n+1}(\boldsymbol{\theta}_{n+1})-\bar{g}_{n+1}(\boldsymbol{\theta}_{n})+\bar{g}_{n+1}(\boldsymbol{\theta}_{n})-\bar{g}_{n}(\boldsymbol{\theta}_{n})$
(66)
$\displaystyle\quad=\bar{g}_{n+1}(\boldsymbol{\theta}_{n+1})-\bar{g}_{n+1}(\boldsymbol{\theta}_{n})+(1-w_{n+1})\bar{g}_{n}(\boldsymbol{\theta}_{n})+w_{n+1}\ell(\mathbf{x}_{n+1},\boldsymbol{\theta}_{n})-\bar{g}_{n}(\boldsymbol{\theta}_{n})$
(67)
$\displaystyle\quad=\bar{g}_{n+1}(\boldsymbol{\theta}_{n+1})-\bar{g}_{n+1}(\boldsymbol{\theta}_{n})+w_{n+1}(\ell(\mathbf{x}_{n+1},\boldsymbol{\theta}_{n})-\bar{f}_{n}(\boldsymbol{\theta}_{n}))+w_{n+1}(\bar{f}_{n}(\boldsymbol{\theta}_{n})-\bar{g}_{n}(\boldsymbol{\theta}_{n}))$
(68) $\displaystyle\hskip
142.26378pt+w_{n+1}(g_{n+1}(\boldsymbol{\theta}_{n})-\ell(\mathbf{x}_{n+1},\boldsymbol{\theta}_{n})).$
(69)

Now note that
$\bar{g}_{n+1}(\boldsymbol{\theta}_{n+1})-\bar{g}_{n+1}(\boldsymbol{\theta}_{n})\leq
0$ by Lemma 6.1 (i). Also, recall that $\bar{g}_{n}-\bar{f}_{n}\geq-
w_{n}\sum_{k=1}^{n}\varepsilon_{k}$ from (63) and $0\leq
g_{n+1}(\boldsymbol{\theta}_{n})-\ell(\mathbf{x}_{n+1},\boldsymbol{\theta}_{n})\leq\varepsilon_{n+1}$
since $g_{n+1}$ is a multi-convex $\varepsilon_{n}$-approximate surrogate of
$f_{n+1}=\ell(\mathbf{x}_{n+1},\cdot)$ at $\boldsymbol{\theta}_{n}$. Thus the
inequalities in both (i) and (ii) follow using $w_{n+1}\leq w_{n}$. ∎

######  Lemma 6.3.

Let $(a_{n})_{n\geq 0}$ and $(b_{n})_{n\geq 0}$ be sequences of nonnegative
real numbers such that $\sum_{n=0}^{\infty}a_{n}b_{n}<\infty$. Then the
following hold.

(i)

$\displaystyle\min_{1\leq k\leq
n}b_{k}\leq\frac{\sum_{k=0}^{\infty}a_{k}b_{k}}{\sum_{k=1}^{n}a_{k}}=O\left(\left(\sum_{k=1}^{n}a_{k}\right)^{-1}\right)$.

(ii)

Further assume $\sum_{n=0}^{\infty}a_{n}=\infty$ and
$|b_{n+1}-b_{n}|=O(a_{n})$. Then $\lim_{n\rightarrow\infty}b_{n}=0$.

###### Proof.

(i) follows from noting that

$\displaystyle\left(\sum_{k=1}^{n}a_{k}\right)\min_{1\leq k\leq
n}b_{k}\leq\sum_{k=1}^{n}a_{k}b_{k}\leq\sum_{k=1}^{\infty}a_{k}b_{k}<\infty.$
(70)

The proof of (ii) is omitted and can be found in [35, Lem. A.5]. ∎

##  7\. Key lemmas and proofs of main results

In this section, we state all key lemmas without proofs and derive the main
results (Theorems 4.1, 4.2, and 4.3) assuming them. The key lemmas stated in
this section will be proved in the subsequent sections.

###  7.1. Key Lemmas

In this subsection, we state all key lemmas that are sufficient to derive the
main results in this paper.

First, Lemma 7.1 states some general concentration inequalities of recursively
defined functions similarly as the empirical and the surrogate loss functions
in (1) and (8). One may regard them as the classical Glivenko-Cantelli theorem
for a general weighting scheme.

######  Lemma 7.1.

Suppose $w_{n}\in(0,1]$, $n\geq 1$ are non-increasing and satisfy
$w_{n}^{-1}-w_{n-1}^{-1}\leq 1$ for all sufficiently large $n\geq 1$. Fix a
continuous and bounded function
$\psi:\mathfrak{X}\times\boldsymbol{\Theta}\rightarrow\mathbb{R}^{d}$. Define
functions
$\bar{\psi}(\cdot):=\mathbb{E}_{\mathbf{y}\sim\pi}\left[\psi(\varphi(\mathbf{y}),\cdot)\right]$
and $\bar{\psi}_{n}:\boldsymbol{\Theta}\rightarrow\mathbb{R}^{d}$ recursively
as $\bar{\psi}_{0}\equiv\mathbf{0}$ and

$\displaystyle\bar{\psi}_{n}(\cdot)=(1-w_{n})\bar{\psi}_{n-1}(\cdot)+w_{n}\psi(\mathbf{x}_{n},\cdot).$
(71)

Then under Assumptions (A1)-(A3), there exists a constant $C>0$ such that for
all $n\geq 1$,

$\displaystyle\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\lVert\bar{\psi}(\boldsymbol{\theta})-\mathbb{E}[\bar{\psi}_{n}(\boldsymbol{\theta})]\right\rVert\leq
Cw_{n},\quad\mathbb{E}\left[\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\lVert\bar{\psi}(\boldsymbol{\theta})-\bar{\psi}_{n}(\boldsymbol{\theta})\right\rVert\right]\leq
Cw_{n}\sqrt{n}.$ (72)

Furthermore, if $w_{n}\sqrt{n}=O(1/(\log n)^{1+\varepsilon})$ for some
$\varepsilon>0$, then
$\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\lVert\bar{\psi}(\boldsymbol{\theta})-\bar{\psi}_{n}(\boldsymbol{\theta})\right\rVert\rightarrow
0$ as $t\rightarrow\infty$ almost surely.

We remark that a similar result appeard in [35, Lem B.7], which states the
second inequality in (72) as well as the last almost sure convergence
statement under the stronger condition of $\sum_{n=1}^{\infty}w_{n}=\infty$
and $\sum_{n=1}^{\infty}w_{n}^{2}\sqrt{n}<\infty$. We will refer to the latter
condition as the ‘strong square-summability’ condition. These conditions were
necessary in order to use Lemma 6.3 (ii) to deduce the almost sure
convergence. We give a directly argument using Borel-Cantelli lemma without
these additional condition, for which the weaker condition of
$w_{n}\sqrt{n}=O(1/(\log n)^{1+\varepsilon})$ for some $\varepsilon>0$ is
sufficient.

Next, Lemma 7.2 states a series of finite variation statements that provide a
basis for the forthcoming arguments. Most of the statements also appeared in
the literature (see, e.g., [37, 35, 30, 32]) in some special cases, where the
strong square summability condition was used. We give an improved argument
only using the square summability condition except for the last item.

######  Lemma 7.2.

Let $(\boldsymbol{\theta}_{n})_{n\geq 1}$ be an output of Algorithm 1. Suppose
(A1)-(A4). Then the following hold:

(i)

$\displaystyle\sum_{n=1}^{\infty}\mathbb{E}\left[w_{n+1}\left(\ell(\mathbf{x}_{n+1},\boldsymbol{\theta}_{n})-\bar{f}_{n}(\boldsymbol{\theta}_{n})\right)\right]^{+}<\infty$;

(ii)

$\displaystyle\sum_{n=0}^{\infty}\left(\mathbb{E}\left[\bar{g}_{n+1}(\boldsymbol{\theta}_{n+1})-\bar{g}_{n}(\boldsymbol{\theta}_{n})\right]\right)^{+}<\infty$;

(iii)

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

(iv)

$\displaystyle\sum_{n=0}^{\infty}w_{n+1}\left|\mathbb{E}\left[f(\boldsymbol{\theta}_{n})-\bar{f}_{n}(\boldsymbol{\theta}_{n})\right]\right|<\infty$;

(v)

Suppose the optional condition in (A4) holds. Then
$\displaystyle\mathbb{E}\left[\sum_{n=0}^{\infty}w_{n+1}\left|f(\boldsymbol{\theta}_{n})-\bar{f}_{n}(\boldsymbol{\theta}_{n})\right|\right]<\infty$.

The following lemma is one of the key innovations in the present work that
allows us to obtain convergence rate bounds stated in Theorems 4.2 and 4.3.
Roughly speaking, it gives gradient versions of the finite variation
statements in Lemma 7.2.

######  Lemma 7.3 (Variation of gradients).

Let $(\boldsymbol{\theta}_{n})_{n\geq 1}$ be the output of Algorithm 1.
Suppose (A1)-(A4). Then there exist a constant $M>0$ independent of
initialization $\boldsymbol{\theta}_{0}$ such that the following hold:

(i)

We have

$\displaystyle\mathbb{E}\left[\sum_{n=1}^{\infty}w_{n}\,\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})-\nabla\bar{f}_{n}(\boldsymbol{\theta}_{n})\rVert^{2}+\left\lVert\mathbb{E}\left[\nabla\bar{f}_{n}(\boldsymbol{\theta}_{n})\right]-\mathbb{E}\left[\nabla
f(\boldsymbol{\theta}_{n})\right]\right\rVert^{2}\right]<\infty.$ (73)

(ii)

If in addition the optional condition in (A4) holds, then

$\displaystyle\mathbb{E}\left[\sum_{n=1}^{\infty}w_{n+1}\,\left(\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})-\nabla\bar{f}_{n}(\boldsymbol{\theta}_{n})\rVert^{2}+\lVert\nabla\bar{f}_{n}(\boldsymbol{\theta}_{n})-\nabla
f(\boldsymbol{\theta}_{n})\rVert^{2}+\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})-\nabla
f(\boldsymbol{\theta}_{n})\rVert^{2}\right)\right]<\infty.$ (74)

Based on Lemmas 7.1, 7.2, and 7.3, one can deduce Theorem 4.2 as well as
Theorem 4.1 for case C1. A nice property of Algorithm 1 with convex surrogates
without radius restriction ($c^{\prime}=\infty$) is that the iterates
$\boldsymbol{\theta}_{n}$ are exact minimizers of the convex averaged
surrogates $\bar{g}_{n}$ over $\boldsymbol{\Theta}$. Hence
$-\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})$ lies in the normal cone of
$\boldsymbol{\Theta}$ at $\boldsymbol{\theta}_{n}$ for each $n\geq 1$. Thus,
in order for asymptotic stationarity with respect to the empricial loss
$\bar{f}_{n}$ or the expected loss $f$, one only needs to show that every
convergent subsequence of
$\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})-\nabla
h_{n}(\boldsymbol{\theta}_{n})\rVert$ vanish almost surely, where
$h_{n}\in\\{\bar{f}_{n},f\\}$ for $n\geq 1$. However, when $\bar{g}_{n}$ is
only approximately minimized in each iteration $n$,
$-\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})$ may be close to but not within
the normal cone of $\boldsymbol{\Theta}$ at $\boldsymbol{\theta}_{n}$.
Moreover, although the technique of radius restriction in Algorithm 1 gives
stability of estimates (see Lemma 6.1 (ii)) in this inexact case, it
introduces additional radius constraints so that
$-\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})$ may be normal the trust region
boundary, which is the sphere of distance $c^{\prime}w_{n}$ centered at
$\boldsymbol{\theta}_{n-1}$. In order to handle these issues, in Lemma 7.4, we
will show that despite $\boldsymbol{\theta}_{n}$ being an inexact minimizer of
$\bar{g}_{n}$ over $\boldsymbol{\Theta}$ intersected with a trust region, the
sequence $\boldsymbol{\theta}_{n}$, $n\geq 1$ still verifies stationarity for
$\bar{g}_{n}$ in an asymptotic sense. For this, we need to take ‘convergent
subsequences‘ of the averaged surrogates $\bar{g}_{n}$, which can be easily
done under the compact parameterization assumption in (A6).

######  Lemma 7.4 (Asymptotic Surrogate Stationarity).

Assume (A1)-(A6) and let $(\boldsymbol{\theta}_{n})_{n\geq 1}$ be an output of
Algorithm 1 with $c^{\prime}<\infty$. Let $(t_{k})_{k\geq 1}$ be any sequence
such that $\boldsymbol{\theta}_{t_{k}}$ and $\bar{g}_{t_{k}}$ converges almost
surely. Then
$\boldsymbol{\theta}_{\infty}=\lim_{k\rightarrow\infty}\boldsymbol{\theta}_{t_{k}}$
is almost surely a stationary point of
$\bar{g}_{\infty}:=\lim_{k\rightarrow\infty}\bar{g}_{t_{k}}$ over
$\boldsymbol{\Theta}$.

Finally, Lemma 7.5 is the last ingredient for proving Theorem 4.3, which
states that the optimality for the averaged surrogates $\bar{g}_{n}$ is
summable with weights $w_{n}$.

######  Lemma 7.5 (Finite variation of surrogate optimality).

Assume (A1)-(A5). Then there exist a constant $M>0$ independent of
initialization $\boldsymbol{\theta}_{0}$ such that

$\displaystyle\mathbb{E}\left[\sum_{n=1}^{\infty}w_{n+1}\left|-\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\langle\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n}),\,\frac{\boldsymbol{\theta}-\boldsymbol{\theta}_{n}}{\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert}\right\rangle\right|\right]<\infty.$
(75)

###  7.2. Proof of Theorem 4.2

In this subsection, we prove Theorem 4.2 assuming Lemmas 7.1, 7.2, and 7.3.

###### Proof of Theorem 4.2.

Suppose (A1)-(A5) hold. Throughout the proof, we let $M>0$ denote a constant
independent of the initial estimate $\boldsymbol{\theta}_{0}$. By Lemma 7.2
(iii) and Lemma 7.3 (i), we get

$\displaystyle\mathbb{E}\left[\sum_{n=1}^{\infty}w_{n+1}\left(\left|\bar{g}_{n}(\boldsymbol{\theta}_{n})-\bar{f}_{n}(\boldsymbol{\theta}_{n})\right|+\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})-\nabla\bar{f}_{n}(\boldsymbol{\theta}_{n})\rVert^{2}\right)\right]<\infty.$
(76)

Recall that $\mathbb{E}[|X|]<\infty$ implies $|X|<\infty$ almost surely for
any random variable $X$. It follows that the summation above is almost surely
finite. Then (19) follows from Lemma 6.3.

Next by Lemma 7.2 (iii)-(iv) and Jensen’s inequality, we get

$\displaystyle\mathbb{E}\left[\sum_{n=1}^{\infty}w_{n}\left|\mathbb{E}[\bar{g}_{n}(\boldsymbol{\theta}_{n})]-\mathbb{E}[f(\boldsymbol{\theta}_{n})]\right|\right]$
(77)
$\displaystyle\qquad\leq\mathbb{E}\left[\sum_{n=1}^{\infty}w_{n}\left(\mathbb{E}\left[\left|\bar{g}_{n}(\boldsymbol{\theta}_{n})-\bar{f}_{n}(\boldsymbol{\theta}_{n})\right|\right]+\left|\mathbb{E}[\bar{f}_{n}(\boldsymbol{\theta}_{n})]-\mathbb{E}[f(\boldsymbol{\theta}_{n})]\right|\right)\right]<\infty.$
(78)

Similarly, using Lemma 7.3 (i), we can deduce

$\displaystyle\sum_{n=1}^{\infty}w_{n}\,\left\lVert\mathbb{E}\left[\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})\right]-\mathbb{E}\left[\nabla
f(\boldsymbol{\theta}_{n})\right]\right\rVert^{2}<\infty.$ (79)

Then the assertion follows from combining the above inequalities and using
Lemma 6.3.

Recall that by Lemma 7.1, there exists a constant $C>0$ such that for all
$n\geq 1$,

$\displaystyle\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\lvert
f(\boldsymbol{\theta})-\mathbb{E}[\bar{f}_{n}(\boldsymbol{\theta})]\right\rvert\leq
Cw_{n},\qquad\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\lVert\nabla
f(\boldsymbol{\theta})-\mathbb{E}[\nabla\bar{f}_{n}(\boldsymbol{\theta})]\right\rVert\leq
Cw_{n}.$ (80)

Using the first bound in (80), we get

$\displaystyle\left\lvert\mathbb{E}[\bar{g}_{n}(\boldsymbol{\theta}_{n})]-f(\boldsymbol{\theta}_{n})\right\rvert^{2}$
$\displaystyle\leq\left(\lvert\mathbb{E}[\bar{g}_{n}(\boldsymbol{\theta}_{n})]-\mathbb{E}[\bar{f}_{n}(\boldsymbol{\theta}_{n})]\rvert+\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\lvert\mathbb{E}[\bar{f}_{n}(\boldsymbol{\theta})]-f(\boldsymbol{\theta})\rvert\right)^{2}$
(81)
$\displaystyle\leq\left\lvert\mathbb{E}[\bar{g}_{n}(\boldsymbol{\theta}_{n})]-\mathbb{E}[\bar{f}_{n}(\boldsymbol{\theta}_{n})]\right\rvert^{2}+Cw_{n+1}\left\lvert\mathbb{E}[\bar{g}_{n}(\boldsymbol{\theta}_{n})]-\mathbb{E}[\bar{f}_{n}(\boldsymbol{\theta}_{n})]\right\rvert+C^{2}w_{n+1}^{2}.$
(82)

Recalling that $\bar{g}_{n}$ and $\bar{f}_{n}$ are uniformly bounded over
$\boldsymbol{\Theta}$ by (A1), using (77) we get

$\displaystyle\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\,\sum_{n=1}^{\infty}w_{n+1}\left\lvert\mathbb{E}[\bar{g}_{n}(\boldsymbol{\theta}_{n})]-f(\boldsymbol{\theta}_{n})\right\rvert^{2}<\infty.$
(83)

An identical argument using the second bound in (80) as well as (79) shows

$\displaystyle\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\,\sum_{n=1}^{\infty}w_{n+1}\left\lVert\mathbb{E}[\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})]-\nabla
f(\boldsymbol{\theta}_{n})\right\rVert^{2}<\infty.$ (84)

Combining the above two bounds and using from Lemma 6.3 give (21).

It remains to show (22) and (20). To this end, suppose the optional condition
in (A4) holds. Then by Lemma 7.2 (v) and Lemma 7.3 (ii), we deduce

$\displaystyle\mathbb{E}\left[\sum_{n=1}^{\infty}w_{n+1}\left(\left|\bar{g}_{n}(\boldsymbol{\theta}_{n})-f(\boldsymbol{\theta}_{n})\right|+\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})-\nabla
f(\boldsymbol{\theta}_{n})\rVert^{2}\right)\right]<\infty.$ (85)

Then (22) follows from Lemma 6.3.

Lastly, recall that by Lemma 7.1, there exists a constant $C>0$ such that for
all $n\geq 1$,

$\displaystyle\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\lvert
f(\boldsymbol{\theta})-\bar{f}_{n}(\boldsymbol{\theta})\right\rvert\leq
Cw_{n}\sqrt{n},\qquad\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\lVert\nabla
f(\boldsymbol{\theta})-\nabla\bar{f}_{n}(\boldsymbol{\theta})\right\rVert\leq
Cw_{n}\sqrt{n}.$ (86)

Then note that

$\displaystyle\lvert\mathbb{E}[\bar{g}_{n}(\boldsymbol{\theta}_{n})]-\bar{f}_{n}(\boldsymbol{\theta}_{n})\rvert^{2}$
$\displaystyle\leq\left(\lvert\mathbb{E}[\bar{g}_{n}(\boldsymbol{\theta}_{n})]-f(\boldsymbol{\theta}_{n})]\rvert+\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\,\lvert\bar{f}_{n}(\boldsymbol{\theta})-f(\boldsymbol{\theta})\rvert\right)^{2}$
(87)
$\displaystyle\leq\left\lvert\mathbb{E}[\bar{g}_{n}(\boldsymbol{\theta}_{n})]-f(\boldsymbol{\theta}_{n})\right\rvert^{2}+Cw_{n+1}\sqrt{n+1}\left\lvert\mathbb{E}[\bar{g}_{n}(\boldsymbol{\theta}_{n})]-\mathbb{E}[\bar{f}_{n}(\boldsymbol{\theta}_{n})]\right\rvert+C^{2}(n+1)w_{n+1}^{2}.$
(88)

Note that by (A4), $w_{n}\sqrt{n}=o(1)$, so we have

$\displaystyle\sum_{n=1}^{\infty}(n+1)w_{n+1}^{3}=\sum_{n=1}^{\infty}\left(w_{n+1}\sqrt{n+1}\right)w_{n+1}^{2}\leq
C^{\prime}\sum_{n=1}^{\infty}w_{n+1}^{2}<\infty$ (89)

for some constant $C^{\prime}>0$. Then again using the fact that $\bar{g}_{n}$
and $\bar{f}_{n}$ are uniformly bounded over $\boldsymbol{\Theta}$, it follows
that

$\displaystyle\sum_{n=1}^{\infty}w_{n+1}\left\lvert\mathbb{E}[\bar{g}_{n}(\boldsymbol{\theta}_{n})]-\bar{f}_{n}(\boldsymbol{\theta}_{n})\right\rvert^{2}<\infty.$
(90)

An identical argument using the second bound in (86) as well as (79) shows

$\displaystyle\sum_{n=1}^{\infty}w_{n+1}\left\lVert\mathbb{E}[\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})]-\nabla\bar{f}_{n}(\boldsymbol{\theta}_{n})\right\rVert^{2}<\infty.$
(91)

Combining the above two bounds gives (20). This completes the proof. ∎

###  7.3. Proof of Theorem 4.1

Next, we prove Theorem 4.1 assuming Lemmas 7.1, 7.2, 7.3, and 7.4.

###### Proof of Theorem 4.1.

Suppose (A1)-(A5) hold. Assume cases C1 or C2 in Theorem 4.1. We will first
show the assertion except the asymptotic stationarity statements.

We first show the first part of (i). By Lemma 7.2 (iii) and Lemma 7.3 (i), we
get

$\displaystyle\mathbb{E}\left[\sum_{n=1}^{\infty}w_{n+1}\left(\left|\bar{g}_{n}(\boldsymbol{\theta}_{n})-\bar{f}_{n}(\boldsymbol{\theta}_{n})\right|+\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})-\nabla\bar{f}_{n}(\boldsymbol{\theta}_{n})\rVert^{2}\right)\right]<\infty.$
(92)

Denote $h^{(1)}_{n}:=\bar{g}_{n}-\bar{f}_{n}$ and
$h^{(2)}_{n}:=\lVert\nabla\bar{g}_{n}-\nabla\bar{f}_{n}\rVert^{2}$. We claim
that
$|h^{(i)}_{n+1}(\boldsymbol{\theta}_{n+1})-h^{(i)}_{n}(\boldsymbol{\theta}_{n})|=O(w_{n+1})$
for $i=1,2$. Then by Lemma 6.3, it holds that
$|h^{(i)}_{n}(\boldsymbol{\theta}_{n})|\rightarrow 0$ almost surely as
$n\rightarrow\infty$ for $i=1,2$.

To show the claim, recall the recursive definitions of $\bar{g}_{n}$ and
$\bar{f}_{n}$ and that
$\boldsymbol{\theta}\mapsto\ell(\mathbf{x},\boldsymbol{\theta})$ is
$R$-Lipschitz continuous by (A1) for each $\mathbf{x}\in\mathfrak{X}$. Also,
$h^{(1)}_{n}(\boldsymbol{\theta})=(1-w_{n})h^{(1)}_{n-1}(\boldsymbol{\theta})+w_{n}(g_{n}(\boldsymbol{\theta})-f_{n}(\boldsymbol{\theta}))$,
where $f_{n}(\cdot)=\ell(\mathbf{x}_{n},\cdot)$ and
$g_{n}\in\operatorname{\textup{{Srg}}}_{L,\rho}^{\mathbb{L}}(f_{n},\boldsymbol{\theta}_{n-1})$
(see Definition 3.1) and $\nabla(g_{n}-f_{n})$ is $L$-Lipschitz continuous and
has norm $\leq\varepsilon_{n}$ at $\boldsymbol{\theta}_{n-1}$. Since the
parameter space $\boldsymbol{\Theta}$ is compact by (A3), it follws that
$\nabla(g_{n}-f_{n})$ is uniformly bounded by some constant $L^{\prime}>0$
over $\boldsymbol{\Theta}$. Then $g_{n}-f_{n}$ is $L^{\prime}$-Lipschitz, and
by induction, $h^{(1)}_{n}$ is also $L$-Lipschitz for $n\geq 1$. Now note that

$\displaystyle|h^{(1)}_{n+1}(\boldsymbol{\theta}_{n+1})-h^{(1)}_{n}(\boldsymbol{\theta}_{n})|$
$\displaystyle\leq|h^{(1)}_{n+1}(\boldsymbol{\theta}_{n+1})-h^{(1)}_{n+1}(\boldsymbol{\theta}_{n})|+|h^{(1)}_{n+1}(\boldsymbol{\theta}_{n})-h^{(1)}_{n}(\boldsymbol{\theta}_{n})|$
(93) $\displaystyle\leq
L^{\prime}\lVert\boldsymbol{\theta}_{n+1}-\boldsymbol{\theta}_{n}\rVert+\left|\left(\bar{g}_{n+1}(\boldsymbol{\theta}_{n})-\bar{g}_{n}(\boldsymbol{\theta}_{n})\right)-\left(\bar{f}_{n+1}(\boldsymbol{\theta}_{n})-\bar{f}_{n}(\boldsymbol{\theta}_{n})\right)\right|$
(94) $\displaystyle\leq
L^{\prime}c^{\prime}w_{n}+w_{n+1}\left|g_{n+1}(\boldsymbol{\theta}_{n})-\bar{g}_{n}(\boldsymbol{\theta}_{n+1})-f_{n+1}(\boldsymbol{\theta}_{n})+\bar{f}_{n}(\boldsymbol{\theta}_{n})\right|$
(95)
$\displaystyle=L^{\prime}\lVert\boldsymbol{\theta}_{n+1}-\boldsymbol{\theta}_{n}\rVert+w_{n+1}\left|\bar{g}_{n}(\boldsymbol{\theta}_{n+1})-\bar{f}_{n}(\boldsymbol{\theta}_{n})\right|,$
(96)

where the third inequality uses Lemma 6.1 (iii) as well as the recursive
definitions of $\bar{g}_{n+1}$ and $\bar{f}_{n+1}$. Then note that
$\ell:\mathcal{X}\times\boldsymbol{\Theta}\rightarrow\mathbb{R}$ is bounded by
(A1), so $\bar{g}_{n}$ and $\bar{f}_{n}$ are also uniformly bounded in $n$.
This shows
$|h^{(1)}_{n+1}(\boldsymbol{\theta}_{n+1})-h^{(1)}_{n}(\boldsymbol{\theta}_{n})|=O(w_{n+1})$.

Next, we verify
$|h^{(2)}_{n}(\boldsymbol{\theta}_{n})-h^{(3)}_{n}(\boldsymbol{\theta}_{n-1})|=O(w_{n+1})$.
Note that by linearity of gradients, $\nabla\bar{g}_{n}$ and
$\nabla\bar{f}_{n}$ satisfy the same recursion as $\bar{g}_{n}$ and
$\bar{f}_{n}$. Moreover, $\nabla
h^{(2)}_{n}(\boldsymbol{\theta})=(1-w_{n})\nabla
h^{(2)}_{n-1}(\boldsymbol{\theta})+w_{n}\nabla(g_{n}-f_{n})(\boldsymbol{\theta})$
and since $\nabla(g_{n}-f_{n})$ is $L$-Lipschitz, so is $\nabla h^{(2)}_{n}$
for all $n\geq 1$. Then by using a similar argument as above,

$\displaystyle\left|\sqrt{h^{(2)}_{n+1}(\boldsymbol{\theta}_{n+1})}-\sqrt{h^{(2)}_{n}(\boldsymbol{\theta}_{n})}\right|$
(97) $\displaystyle\qquad=\lVert\nabla
h_{n+1}^{(2)}(\boldsymbol{\theta}_{n+1})-\nabla
h_{n+1}^{(2)}(\boldsymbol{\theta}_{n})\rVert+\lVert\nabla
h_{n+1}^{(2)}(\boldsymbol{\theta}_{n})-\nabla
h_{n}^{(2)}(\boldsymbol{\theta}_{n})\rVert$ (98) $\displaystyle\qquad\leq
L\lVert\boldsymbol{\theta}_{n+1}-\boldsymbol{\theta}_{n}\rVert+w_{n+1}\lVert\nabla(g_{n+1}-f_{n+1})(\boldsymbol{\theta}_{n})+\nabla\bar{f}_{n}(\boldsymbol{\theta}_{n})-\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})\rVert$
(99) $\displaystyle\qquad\leq Lc^{\prime}w_{n+1}+w_{n+1}\lVert\nabla
h^{(2)}_{n}(\boldsymbol{\theta}_{n})\rVert.$ (100)

But recall that $\nabla h^{(2)}_{n}$ is uniformly bounded by $L^{\prime}$ over
$\boldsymbol{\Theta}$. To conclude, write

$\displaystyle|h^{(2)}_{n}(\boldsymbol{\theta}_{n})-h^{(2)}_{n}(\boldsymbol{\theta}_{n-1})|=\left|\sqrt{h^{(2)}_{n}(\boldsymbol{\theta}_{n})}-\sqrt{h^{(2)}_{n-1}(\boldsymbol{\theta}_{n-1})}\right|\cdot\left|\sqrt{h^{(2)}_{n}(\boldsymbol{\theta}_{n})}+\sqrt{h^{(2)}_{n-1}(\boldsymbol{\theta}_{n-1})}\right|.$
(101)

Noting that $h^{(2)}_{n}$ is uniformly bounded for all $n\geq 1$, we conclude
that
$|h^{(2)}_{n}(\boldsymbol{\theta}_{n})-h^{(3)}_{n}(\boldsymbol{\theta}_{n-1})|=O(w_{n+1})$,
as desired.

Next, we show the first part of (ii). Recall that under (A4) (without the
optional condition), by Lemma 7.1 (se also Lemma B.4), we have
$\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}|\bar{f}_{n}(\boldsymbol{\theta})-f(\boldsymbol{\theta})|\rightarrow
0$ and
$\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\lVert\nabla\bar{f}_{n}(\boldsymbol{\theta})-\nabla
f(\boldsymbol{\theta})\rVert\rightarrow 0$ almost surely as
$n\rightarrow\infty$. Then we by (i), we have

$\displaystyle|\bar{g}_{n}(\boldsymbol{\theta}_{n})-f(\boldsymbol{\theta}_{n})|\leq|\bar{g}_{n}(\boldsymbol{\theta}_{n})-\bar{f}_{n}(\boldsymbol{\theta}_{n})|+\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}|\bar{f}_{n}(\boldsymbol{\theta})-f(\boldsymbol{\theta})|\rightarrow
0\quad\textup{a.s. as $n\rightarrow\infty$}.$ (102)

Similarly, by triangle inequality and (i),

$\displaystyle\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})-\nabla
f(\boldsymbol{\theta}_{n})\rVert\leq\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})-\nabla\bar{f}_{n}(\boldsymbol{\theta}_{n})\rVert+\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\lVert\nabla\bar{f}_{n}(\boldsymbol{\theta})-\nabla
f(\boldsymbol{\theta})\rVert\rightarrow 0\quad\textup{a.s. as
$n\rightarrow\infty$}.$ (103)

Finally, recall that under case C1,
$-\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})$ is in the normal cone of
$\boldsymbol{\Theta}$ at $\boldsymbol{\theta}_{n}$ for each $n\geq 1$. Hence
from the last result, we also deduce that $-\nabla f(\boldsymbol{\theta}_{n})$
lies in the normal cone of $\boldsymbol{\Theta}$ at $\boldsymbol{\theta}_{n}$
asymptotically as $n\rightarrow\infty$.

Next, we show the second parts of (i) and (ii). Assume case C1. Then each
$\boldsymbol{\theta}_{n}$ is an exact minimizer of $\bar{g}_{n}$ over
$\boldsymbol{\Theta}$, so $-\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})$ is in
the normal cone of $\boldsymbol{\Theta}$ at $\boldsymbol{\theta}_{n}$ for each
$n\geq 1$. Hence, $h^{(i)}_{n}\rightarrow 0$ almost surely for $i=2,4$ implies
that both $-\nabla\bar{f}_{n}(\boldsymbol{\theta}_{n})$ and $-\nabla
f(\boldsymbol{\theta}_{n})$ belong to the normal cone of $\boldsymbol{\Theta}$
at $\boldsymbol{\theta}_{n}$ asymptotically almost surely as
$n\rightarrow\infty$.

Lastly, assume case C2. Further assume (A6). Let
$\boldsymbol{\theta}_{\infty}\in\boldsymbol{\Theta}$ be an arbitrary limit
point of the sequence $(\boldsymbol{\theta}_{t})_{t\geq 1}$ and let
$(t_{k})_{k\geq 1}$ be a (random) sequence such that
$\boldsymbol{\theta}_{t_{k}}\rightarrow\boldsymbol{\theta}_{\infty}$ almost
surely as $n\rightarrow\infty$. Since $\boldsymbol{\Theta}\times\mathcal{K}$
is compact, we may choose a further sequence of $(t_{k})_{k\geq 1}$, which we
will denote the same, so that $\kappa_{t_{k}}$ converges to some element
$\kappa_{\infty}\in\mathcal{K}$. Hence
$\bar{g}_{\infty}:=\lim_{k\rightarrow\infty}\bar{g}_{t_{k}}$ is well-defined
almost surely. It is important to note that $\boldsymbol{\theta}_{\infty}$ is
a stationary point of $\bar{g}_{\infty}$ over $\boldsymbol{\Theta}$ by Lemma
7.4. Hence $-\nabla\bar{g}_{\infty}(\boldsymbol{\theta}_{\infty})$ is in the
normal cone of $\boldsymbol{\Theta}$ at $\boldsymbol{\theta}_{\infty}$. But
since we have shown that $h^{(i)}\rightarrow 0$ almost suresly as
$n\rightarrow\infty$ for $i=2,4$, we must have

$\displaystyle\lim_{k\rightarrow\infty}\nabla\bar{f}_{t_{k}}(\boldsymbol{\theta}_{\infty})=\nabla\bar{g}_{\infty}(\boldsymbol{\theta}_{\infty})=\nabla
f(\boldsymbol{\theta}_{\infty}).$ (104)

Since $\boldsymbol{\theta}_{\infty}$ was an arbitrary limit point of
$(\boldsymbol{\theta}_{n})_{n\geq 1}$, this completes the proof. ∎

######  Remark 7.1.

Without appealing to Lemma 7.1 and use the convergence results for the
empirical loss minimization stated in Theorem 4.1, one can directly prove
Theorem 4.1 (ii) using a similar argument as in the proof of (i). This amounts
to use the finite sum (85) in place of (76) and showing the bounds
$|h^{(i)}_{n+1}(\boldsymbol{\theta}_{n+1})-h^{(i)}_{n}(\boldsymbol{\theta}_{n})|=O(w_{n+1})$
for $i=3,4$, where $h^{(3)}=\bar{g}_{n+1}-f$ and
$h^{(4)}_{n}=\lVert\nabla\bar{g}_{n}-\nabla f\rVert^{2}$. However, (85) holds
under the additional optional condition in (A4), which essentially says
$\sum_{n=1}^{\infty}w_{n}^{2}\sqrt{n}<\infty$. While this condition is
standard in the literature (see, e.g., [35, (E)]), [30, (M2)], and [32, (A3)],
our proof above avoids it and only relies on the more standard square
summability condition $\sum_{n=1}^{\infty}(\log n)w^{2}<\infty$ (see (A4)).
However, the stronger condition $\sum_{n=1}^{\infty}w_{n}^{2}\sqrt{n}<\infty$
as well as the finite sum (85) is crucial in our proof of Theorem 4.3.

###  7.4. Proof of Theorem 4.3

Next, we prove Theorem 4.3 assuming Lemmas 7.2, 7.3, and 7.5.

###### Proof of Theorem 4.3.

We first show (i). Note that (23) follows immediately from Lemmas 7.5 and 6.3.
Next we show (24). By Cauchy-Schwarz inequality, for all
$\boldsymbol{\theta}\in\boldsymbol{\Theta}$,

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

It follows that for all $n\geq 1$,

$\displaystyle\left|-\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\langle\nabla\bar{f}_{n}(\boldsymbol{\theta}_{n}),\,\frac{\boldsymbol{\theta}-\boldsymbol{\theta}_{n}}{\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert}\right\rangle\right|\leq\left|-\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\langle\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n}),\,\frac{\boldsymbol{\theta}-\boldsymbol{\theta}_{n}}{\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert}\right\rangle\right|+\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})-\nabla\bar{f}_{n}(\boldsymbol{\theta}_{n})\rVert.$
(106)

On the other hand, by Lemmas 7.5 and 7.3 (i), we have

$\displaystyle\sum_{n=1}^{\infty}w_{n+1}\left[\left|-\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\langle\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n}),\,\frac{\boldsymbol{\theta}-\boldsymbol{\theta}_{n}}{\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert}\right\rangle\right|+\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})-\nabla\bar{f}_{n}(\boldsymbol{\theta}_{n})\rVert^{2}\right]<\infty.$
(107)

Then by Lemma 6.3, we have

$\displaystyle\min_{1\leq k\leq
n}\left[\left|-\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\langle\nabla\bar{g}_{k}(\boldsymbol{\theta}_{k}),\,\frac{\boldsymbol{\theta}-\boldsymbol{\theta}_{k}}{\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{k}\rVert}\right\rangle\right|+\lVert\nabla\bar{g}_{k}(\boldsymbol{\theta}_{k})-\nabla\bar{f}_{k}(\boldsymbol{\theta}_{k})\rVert^{2}\right]=O\left(\left(\sum_{k=1}^{n}w_{k}\right)^{-1}\right).$
(108)

Let $t_{n}\in\\{1,\dots,n\\}$ for $n\geq 1$ be such that the minimum above is
achieved. Namely, denoting the term in the minimum above by $A_{k}$, we have
$A_{t_{n}}=O\left(\left(\sum_{k=1}^{t_{n}}w_{k}\right)^{-1}\right)$. Since all
terms in $A_{k}$ are nonnegative, it follows that there exists a constant
$c_{1},c_{2}>0$ such that fro all $n\geq 1$,

$\displaystyle\left|-\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\langle\nabla\bar{g}_{t_{n}}(\boldsymbol{\theta}_{t_{n}}),\,\frac{\boldsymbol{\theta}-\boldsymbol{\theta}_{t_{n}}}{\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{t_{n}}\rVert}\right\rangle\right|\leq\frac{c_{1}}{\sum_{k=1}^{t_{n}}w_{k}},\quad\lVert\nabla\bar{g}_{t_{n}}(\boldsymbol{\theta}_{t_{n}})-\nabla
f(\boldsymbol{\theta}_{t_{n}})\rVert\leq\frac{c_{2}}{\sqrt{\sum_{k=1}^{t_{n}}w_{k}}}.$
(109)

Hence from (106), it follows that there exists some constant $c_{3}>0$ such
that for all $n\geq 1$,

$\displaystyle\left|-\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\langle\nabla\bar{f}_{t_{n}}(\boldsymbol{\theta}_{t_{n}}),\,\frac{\boldsymbol{\theta}-\boldsymbol{\theta}_{t_{n}}}{\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{t_{n}}\rVert}\right\rangle\right|\leq\frac{c_{3}}{\sqrt{\sum_{k=1}^{t_{n}}w_{k}}}.$
(110)

This completes the proof of (i).

It remains to show (25) in (ii). Assume the optional condition in (A4) holds.
The argument is similar as before. We first observe

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

It follows that for all $n\geq 1$,

$\displaystyle\left|-\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\langle\nabla
f(\boldsymbol{\theta}_{n}),\,\frac{\boldsymbol{\theta}-\boldsymbol{\theta}_{n}}{\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert}\right\rangle\right|\leq\left|-\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\langle\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n}),\,\frac{\boldsymbol{\theta}-\boldsymbol{\theta}_{n}}{\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert}\right\rangle\right|+\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})-\nabla
f(\boldsymbol{\theta}_{n})\rVert.$ (112)

Also, by Lemmas 7.5 and 7.3 (ii), we have

$\displaystyle\sum_{n=1}^{\infty}w_{n+1}\left[\left|-\inf_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\langle\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n}),\,\frac{\boldsymbol{\theta}-\boldsymbol{\theta}_{n}}{\lVert\boldsymbol{\theta}-\boldsymbol{\theta}_{n}\rVert}\right\rangle\right|+\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})-\nabla
f(\boldsymbol{\theta}_{n})\rVert^{2}\right]<\infty.$ (113)

The rest of the argument is identical to the previous case. ∎

###  7.5. Proofs of Corollaries 4.4 and 4.5

Lastly in this section, we prove Corollaries 4.4 and 4.5.

###### Proof of Corollary 4.4.

Suppose the iterates $\boldsymbol{\theta}_{n}$ are all in the interior of
$\boldsymbol{\Theta}$. Then in Theorem 4.3 (23), we may choose
$\boldsymbol{\theta}\in\boldsymbol{\Theta}$ in a way that
$\boldsymbol{\theta}-\boldsymbol{\theta}_{n}$ is a positive scalar multiple of
$\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})$. Then (23) reduces to the first
bound in (26). The same argument shows the second bound in (26). For the last
bound in (26), we use (84) we derived in the proof of Theorem 4.2. Combining
it with Lemma 7.5 and using Fubini’s theorem, we get

$\displaystyle\sum_{n=1}^{\infty}w_{n+1}\left(\mathbb{E}\left[\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})\rVert\right]+\left\lVert\mathbb{E}[\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})]-\nabla\bar{f}_{n}(\boldsymbol{\theta}_{n})\right\rVert^{2}\right)<\infty.$
(114)

Then by Lemma 6.3, we have

$\displaystyle\min_{1\leq k\leq
n}\left[\mathbb{E}\left[\lVert\nabla\bar{g}_{n}(\boldsymbol{\theta}_{n})\rVert\right]+\left\lVert\nabla\mathbb{E}[\bar{g}_{k}(\boldsymbol{\theta}_{k})]-\nabla
f(\boldsymbol{\theta}_{k})\right\rVert^{2}\right]=O\left(\left(\sum_{k=1}^{n}w_{k}\right)^{-1}\right)$
(115)

Choose a subsequence $t_{n}$ along which the above minimum is achieved for all
$n\geq 1$. Then we have

$\displaystyle\mathbb{E}\left[\lVert\nabla\bar{g}_{t_{n}}(\boldsymbol{\theta}_{t_{n}})\rVert\right]=O\left(\left(\sum_{k=1}^{t_{n}}w_{k}\right)^{-1}\right),\quad\left\lVert\nabla\mathbb{E}[\bar{g}_{t_{n}}(\boldsymbol{\theta}_{t_{n}})]-\nabla
f(\boldsymbol{\theta}_{t_{n}})\right\rVert=O\left(\left(\sum_{k=1}^{t_{n}}w_{k}\right)^{-1/2}\right).$
(116)

Therefore by using Jensen’s inequality and the above bound, we deduce

$\displaystyle\left\lVert\nabla
f(\boldsymbol{\theta}_{t_{n}})\right\rVert\leq\mathbb{E}\left[\left\lVert\nabla\bar{g}_{t_{n}}(\boldsymbol{\theta}_{t_{n}})\right\rVert\right]+\left\lVert\nabla\mathbb{E}[\bar{g}_{t_{n}}(\boldsymbol{\theta}_{t_{n}})]-\nabla
f(\boldsymbol{\theta}_{t_{n}})\right\rVert=O\left(\left(\sum_{k=1}^{t_{n}}w_{k}\right)^{-1/2}\right).$
(117)

This completes the proof. ∎

###### Proof of Corollary 4.5 .

First, suppose $w_{n}=n^{-1/2}(\log n)^{\delta}$ for some $\delta>1$. Then the
upper bound on the rate of convergence in Theorem 4.3 is of order
$O(n^{-1/2}(\log n)^{\delta})$. Then one can conclude by using the fact that
$n\geq\varepsilon^{-2}(3\log\varepsilon^{-1})^{2\delta}$ implies
$n^{-1/2}(\log n)^{\delta}\leq\varepsilon$ for all sufficiently small
$\varepsilon>0$. Indeed, assuming
$n\geq\varepsilon^{-2}(3\log\varepsilon^{-1})^{2\delta}$,

$\displaystyle n^{-1/2}(\log
n)^{\delta}\leq\frac{\varepsilon}{(3\log\varepsilon^{-1})^{\delta}}(2\log\varepsilon^{-1}+2\delta\log(3\log\varepsilon^{-1}))^{\delta}$
(118)

and the last expression is at most $\varepsilon$ for all sufficiently small
$\varepsilon>0$. This shows (i). A similar argument using Theorem 4.3 (ii) as
well as Corollary 4.4 shows (ii). ∎

##  8\. Proof of Lemmas 7.1, 7.2, and 7.3

For each $n\geq 0$, we let $\mathcal{F}_{n}$ denote the $\sigma$-algebra
generated by the history of underlying Markov chain $Y_{0},Y_{1},\dots,Y_{n}$
as well as the possible randomness in Algorithm 1 up to iteration $n$. Note
that $\mathbf{x}_{n}=\varphi(Y_{n})$ for $n\geq 1$ by (A2).

We first prove Lemma 7.1 below. Our proof uses an auxiliary lemma, Lemma B.4
in the appendix.

###### Proof of Lemma 7.1.

Recall that $\mathbf{x}_{k}=\varphi(Y_{k})$ under (A2). Let $\pi_{k}$ denote
the distribution of $Y_{k}$. Let
$M:=\sup_{\mathbf{x},\boldsymbol{\theta}}\lVert\psi(\mathbf{x},\boldsymbol{\theta})\rVert<\infty$.
Note that, by a change of measure,

$\displaystyle\mathbb{E}\left[\psi(\mathbf{x}_{n},\boldsymbol{\theta})\right]$
$\displaystyle=\sum_{\mathbf{y}\in\Omega}\psi(\varphi(\mathbf{y}),\boldsymbol{\theta})\,P^{n}(Y_{0},\mathbf{y})$
(119)
$\displaystyle=\sum_{\mathbf{y}\in\Omega}\psi(\varphi(\mathbf{y}),\boldsymbol{\theta})\,\pi(\mathbf{y})+\sum_{\mathbf{y}\in\Omega}\psi(\varphi(Y_{0}),\boldsymbol{\theta})(P^{n}(Y_{0},\mathbf{y})-\pi(\mathbf{y}))$
(120)
$\displaystyle=\bar{\psi}(\boldsymbol{\theta})+\sum_{\mathbf{y}\in\Omega}\psi(\varphi(Y_{0}),\boldsymbol{\theta})\,(P^{n}(Y_{0},\mathbf{y})-\pi(\mathbf{y})).$
(121)

By triangle inequality, it follows that

$\displaystyle\left\lVert\mathbb{E}\left[\psi(\mathbf{x}_{n},\boldsymbol{\theta})\right]-\bar{\psi}(\boldsymbol{\theta})\right\rVert$
$\displaystyle\leq\sum_{\mathbf{y}\in\Omega}\left\lVert\psi(\varphi(Y_{0}),\boldsymbol{\theta})\right\rVert\,\,\left|P^{n}(Y_{0},\mathbf{y})-\pi(\mathbf{y})\right|$
(122) $\displaystyle\leq
M\sum_{y\in\Omega}\left|P^{n}(Y_{0},\mathbf{y})-\pi(\mathbf{y})\right|$ (123)
$\displaystyle\leq 2M\lVert P^{n}(Y_{0},\cdot)-\pi\rVert_{TV},$ (124)

where the last inequality follows from the relation
$2\lVert\mu-\nu\rVert_{TV}=\sum_{x}|\mu(x)-\nu(x)|$ where $\mu,\nu$ are
probability distributions on the same sample space (see [31, Prop. 4.2]).

Now recall that under the hypothesis in (A4),
$w^{n}_{k}=w_{k}\prod_{i=k+1}^{n}(1-w_{i})$ is non-decreasing in
$k\in\\{1,\dots,n\\}$. So, $w^{n}_{1}\leq\dots\leq w^{n}_{n}=w_{n}$. Then
using (A2), we have

$\displaystyle\left\lVert\bar{\psi}(\boldsymbol{\theta})-\mathbb{E}\left[\bar{\psi}_{n}(\mathbf{x},\boldsymbol{\theta})\right]\right\rVert$
$\displaystyle\leq\left\lVert\sum_{k=1}^{n}\left(\bar{\psi}(\boldsymbol{\theta})-\mathbb{E}[\psi(\mathbf{x}_{k},\boldsymbol{\theta})]\right)w^{n}_{k}\right\rVert$
(125)
$\displaystyle\leq\sum_{k=1}^{n}\left\lVert\bar{\psi}(\boldsymbol{\theta})-\mathbb{E}[\psi(\mathbf{x}_{k},\boldsymbol{\theta})]\right\rVert
w^{n}_{k}$ (126) $\displaystyle\leq 2M\sum_{k=1}^{n}\lVert
P^{k}(Y_{0},\cdot)-\pi\rVert_{TV}w^{n}_{k}$ (127) $\displaystyle\leq
2Mw_{n}\sum_{k=1}^{n}\sup_{\mathbf{y}\in\mathfrak{X}}\,\lVert
P^{k}(\mathbf{y},\cdot)-\pi\rVert_{TV}\leq\frac{2L\lambda w_{n}}{1-\lambda}.$
(128)

This shows the first bound in (72).

Finally, for $d=1$, the second inequality in (72) as well as the last part ff
the statement are direct consequenes of Lemma B.4, nothing that
$w^{n}_{1}\leq\dots\leq w^{n}_{n}=w_{n}$. Applying this to each of the $d$
coordinates of $\psi$ and applying triangle inequality will imply the
assertion. Namely, write $\psi=(\psi^{(1)},\dots,\psi^{(d)})^{T}$. Then noting
that $\psi,\bar{\psi}_{n}$ are uniformly bounded by some constant $M>0$, we
have

$\displaystyle\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\lVert\bar{\psi}(\boldsymbol{\theta})-\bar{\psi}_{n}(\boldsymbol{\theta})\right\rVert^{2}$
$\displaystyle\leq\sum_{i=1}^{d}\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\lvert\bar{\psi}^{(i)}(\boldsymbol{\theta})-\bar{\psi}^{(i)}_{n}(\boldsymbol{\theta})\right\rvert^{2}\leq
2M\sum_{i=1}^{d}\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left\lvert\bar{\psi}^{(i)}(\boldsymbol{\theta})-\bar{\psi}^{(i)}_{n}(\boldsymbol{\theta})\right\rvert.$
(129)

Note that by Lemma B.4, the expectation of each summand in the last expression
is of order $O(w_{n}\sqrt{n})$. Hence the second inequality in (72) follows.
Lastly, the summands in the last expresssion above converges to zero almost
surely if $w^{n}_{n}\sqrt{n}=O(1/(\log n)^{1+\varepsilon})$ for some
$\varepsilon>0$ by Lemma B.4. Note that $w^{n}_{n}=w_{n}$. Hence the left hand
side above converges to zero almost surely $w_{n}\sqrt{n}=O(1/(\log
n)^{1+\varepsilon})$. This completes the proof. ∎

To handle the issue of dependence in signals, we adopt the strategy developed
in [30] in order to handle a similar issue for vector-valued signals (or
matrix factorization). The key insight in [30] is that, while the 1-step
conditional distribution $P(Y_{t-1},\cdot)$ may be far from the stationary
distribution $\pi$, the $N$-step conditional distribution
$P^{N}(Y_{t-N},\cdot)$ is exponentially close to $\pi$ under mild conditions.
Hence we can condition much early on – at time $t-N$ for some suitable
$N=N(t)$. Then the Markov chain runs $N+1$ steps up to time $t+1$, so if $N$
is large enough for the chain to mix to its stationary distribution $\pi$,
then the distribution of $Y_{t+1}$ conditional on $\mathcal{F}_{t-N}$ is close
to $\pi$. The error of approximating the stationary distribution by the $N+1$
step distribution can be controlled using total variation distance and Markov
chain mixing bound. We refine the analysis for [30, Lem. 12] and obtain a more
improved version below.

######  Proposition 8.1.

Assumptions (A1)-(A4) hold. Let $(\boldsymbol{\theta}_{n})_{n\geq 1}$ be an
output of Algorithm and denote $L=\lVert\ell(\cdot,\cdot)\rVert_{\infty}$.
Then there exists a constant $C>0$ such that for all $0\leq M\leq N\leq n$, it
holds that

$\displaystyle\mathbb{E}\left[\mathbb{E}\left[\ell(\mathbf{x}_{n+1},\boldsymbol{\theta}_{n-N})-\bar{f}_{n}(\boldsymbol{\theta}_{n-N})\,\bigg{|}\,\mathcal{F}_{n-N}\right]^{+}\right]\leq
Cw_{n-N}+2L\lambda^{N+1},$ (130)

where $\lambda\in[0,1)$ is the exponential mixing rate of the underlying
Markov chain in (A2).

###### Proof.

Fix $\mathbf{y}\in\Omega$ and suppose $Y_{n-N}=\mathbf{y}$. By the Markov
property, the distribution of $Y_{n+1}$ conditional on $\mathcal{F}_{n-N}$
equals $P^{N+1}(Y_{n-N},\cdot)$, where $P$ denotes the transition kernel of
the chain $(Y_{n})_{n\in\mathbb{N}}$. Then by the Markov property, the
distribution of $Y_{n+1}$ conditional on $\mathcal{F}_{n-N}$ equals
$P^{N+1}(\mathbf{y},\cdot)$, where $P$ denotes the transition kernel of the
chain $(Y_{n})_{t\in\mathbb{N}}$. Denote
$L:=\lVert\ell(\cdot,\cdot)\rVert_{\infty}$. Denote
$\boldsymbol{\theta}:=\boldsymbol{\theta}_{n-N}$, which is deterministic with
respect to $\mathcal{F}_{n-N}$.

We first claim that

$\displaystyle\left|\mathbb{E}\left[\ell(\mathbf{x}_{n+1},\boldsymbol{\theta})\,\bigg{|}\,\mathcal{F}_{n-N}\right]-f(\boldsymbol{\theta})\right|\leq
2L\lVert P^{N+1}(\mathbf{y},\cdot)-\pi\rVert_{TV}\leq 2L\lambda^{N+1}.$ (131)

Then this would yield

$\displaystyle\mathbb{E}\left[\ell(\mathbf{x}_{n+1},\boldsymbol{\theta})\,\bigg{|}\,\mathcal{F}_{n-N}\right]\leq
f(\boldsymbol{\theta})+\left|\mathbb{E}\left[\ell(\mathbf{x}_{n+1},\boldsymbol{\theta})\,\bigg{|}\,\mathcal{F}_{n-N}\right]-f(\boldsymbol{\theta})\right|\leq
f(\boldsymbol{\theta})+2L\lambda^{N+1}.$ (132)

To justify (131), note that by a change of measure,

$\displaystyle\mathbb{E}\left[\ell(\mathbf{x}_{n+1},\boldsymbol{\theta})\,\bigg{|}\,\mathcal{F}_{n-N}\right]$
$\displaystyle=\sum_{\mathbf{y}^{\prime}\in\Omega}\ell(\varphi(\mathbf{y}^{\prime}),\boldsymbol{\theta})\,P^{N+1}(\mathbf{y},\mathbf{y}^{\prime})$
(133)
$\displaystyle=\sum_{\mathbf{y}^{\prime}\in\Omega}\ell(\varphi(\mathbf{y}^{\prime}),\boldsymbol{\theta})\,\pi(\mathbf{y}^{\prime})+\sum_{\mathbf{y}^{\prime}\in\Omega}\ell(\varphi(\mathbf{y}^{\prime}),\boldsymbol{\theta})(P^{N+1}(\mathbf{y},\mathbf{y}^{\prime})-\pi(\mathbf{y}^{\prime}))$
(134)
$\displaystyle=f(\boldsymbol{\theta})+\sum_{\mathbf{y}^{\prime}\in\Omega}\ell(\varphi(\mathbf{y}^{\prime}),\boldsymbol{\theta})(P^{N+1}(\mathbf{y},\mathbf{y}^{\prime})-\pi(\mathbf{y}^{\prime})).$
(135)

Then the claim follows by bouding
$\ell(\varphi(\mathbf{y}^{\prime}),\boldsymbol{\theta})\leq\lVert\cdot,\boldsymbol{\theta}\rVert_{\infty}$
for all $\mathbf{y}^{\prime}\in\Omega$ and
$\sum_{\mathbf{y}^{\prime}\in\Omega}\left|P^{N+1}(\mathbf{y},\mathbf{y}^{\prime})-\pi(\mathbf{y}^{\prime})\right|=2\lVert
P^{N+1}(\mathbf{y},\cdot)-\pi\rVert_{TV}$.

Next, we analyze
$\mathbb{E}\left[-\bar{f}_{n}(\boldsymbol{\theta})\,|\,\mathcal{F}_{n-N}\right]$
using a different approach. Namely, we decompose the times interval $[0,n]$
into two intervals $[1,n-N]$ and $[n-N+1,t]$ and consider the trajectory
$(\mathbf{x}_{k})_{1\leq k\leq n}$ restricted onto each of these intervals.
Since we are conditioning on $\mathcal{F}_{n-N}$, the trajectory
$(\mathbf{x}_{k})_{1\leq k\leq n-N}$ is fully observed and may be far from the
average behavior. The distribution of $\mathbf{x}_{k}$ then becomes close to
the stationary distribution $\pi$ during the second interval of length at an
exponential mixing rate of $\lambda$ (see (A2)). We first compare the
trajectory on the second interval using the Markov chain mixing, and then
compare the mean of the trajectory of the first interval with the stationary
trajectory by integrating out the conditioning on $\mathcal{F}_{n-N}$. See
Figure 1 for illustration.

Figure 1.  Illustration of the decomposition of the process used to bound the
positive variation.

We proceed with the sketch given above. Denote
$\hat{f}_{n:N}:=\sum_{k=t-N+1}^{n}\ell(\mathbf{x}_{k},\cdot)w^{n}_{k}$. Write

$\displaystyle\mathbb{E}\left[-\bar{f}_{n}(\boldsymbol{\theta})\,\bigg{|}\,\mathcal{F}_{n-N}\right]$
$\displaystyle=-\mathbb{E}\left[\sum_{k=1}^{n}\ell(\mathbf{x}_{k},\boldsymbol{\theta})w^{n}_{k}\,\bigg{|}\,\mathcal{F}_{n-N}\right]=-\left(\sum_{k=1}^{n-N}\ell(\mathbf{x}_{k},\boldsymbol{\theta})w^{n}_{k}\right)-\mathbb{E}\left[\hat{f}_{n:N}\,\bigg{|}\,\mathcal{F}_{n-N}\right]$
(136)

Next, denote $C_{N}=\sum_{k=n-N+1}^{n}w^{n}_{k}$. Using Markov property, we
have

$\displaystyle\left|C_{N}f(\boldsymbol{\theta})-\mathbb{E}\left[\hat{f}_{n:N}(\boldsymbol{\theta})\,\bigg{|}\,\mathcal{F}_{n-N}\right]\right|$