\section{Stochastic gradient descent with idealized \texttt{MICE}}
\label{sc:ideal_mice}

This section describes the stochastic gradient descent with Gaussian noise in gradients satisfying Assumption \ref{as:grad_rel_err}, which we call stochastic gradient descent with idealized \texttt{MICE}.
We focus on the case where Assumptions \ref{as:lipschitz} and \ref{as:strongly.convex} are also satisfied.
This method is not a practical stochastic optimization method since we use the exact mean gradient to generate a noisy gradient estimate.
The goal of this artificial method is to validate Proposition \ref{prp:convergence.mice} and Corollary \ref{cor:as_rate_conv} on, respectively, the $L^2$ and a.s. convergence rates of \texttt{SGD-MICE}.
If the bias of \texttt{MICE} and the errors in gradient norm and variances are negligible, \texttt{SGD-MICE} should behave like idealized \texttt{SGD-MICE}.
Thus, by comparing \texttt{SGD-MICE} with idealized \texttt{SGD-MICE}, we can verify how close \texttt{SGD-MICE} is to achieving the $L^2$ and a.s. convergences in Proposition \ref{prp:convergence.mice} and Corollary \ref{cor:as_rate_conv}.

To implement idealized \texttt{SGD-MICE}, we add Gaussian noise to the true gradient and use it as an estimate of the gradient in \texttt{SGD} in such a way that Assumption \ref{as:grad_rel_err} is satisfied.
For a true gradient $\nabla_{\bs{\xi}} F(\bs{\xi}_k)$ and a covariance matrix $\bs{\Sigma}$ of $\nabla_{\bs{\xi}} f(\bs{\xi}, \bs{\theta})$, we define the mean gradient estimate $\tilde{\bs{\upsilon}}$ as
\begin{equation}
  \tilde{\bs{\upsilon}} = \nabla_{\bs{\xi}} F(\bs{\xi}) +
  \frac{\epsilon\|\nabla_{\bs{\xi}} F(\bs{\xi})\|}{\sqrt{\hbox{\textrm{tr}}\mskip2mu(\bs{\Sigma})}}
  \bs{\chi},
\end{equation}
where $\bs{\chi} \sim \cl{N}(\bs{0}, \bs{\Sigma})$.
Taking the expectation of $\tilde{\bs{\upsilon}}$ conditioned on $\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}$,
\begin{equation}
\bb{E}[\tilde{\bs{\upsilon}}_k |\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}] = \nabla_{\bs{\xi}} F(\bs{\xi}_k),
\end{equation}
and taking the expectation of mean squared error conditioned on $\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}$ we see that
\begin{equation}
 \bb{E}[\norm{\nabla_{\bs{\xi}} F(\bs{\xi}_k)-\tilde{\bs{\upsilon}}_k}{}^2|\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}]
 =\epsilon^2\norm{\nabla_{\bs{\xi}} F(\bs{\xi}_k)}{}^2,
 \end{equation}
 hence, Assumption \ref{as:grad_rel_err} holds.
 Consequently, idealized \texttt{SGD-MICE} $L^2$ and a.s. convergences follow Proposition \ref{prp:convergence.mice} and Corollary \ref{cor:as_rate_conv}, respectively.

\section{Multi-iteration stochastic optimizers} \label{sc:algorithms}
In this section, we present the detailed algorithms for the multi-iteration stochastic optimizers using the \texttt{MICE} estimator for the mean gradient.
In Algorithms \ref{alg:sgd}, \ref{alg:adam}, we respectively describe the pseudocodes for \texttt{SGD-MICE} and \texttt{Adam-MICE}.

\begin{algorithm}[h]
\setstretch{1.5}
\caption{Pseudocode for the \texttt{SGD-MICE} with fixed step-size.
\texttt{SGD-MICE} requires an unbiased estimator of the true gradient, $\nabla f$; a starting point, $\bs{\xi}_0$; a step-size, $\eta$; and a distribution from which $\bs{\theta}$ can be sampled, $\pi$.
}
\label{alg:sgd}
\begin{algorithmic}[1]
\Procedure{SGD-MICE}{$\nabla_{\bs{\xi}} f$, $\bs{\xi}_0, \eta, \pi$}
\For{$k=0,1,2...K$}
\State Evaluate $\nabla_{\bs{\xi}} \cl{F}_k$ using Algorithm \ref{alg:mice}
\State $\bs{\xi}_{k+1} \gets \bs{\xi}_k - \eta \nabla_{\bs{\xi}} \cl{F}_k$
\EndFor
\State \textbf{return} optimal point approximation $\bs{\xi}_{K+1}$
\EndProcedure
\end{algorithmic}
\end{algorithm}

\begin{algorithm}[h]
\setstretch{1.5}
\caption{Pseudocode for the \texttt{Adam-MICE} with fixed step-size.
\texttt{Adam-MICE} requires an unbiased estimator of the true gradient, $\nabla_{\bs{\xi}} f$; a starting point, $\bs{\xi}_0$; a step-size, $\eta$; and a distribution from which $\bs{\theta}$ can be sampled, $\pi$.
Moreover, \texttt{Adam-MICE} requires the constants $\beta_1$, $\beta_2$, and $\epsilon_{\text{Adam}}$.
We use the values recommended by \cite{kingma2014adam}, $\beta_1=0.9$, $\beta_2=0.999$, and $\epsilon_{\text{Adam}}=10^{-8}$.
}
\label{alg:adam}
\begin{algorithmic}[1]
\Procedure{Adam-MICE}{$\nabla_{\bs{\xi}} f$, $\bs{\xi}_0, \eta, \pi$}
\State Initialize $\bs{m}_0$ and $\bs{v}_0$ as zero-vectors
\For{$k=1,2...K$}
\State Evaluate $\nabla_{\bs{\xi}} \cl{F}_k$ using Algorithm \ref{alg:mice}
\State $\bs{m}_k \gets \beta_1 \bs{m}_{k-1} + (1 - \beta_1) \nabla_{\bs{\xi}} \cl{F}_k$
\State $\bs{v}_k \gets \beta_2 \bs{v}_{k-1} + (1 - \beta_2) \nabla_{\bs{\xi}} \cl{F}_k^2$ \Comment{The gradient estimates are squared element-wise}
\State $\bs{\hat{m}}_k \gets \bs{m}_{k-1} / (1- \beta_1^k)$
\State $\bs{\hat{v}}_k \gets \bs{v}_{k-1} / (1- \beta_2^k)$
\State $\bs{\xi}_{k+1} \gets \bs{\xi}_k - \eta \bs{\hat{m}}_k / (\sqrt{\bs{\hat{v}}_k} + \epsilon_{\text{Adam}})$
\EndFor
\State \textbf{return} optimal point approximation $\bs{\xi}_{K+1}$
\EndProcedure
\end{algorithmic}
\end{algorithm}

\section{Acknowledgments}

This work was partially supported by the KAUST Office of Sponsored Research (OSR) under Award numbers URF$/1/2281-01-01$, URF$/1/2584-01-01$ in the KAUST Competitive Research Grants Program Round 8, the Alexander von Humboldt Foundation, and Coordination for the Improvement of Higher Education Personnel (CAPES).

Last but not least, we want to thank Dr. S\"{o}ren Wolfers and Prof. Jesper Oppelstrup. They both provided us with valuable ideas and constructive comments.

\section{Conclusion}

We propose the \emph{Multi-Iteration Stochastic Optimizers}, a novel class of first-order stochastic optimizers. To build this class, we presented an estimator of the mean gradient \emph{Multi-Iteration stochastiC Estimator}---\texttt{MICE}---based on successive control variates along the path of iterations. \texttt{MICE} uses a hierarchy based on previously computed sample gradients from iterations in the optimization path.
At each iteration, \texttt{MICE} samples new gradients to achieve a desired relative tolerance on the gradient's statistical error with minimal computational cost.
We coupled \texttt{MICE} with techniques to automatically update its hierarchy. These updates selectively perform restarting, dropping iterations, or clipping the hierarchy to minimize the computational work.
Furthermore, to robustify the use of \texttt{MICE}, we use a jackknife resampling technique to estimate the norm of the mean gradient. This norm estimate is useful for picking the optimizer's step length, thus controlling \texttt{MICE}'s relative statistical error.
Due to its nonintrusive nature, it is relatively simple to couple \texttt{MICE} with many existing first-order optimization methods, resulting in a new family of optimizers.

An idealized theoretical analysis for Stochastic Gradient Descent coupled with the \texttt{MICE} estimator (\texttt{SGD-MICE}) is provided for $L$-smooth and strongly-convex functions.
In this context, we motivate the relative error control's usefulness in the mean gradient estimates, proving exponential convergence in $L^2$ sense with respect to the number of iterations using a constant step-size.
We also estimate the resulting mean computational cost in terms of the number of sampled gradients.

To test the performance of this new family of optimizers, we present three numerical examples.
In all three examples, with either finite and infinite populations, we use the same parameters for \texttt{MICE}, changing only the tolerance on the statistical error.
The first example is a quadratic function with a stochastic Hessian, used to verify our analysis. We compare \texttt{SGD-MICE} with the expected convergence rates in terms of iterations and gradient evaluations.
The second example is an adaptation of the Rosenbrock function to the stochastic setting, where we use \texttt{Adam-MICE} and \texttt{Adam} to perform optimization.
The third example consists of training a logistic regression model over three different datasets, one of which is $11 \times 10^6$, showing that \texttt{SGD-MICE} can outperform common optimization methods, such as \texttt{SGD}, \texttt{SAG}, \texttt{SAGA}, \texttt{SRVG}, and \texttt{SARAH}, in these supervised machine learning problems with large datasets.

Finally, although in this work we only address unconstrained optimization,
one can use similar $\texttt{MICE}$ estimators in constrained optimization, reusing all the
standard techniques, for instance, projected gradients, active sets, among others.

\section{Introduction}

We focus on the stochastic optimization problem of minimizing the objective function $\bb{E} [f(\bs{\xi},\bs{\theta})|\bs{\xi}],$ where $f$ is a given real-valued function, $\bs{\xi}$ is the design variable vector, $\bs{\theta}$ is a random vector, and $\bb{E}[\cdot|\bs{\xi}]$ is the expectation conditioned on $\bs{\xi}$.
Stochastic optimization problems \cite{Marti,Pardalos,Birge2011} are relevant to different fields, such as Machine Learning \cite{Lan2020}, Stochastic Optimal Control \cite{wallace2005applications,fleming2012deterministic}, Computational Finance \cite{Ziemba1975,Conejo,bayer2018pricing}, Economics \cite{Chang}, Insurance \cite{Azcue}, Communication Networks \cite{ding2012stochastic}, Queues and Supply Chains \cite{Yao}, and  Bayesian Optimal Design of Experiments \cite{ryan2016review,carlon2020nesterov}, among many others.

In the same spirit and inspired by the work by Heinrich \cite{heinrich2000multilevel} and Giles \cite{giles:MLMC} on \emph{Multilevel Monte Carlo} methods, we propose the \emph{Multi-Iteration stochastiC Estimator}---\texttt{MICE}---to obtain a computationally efficient and affordable approximation of the mean gradient, $\nabla_{\bs{\xi}} \bb{E}[f(\bs{\xi},\bs{\theta})|\bs{\xi}]$, which may be coupled with any {first order} stochastic optimizer in a non-intrusive fashion.
Combining \texttt{MICE} with any stochastic optimizer furnishes \emph{Multi-Iteration Stochastic Optimizers}, a novel class of efficient and robust stochastic optimizers.
In this class of stochastic optimizers, the mean gradient estimator's relative variance is controlled using successive control variates based on previous iterations' available information.
This procedure results in a more accurate yet affordable estimation of the mean gradient.
In approximating the mean gradient, \texttt{MICE} constructs a hierarchy of iterations and performs control variates for every pair of nested elements of this hierarchy.
As the stochastic optimization evolves, we increase the number of samples along the hierarchy while keeping the previously sampled gradients, i.e., we use gradient information from previous iterations to reduce the variance in the current gradient estimate, which is a crucial feature to make \texttt{MICE} competitive.
We design \texttt{MICE} to achieve a given relative statistical error for the mean gradient with minimum additional gradient sampling cost.
Indeed, in the \texttt{MICE} hierarchy constructed along the stochastic optimization path $\{\bs{\xi}_\ell\}_{\ell=0}^k$, our generic optimizer optimally decides whether to drop a particular iteration $\ell$ out of the hierarchy or restart it in order to reduce the total optimization work.
Moreover, it can decide if it is advantageous, from the computational work perspective, to clip the hierarchy at some point $\ell$, discarding the iterations before $\ell$.
Since we control the gradients' statistical error using an estimate of the gradient norm, we employ a resampling technique to get a gradient norm estimate, reducing the effect of sampling error and resulting in robust optimizers.
We note in passing that \texttt{MICE} can be adjusted to the case of finite populations; see Remark \ref{rem:mice.optimal.setting.discrete}, for optimization problems arising in supervised machine learning.

Generally speaking, in first-order stochastic optimization algorithms that produce convergent iterates, the mean gradient converges to zero as the number of iterations, $k$, goes to infinity, that is $\norm{\bb{E} [\nabla_{\bs{\xi}}f(\bs{\xi}_k,\bs{\theta})|\bs{\xi}_k]}{}\to0$; however, the gradient covariance,
$\bb{C} [\nabla_{\bs{\xi}}f(\bs{\xi}_k,\bs{\theta}),\nabla_{\bs{\xi}}f(\bs{\xi}_k,\bs{\theta})|\bs{\xi}_k]$,
does not.
Thus, to ensure convergence of the iterates $\bs{\xi}_k$, in the literature it is customary to use decreasing step-size (learning rate) schedules, reducing the effect of the statistical error in the gradient onto the iterates $\bs{\xi}_k$.
However, this approach also results in sublinear convergence rates \cite{Rup88}. Another approach to deal with the gradient's statistical error is to increase the sample sizes (batch sizes) while keeping the step-size fixed, thus avoiding worsening the convergence.
Byrd et al. \cite{Byr12} propose to adaptively increase the sample sizes to guarantee that the trace of the covariance matrix of the mean gradient is proportional to its norm.
This approach forces the statistical error to decrease as fast as the gradient norm.
Balles et al. \cite{balles2016coupling} use a similar approach; however, instead of setting a parameter to control the statistical error, they set a step-size and find the parameter that guarantees the desired convergence.
Bollapragada et al. \cite{bollapragada2018adaptive} propose yet another approach to control the variance of gradient estimates in stochastic optimization, which they call the inner product test.
Their approach ensures that descent directions are generated sufficiently often.

Instead of increasing the sample size, some methods rely on using control variates with respect to previously sampled gradients to reduce the variance in current iterations and thus be able to keep a fixed step-size.
Pioneering ideas of control variates in stochastic optimization, by Johnson \& Zhang \cite{Joh13}, profit on an accurate mean gradient estimation at the initial guess $\bs{\xi}_0$, $\nabla_{\bs{\xi}}\bb{E} [f(\bs{\xi}_0,\bs{\theta})|\bs{\xi}_{0}]$, to update and compute, via single control variates, an inexpensive and accurate version of the mean gradient at the iteration $k$, $\nabla_{\bs{\xi}}\bb{E} [f(\bs{\xi}_k,\bs{\theta})|\bs{\xi}_k]$.
Instead of doing control variates with respect to one starting full-gradient, \texttt{SARAH}, by Nguyen et al. \cite{nguyen2017sarah}, computes an estimate of the gradient at the current iteration by using control variates with respect to the last iteration.
An `inexact' version of \texttt{SARAH} is presented in \cite{nguyen2018inexact}, where \texttt{SARAH} is generalized to the minimization of expectations.
In the spirit of successive control variates, \texttt{SPIDER} by Fang et al. \cite{fang2018spider} uses control variates between subsequent iterations; however, it employs the \emph{normalized} gradient descent instead of plain gradient descent.
In a different approach, \texttt{SAGA}, by Defazio et al. \cite{defazio2014saga}, keeps in the memory the last gradient $\nabla_{\bs{\xi}} f$ observed for each data point and computes $\nabla_{\bs{\xi}}\bb{E} [f(\bs{\xi}_k,\bs{\theta})|\bs{\xi}_k]$ using control variates with respect to the average of this memory.
Lastly, many algorithms try to `adapt' the initial batch size of the hierarchy of batches using predefined rules, such as exponential or polynomial growth, as presented by Friedlander \& Schmidt \cite{friedlander2012hybrid}, or based on statistical bounds as discussed by De Soham et al. \cite{de2016big} and Ji et al. \cite{ji2019faster}, to mention a few.

Although \emph{Multi-Iteration Stochastic Optimizers} share similarities with \texttt{SVRG}  \cite{Joh13}, \texttt{SARAH}, and \texttt{SPIDER}, our stochastic optimizers distinctly control the relative variance in gradient estimates.
We achieve this control by sampling the entire hierarchy of iterations,
optimally distributing the samples to minimize the gradient sampling cost.
Moreover, we provide additional flexibility by including dropping, restart, and clipping operations in the \texttt{MICE} hierarchy updates.

Since \texttt{MICE} is non-intrusive and designed for both continuous and discrete random variables, it can be coupled with most available optimizers with ease. For instance, we couple \texttt{MICE} with \texttt{SGD} \cite{robbins1951stochastic} and \texttt{Adam} \cite{kingma2014adam}, showing the robustness of our approach.
The \texttt{Adam} algorithm by Kingma \& Ba \cite{kingma2014adam} does not exploit control variates techniques for variance reduction.
Instead, it reduces the gradient estimator's variance based on iterate history by adaptive estimates of lower-order moments, behaving similarly to a filter.
Thus, the coupling \texttt{Adam-MICE} profits from the information available in the optimizer path in more than one way.
Finally, the reader is referred to the books by Spall \cite{spall2005introduction} and Shapiro, Dentcheva, and Ruszczy{\'n}ski \cite{shapiro2014lectures} for comprehensive overviews on stochastic optimization.

To assess \texttt{MICE}'s applicability, we numerically minimize expectations of continuous and discrete random variables using analytical functions and logistic regression datasets, respectively.
Also, we compare \texttt{SGD-MICE} with \texttt{SVRG}, \texttt{SARAH}, \texttt{SAG}, and \texttt{SAGA} in training the logistic regression model with datasets with different sizes and numbers of features.

The remainder of this work is as follows.
In \S\ref{sc:mice}, we construct the \texttt{MICE} statistical estimator \S\ref{sc:mice.construction}; define its statistical error \S\ref{sc:mice.statistical.error}; compute the optimal  number of samples for the current hierarchy  \S\ref{sc:mice.optimal.setting}; present and derive the criteria for dropping iterations out of the hierarchy \S\ref{sc:mice.dropping}, restarting the hierarchy \S\ref{sc:mice.restart}, and clipping the hierarchy \S\ref{sc:mice.clipping}.
In  \S\ref{sc:mice.conv_and_work}, we present a simplified study of the convergence \S\ref{sc:mice.convergence} and complexity \S\ref{sc:mice.work} of \emph{Multi-Iteration Stochastic Optimizers} for strongly-convex and $L$-smooth functions.
In \S\ref{sc:mice_algorithm}, a pseudocode of the \texttt{MICE} algorithm is presented.
In \S\ref{sc:numerics}, to assess the efficiency of \emph{Multi-Iteration Stochastic Optimizers}, we present some numerical examples, ranging from analytical functions to the training of a logistic regression model over a dataset with millions of data points.
In Appendix \ref{sc:ideal_mice}, we present implementation details of the \emph{idealized} \texttt{MICE} with stochastic gradient descent that is used to validate the theory presented in \S\ref{sc:mice.conv_and_work}.
In Appendix \ref{sc:algorithms}, detailed pseudocodes for the \emph{Multi-Iteration Stochastic Optimizers} used in this work are presented.

\subsection{Optimization of expectations and stochastic optimizers}

To state the stochastic optimization problem, let $\bs{\xi}$ be the design variable in dimension $d_{\bs{\xi}}$ and $\bs{\theta}$ a vector-valued random variable in dimension $d_{\theta}$, whose probability distribution $\pi$ may depend on  $\bs{\xi}.$ Throughout this work we assume that we can produce as many independent, identically distributed samples from $\pi$ as needed.
Here, $\bb{E}[\cdot|\bs{\xi}]$, $\bb{V}[\cdot|\bs{\xi}]$, $\bb{C}[\cdot,\cdot|\bs{\xi}]$ are respectively the expectation, variance, and covariance operators conditioned on $\bs{\xi}$. Aiming at optimizing conditional expectations on $\bs{\xi}$, we state our problem as follows. Find $\bs{\xi}^*$ such that
\begin{equation}\label{eq:underlying.problem}
\bs{\xi}^* = \underset{\bs{\xi} \in \bb{R}^{d_{\bs{\xi}}}}{\arg\min} \,  \bb{E} [f(\bs{\xi},\bs{\theta})|\bs{\xi}],
\end{equation}
where
$f \colon \bb{R}^{d_{\bs{\xi}}} \times \bb{R}^{d_{\theta}} \to \bb{R}$. Through what follows, let the objective function in our problem be denoted by $F(\bs{\xi})\coloneqq\bb{E} [f(\bs{\xi},\bs{\theta})|\bs{\xi}]$.

In minimizing \eqref{eq:underlying.problem} with respect to the design variable $\bs{\xi} \in \bb{R}^{d_{\bs{\xi}}}$, \texttt{SGD} is constructed with the following updating rule
\begin{equation}\label{eq:sgd}
\bs{\xi}_{k+1} = \bs{\xi}_{k} - \eta_k \bs{\upsilon}_k,
\end{equation}
where $\eta_k>0$ is the step-size and $\bs{\upsilon}_k$ is an unbiased estimator of the gradient of $F$ at $\bs{\xi}_k$.
For instance, an unbiased estimator $\bs{\upsilon}_k$ of the gradient of $F$ at $\bs{\xi}$ at the iteration $k$ may be constructed by means of a Monte Carlo estimator, namely
\begin{equation}\label{eq:grad.def.2}
\nabla_{\bs{\xi}} F(\bs{\xi}_k) =  \bb{E}[\nabla_{\bs{\xi}} f(\bs{\xi}_k,\bs{\theta})|\bs{\xi}_k] \approx \bs{\upsilon}_k :=\dfrac{1}{M} \sum_{\alpha\in\cl{I}} \nabla_{\bs{\xi}} f(\bs{\xi}_k,\bs{\theta}_{\alpha}),
\end{equation}
with $M$ independent and identically distributed (iid) random variables $\bs{\theta}_\alpha\sim\pi$ given $\bs{\xi}_k$,  $\alpha\in\cl{I}$, with $\cl{I}$ being an index set with cardinality $M\coloneqq|\cl{I}|$. Bear in mind that an estimator of the type \eqref{eq:grad.def.2} is, in fact, a random variable and its use in optimization algorithms gives rise to the so called \emph{Stochastic Optimizers}. The challenge of computing the gradient of $F$ in an affordable and accurate manner motivated the design of several gradient estimators.

For the sake of brevity, the following  review on control variates techniques for stochastic optimization is not comprehensive. To motivate our approach, we recall the control variates proposed by Johnson \& Zhang \cite{Joh13} (and similarly, by Defazio et al.\ \cite{defazio2014saga}) for the optimization of a function defined by a finite sum of functions. The idea of control variates is to add and subtract the same quantity, that is,
\begin{equation}\label{eq:grad.svrg}
\nabla_{\bs{\xi}} F (\bs{\xi}_k) = \bb{E} \left[\nabla_{\bs{\xi}} f(\bs{\xi}_k,\bs{\theta}) - \nabla_{\bs{\xi}} f(\bs{\xi}_0,\bs{\theta})|\bs{\xi}_k,\bs{\xi}_0 \right] + \bb{E} \left[ \nabla_{\bs{\xi}} f(\bs{\xi}_0,\bs{\theta})|\bs{\xi}_0 \right],
\end{equation}
rendering the following sample-based version
\begin{equation}\label{eq:grad.svrg.discrete}
\nabla_{\bs{\xi}} F (\bs{\xi}_k) \approx \dfrac{1}{M_k} \sum_{\alpha\in\cl{I}_k} (\nabla_{\bs{\xi}} f(\bs{\xi}_k,\bs{\theta}_{\alpha}) - \nabla_{\bs{\xi}} f(\bs{\xi}_0,\bs{\theta}_{\alpha})) + \dfrac{1}{M_0-M_k}\sum_{\alpha\in\cl{I}_0 \!\setminus \cl{I}_k} \nabla_{\bs{\xi}} f(\bs{\xi}_0,\bs{\theta}_{\alpha}),
\end{equation}
where $M_0 \gg M_k$ and $\bs{\theta}_{\alpha}$ are iid samples from the $\pi$ distribution, which does not depend on $\bs{\xi}$ in their setting.
In the original work by Johnson \& Zhang \cite{Joh13}, $M_0$ is the total population and $M_k=1$. Later, Nitanda \cite{Nit14} and Kone\v{c}n\`y et al.\ \cite{Kon16} also used the total populations $M_0$ at $\bs{\xi}_0$, but with $M_k=2,4,8\ldots$, to study the efficiency of the algorithm. Additionally, the work \cite{Joh13} restarts the algorithm after a pre-established number of iterations by setting $\bs{\xi}_0\leftarrow\bs{\xi}_k$. The efficiency of this algorithm relies on the correlation between the components of the gradients $\nabla_{\bs{\xi}} F (\bs{\xi}_0)$ and $\nabla_{\bs{\xi}} F (\bs{\xi}_k)$. If this correlation is high, the variance of the mean gradient estimator \eqref{eq:grad.svrg.discrete} is reduced.

\section{Multi-iteration stochastic optimizers}\label{sc:mice}

\subsection{Multi-iteration gradient estimator}\label{sc:mice.construction}

We now construct an affordable estimator of the mean gradient at the current iteration $k$,
$\nabla_{\bs{\xi}} F(\bs{\xi}_k) = \bb{E}[\nabla_{\bs{\xi}}f(\bs{\xi}_k,\bs{\theta})|\bs{\xi}_k]$. We name the estimator \eqref{eq:mice}
\emph{Multi-Iteration stochastiC Estimator}---\texttt{MICE}, see Definition \ref{df:general.mice.estimator} below.
Profiting from available information already computed in previous iterations, \texttt{MICE} uses multiple control variates between pairs of, possibly non-consecutive, iterations along the optimization path to approximate the mean gradient at the iteration $k$.
Bearing in mind that stochastic optimization algorithms, in a broad sense, create a strongly convergent path where $\norm{\bs{\xi}_k-\bs{\xi}_\ell}{} \to 0$ as $\ell,k \to \infty$,
the gradients evaluated at $\bs{\xi}_{\ell}$ and $\bs{\xi}_k$ should become more and more correlated for $k,\ell \to \infty$. In this scenario, control variates with respect to previous iterations become more efficient, in the sense that one needs fewer and fewer new samples to accurately estimate the mean gradient.

To introduce the \texttt{MICE} gradient estimator below in Definition \ref{df:general.mice.estimator}, we need first some notation. Let $\cl{L}_k$ be an index set, such that, $\cl{L}_k\subset\{0,\ldots,k\}$, where $k$ is the current iteration and $k\in\cl{L}_k$.
{This index set is just  $\cl{L}_0= \{0\}$ at the initial iteration, $k=0$, and for later iterations it contains the indices of the iterations \texttt{MICE} uses to reduce the computational work at the current iteration, $k>0$, via control variates.}

Next, for any $\min\{\cl{L}_k\}<\ell\in \cl{L}_k$, let $p_k(\ell)$ be the \emph{previous} from $\ell$ in $\cl{L}_k$  given by
\begin{equation}\label{eq:previous.operator}
p_k(\ell)\coloneqq\max\{\ell^\prime\in\cl{L}_k\colon\ell^\prime<\ell\}.
\end{equation}
Then, the mean gradient at $\bs{\xi}_k$ conditioned on the sequence of random iterates, $\bs{\xi}$, indexed by the set $\cl{L}_k$ can be decomposed as
\begin{equation}\label{eq:mice.continuous}
\bb{E}[\nabla_{\bs{\xi}} f(\bs{\xi}_k, \bs{\theta})|\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}]=\sum_{\ell\in\cl{L}_k}\bs{\mu}_{\ell,k},\qquad\bs{\mu}_{\ell,k}\coloneqq\bb{E}[\Delta_{\ell,k}|\bs{\xi}_\ell,\bs{\xi}_{p_k(\ell)}],
\end{equation}
with the gradient difference notation
\begin{equation}\label{eq:Delta.continuous}
\Delta_{\ell,k} \coloneqq
\left\{
\aligned
\nabla_{\bs{\xi}} f(\bs{\xi}_\ell, \bs{\theta}) - \nabla_{\bs{\xi}} f(\bs{\xi}_{p_k(\ell)}, \bs{\theta}),& \text{ if } \ell > \min\{\cl{L}_k\}, \\
\nabla_{\bs{\xi}} f(\bs{\xi}_\ell, \bs{\theta}),& \text{ if  } \ell = \min\{\cl{L}_k\}.
\endaligned
\right.
\end{equation}
so the conditional mean $\bs{\mu}_{\ell,k}$ defined in \eqref{eq:mice.continuous} is simply
\begin{equation}
\bs{\mu}_{\ell,k} =
\left\{
\aligned
\bb{E}[\nabla_{\bs{\xi}} f(\bs{\xi}_\ell, \bs{\theta})|\bs{\xi}_\ell] -
\bb{E}[\nabla_{\bs{\xi}} f(\bs{\xi}_{p_k(\ell)}, \bs{\theta})|\bs{\xi}_{p_k(\ell)}] ,& \text{ if } \ell > \min\{\cl{L}_k\}, \\
\bb{E}[\nabla_{\bs{\xi}} f(\bs{\xi}_\ell, \bs{\theta})|\bs{\xi}_\ell]  ,& \text{ if } \ell =\min\{\cl{L}_k\}.
\endaligned
\right.
\end{equation}

In what follows, for readability sake,  we make the assumption that the distribution of $\bs{\theta}$ does not depend on $\bs{\xi}.$
Observe that this assumption is more general than it may seem, see the discussion on Remark \ref{rem:theta_distribution}.
\begin{ass}[Simplified probability distribution of $\bs{\theta}$]\label{as:theta_distribution}
The probability distribution of $\bs{\theta}$, $\pi$, does not depend on $\bs{\xi}$.
\end{ass}
Now we are ready to introduce the \texttt{MICE} gradient estimator.
\begin{set}[\texttt{MICE} gradient estimator]\label{df:general.mice.estimator}
Given an index set $\cl{L}_k$ such that $k\in\cl{L}_k\subset\{0,\ldots,k\}$ and nonzero integer numbers $\{M_{\ell,k}\}_{\ell\in\cl{L}_k}$, we define the \texttt{MICE} gradient estimator for $\nabla_{\bs{\xi}}F(\bs{\xi}_k)$ at iteration $k$ as
\begin{equation}\label{eq:mice}
\nabla_{\bs{\xi}}\cl{F}_k\coloneqq\sum_{\ell\in\cl{L}_k}\dfrac{1}{M_{\ell,k}}\sum_{\alpha\in\cl{I}_{\ell,k}}\Delta_{\ell,k,\alpha},
\end{equation}
denoting as before the difference to the previous gradient as
\begin{equation}\label{eq:Delta.sampled}
\Delta_{\ell,k,\alpha} \coloneqq \nabla_{\bs{\xi}} f(\bs{\xi}_\ell, \bs{\theta}_\alpha) - \nabla_{\bs{\xi}} f(\bs{\xi}_{p_k(\ell)}, \bs{\theta}_\alpha).
\end{equation}
Here, for each index $\ell\in\cl{L}_k$, the set of samples, $\cl{I}_{\ell,k}$, has cardinality $M_{\ell,k}$.
\end{set}
Definition \ref{df:general.mice.estimator} accounts for the case when dropping iterations from the \texttt{MICE} hierarchy (see \S\ref{sc:mice.dropping}) or just restarting the hierarchy (see \S\ref{sc:mice.restart}). For instance, if iteration $2$ is to be dropped, $\Delta_{\ell,k,\alpha}$ in \eqref{eq:Delta.sampled} at iteration $\ell=3$ is computed with the pair $(p_k(\ell=3)=1,\ell=3)$. Also, if a restart is to be performed at the current iteration $k$, then we take $\cl{L}_k=\{k\}$ and the previous hierarchy is deleted.

\begin{rmk}[Cumulative sampling in \texttt{MICE}]
As the stochastic optimization progresses, new additional samples of $\bs{\theta}$ are taken and others, already available from previous iterations, are reused to compute the  \texttt{MICE} estimator at the current iteration.
This sampling procedure is defined by the couples $(M_{\ell,k},\bs{\xi}_{\ell})_{\ell\in\cl{L}_k}$, making $\bs{\xi}_{k+1}$ a deterministic function of all the samples in the hierarchy $\cl{L}_k$.
 For example, using \texttt{SGD-MICE} with step length $\eta>0$ on the first iteration, we have
\begin{equation}\label{eq:sgdb0}
\bs{\xi}_{1} = \bs{\xi}_{0} - \eta \nabla_{\bs{\xi}}\cl{F}_0,
\end{equation}
with
\begin{equation}
\nabla_{\bs{\xi}}\cl{F}_0 = \dfrac{1}{M_{0,0}} \sum_{\alpha\in\cl{I}_{0,0}} \nabla_{\bs{\xi}} f(\bs{\xi}_0,\bs{\theta}_{\alpha}).
\end{equation}
Conversely, we have that the conditional expectation of the $\texttt{MICE}$ estimator given the values of $\bs{\xi}_0$ and $\bs{\xi}_1$ then becomes
\begin{equation}
\bb{E}[ \nabla_{\bs{\xi}}\cl{F}_0 | \bs{\xi}_0,\bs{\xi}_1 ] =
\frac{\bs{\xi}_{0} -\bs{\xi}_{1}}{\eta}.
\end{equation}
However, the $(M_{0,1}-M_{0,0})$ new samples of  $\nabla_{\bs{\xi}} f(\bs{\xi}_0,\bs{\theta})$, taken at iteration $k=1$ once $\bs{\xi}_1$ is known,
will be independent  from $\bs{\xi}_1$.
Therefore, we can see that the samples $\Delta_{\ell,k,\alpha}$ with $\alpha\in\cl{I}_{\ell,k}$ used in \eqref{eq:mice} are in general heterogeneous, since they may be sampled based on different information on the $\bs{\xi}$ iterations, and thus only conditionally independent.
Some of those realizations, corresponding to the underlying random input $\bs{\theta}$,
 are sampled conditioned on the knowledge of the pair $(\bs{\xi}_{p_k(\ell)},\bs{\xi}_\ell)$
and possibly other future iterates $\bs{\xi}_{\ell^\prime},$ with $\ell<\ell^\prime\in\cl{L}_k$.
\end{rmk}

Note that the \texttt{MICE} gradient estimator introduced in Definition \eqref{df:general.mice.estimator} is conditionally unbiased
if all the $\Delta_{\ell,k,\alpha}$ are conditionally sampled based on the knowledge of $\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}$, that is,
\begin{equation}
\bb{E}[ \nabla_{\bs{\xi}}\cl{F}_k |\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}]= \bb{E}[\nabla_{\bs{\xi}} f(\bs{\xi}_k, \bs{\theta})|\bs{\xi}_k]= \nabla_{\bs{\xi}} F(\bs{\xi}_k) .
\end{equation}
This happens for instance at $k=0$ and at any other iteration where the hierarchy is restarted.
As we shall see later in Remark \ref{rem:mice.bias}, the unbiased property is indeed more subtle once we start reusing samples from previous iterations {and may not be fulfilled in that general case}.

As for the construction of the \texttt{MICE} hierarchy at iteration $k$, that is, $\cl{L}_k$, from the previous one, $\cl{L}_{k-1}$, consider the following Definition.
Once the new iterate $\bs{\xi}_k$ is known, the new hierarchy $\cl{L}_k$ is  constructed as follows:
\begin{set}\label{df:add.iteration.consequences}[Construction of the hierarchy $\cl{L}_k$]
For $k=0$, let $\cl{L}_0 = \{0\}$. If $k\ge 1$ then take
\begin{equation}\label{eq:Lk_1ststep}
\cl{L}_k\leftarrow\cl{L}_{k-1}\cup\{k\}.
\end{equation}
After this step, there are three possible cases to finish the construction of $\cl{L}_k$:
\begin{description}
\item[Regular case] \texttt{MICE} keeps the hierarchy $\cl{L}_k$ in \eqref{eq:Lk_1ststep}.
\item[Drop case] \texttt{MICE} drops the last iteration, that is, $k-1$. Thus, the hierarchy is reduced as follows: $\cl{L}_k\leftarrow\cl{L}_k\!\setminus\{k-1\}.$
\item[Restart case] \texttt{MICE} restarts the hierarchy, keeping only iteration $k$, that is, $\cl{L}_k\leftarrow\{k\}$.
\end{description}
\end{set}
For more details, see \S\ref{sc:mice_algorithm} for an algorithmic description.
As for the additional cost in iterating from $\bs{\xi}_{k-1}$ to $\bs{\xi}_k$, consider the following Remark.
\begin{rmk}[\texttt{MICE} estimator  sampling cost  at iteration $k$]\label{rem:add.iteration.cost}
Here, we discuss the three scenarios introduced in Definition \ref{df:add.iteration.consequences}, detailing the associated cost in computing the \texttt{MICE} estimator at iteration $k$.
Denote, with $\Delta_{\ell,k,\alpha}$ as in \eqref{eq:Delta.sampled}, the sample average
\begin{equation}\label{eq:muhat_average_def}
\hat{\bs{\mu}}_{\ell,k} \coloneqq \frac{1}{M_{\ell,k}}\sum_{\alpha\in\cl{I}_{\ell,k}}\Delta_{\ell,k,\alpha}.
\end{equation}
Then, we can rewrite \texttt{MICE} in \eqref{eq:mice} as follows, distinguishing the terms that can be computed using previous available information and those that need fresh gradient evaluations at iteration $k$:
\begin{description}
\item[Regular case]
Rephrasing expression \eqref{eq:mice} to make explicit the difference between the sunk cost and the additional cost, the \texttt{MICE} estimator with the hierarchy $\cl{L}_k$ reads
\begin{equation}\label{eq:regular}
\nabla_{\bs{\xi}}\cl{F}_k=\underbrace{\Bigg\{\sum_{\ell\in\cl{L}_{k-1}}\dfrac{M_{\ell,k-1}}{M_{\ell,k}}\hat{\bs{\mu}}_{\ell,k-1}\Bigg\}}_{\text{Sunk cost, previous iterations}}+\underbrace{\Bigg\{\hat{\bs{\mu}}_{k,k}+\sum_{\ell\in\cl{L}_{k-1}}\dfrac{1}{M_{\ell,k}}\sum_{\alpha\in\cl{I}_{\ell,k}\!\setminus\cl{I}_{\ell, k-1}}\Delta_{\ell,k,\alpha}\Bigg\}}_{\text{Additional \texttt{MICE} cost incurred at iteration $k$}},
\end{equation}
Observe that here only the second and third terms contribute to the increment in the computational cost, namely the cost associated to sampling new gradients.
Thus, the additional cost is proportional to
$|\cl{I}_{k,k}|+\sum_{\ell\in\cl{L}_{k-1}}|\cl{I}_{\ell,k}\!\setminus\cl{I}_{\ell ,k-1}|$;
\item[Dropping case]
Dropping iteration $k-1$ of the hierarchy $\cl{L}_k$, expression \eqref{eq:regular} reads,
\begin{equation}\label{eq:drop}
\nabla_{\bs{\xi}}\cl{F}_k=\underbrace{\Bigg\{\sum_{\ell\in\cl{L}_{k-1}\!\setminus\{k-1\}}\dfrac{M_{\ell,k-1}}{M_{\ell,k}}\hat{\bs{\mu}}_{\ell,k-1}\Bigg\}}_{\text{Sunk cost, previous iterations}}+\underbrace{\Bigg\{\hat{\bs{\mu}}_{k,k}+\sum_{\ell\in\cl{L}_{k-1}\!\setminus\{k-1\}}\dfrac{1}{M_{\ell,k}}\sum_{\alpha\in\cl{I}_{\ell,k}\!\setminus\cl{I}_{\ell ,k-1}}\Delta_{\ell,k,\alpha}\Bigg\}}_{\text{Additional \texttt{MICE} cost incurred at iteration $k$}}.
\end{equation}
Again, only the second and third terms contribute to the increment in the computational cost. Thus, the additional gradient sampling cost is proportional to
$|\cl{I}_{k,k}|+\sum_{\ell\in\cl{L}_{k-1}\!\setminus\{k-1\}}|\cl{I}_{\ell,k}\!\setminus\cl{I}_{\ell,k-1}|$;
\item[Restart case]
In this case, we simply compute the estimator from scratch, namely we base our mean gradient estimator solely on new gradient samples,
\begin{equation}\label{eq:restart}
\nabla_{\bs{\xi}}\cl{F}_k=\hat{\bs{\mu}}_{k,k}.
\end{equation}
Here, the additional cost is the cost of starting from scratch, namely proportional to $|\cl{I}_{k,k}|$.
\end{description}
\end{rmk}
Note that because of the structure of $\cl{L}_k$, when computing $\hat{\bs{\mu}}_{k,k}$, the differences $\Delta_{k,k}$ appearing in \eqref{eq:regular}, \eqref{eq:drop} and \eqref{eq:restart} are not necessarily the same. Indeed, in the regular case we have $\Delta_{k,k} = \nabla_{\bs{\xi}} f(\bs{\xi}_k, \bs{\theta}) - \nabla_{\bs{\xi}} f(\bs{\xi}_{k-1}, \bs{\theta})$
while in the dropping case we have $\Delta_{k,k} =  \nabla_{\bs{\xi}} f(\bs{\xi}_k, \bs{\theta}) - \nabla_{\bs{\xi}} f(\bs{\xi}_{p_{k-1}(k-1)}, \bs{\theta})$. Finally, in the restart case we simply have $\Delta_{k,k} = \nabla_{\bs{\xi}} f(\bs{\xi}_k, \bs{\theta})$. However, the summation over the hierarchy $\cl{L}_k$ and the use of the previous operator \eqref{eq:previous.operator} make the definition of the \texttt{MICE} estimator \eqref{eq:mice} univocal.
We refer the reader to Remark \ref{rem:MICE_cost_work}, which discusses the total cost,  including the aggregation of already known gradients to compute the \texttt{MICE} estimator optimally.
Finally, note that the \texttt{MICE} update, dropping, and restart are performed in an optimal fashion for a given statistical tolerance; see, for instance, \S\ref{sc:mice.optimal.setting}.

\begin{rmk}[About \texttt{MICE} and \texttt{MLMC}]
Note that \texttt{MICE} resembles the estimator obtained in the \emph{Multilevel Monte Carlo} method---\texttt{MLMC} \cite{Heinrich:MLMC,giles:MLMC,Gil15}.
For instance, in the absence of Dropping, Clipping and Restart, \texttt{MICE} reads
\begin{equation}
\nabla_{\bs{\xi}} \cl{F}_k \eqqcolon
  \dfrac{1}{M_{0,k}}\sum_{\alpha\in\cl{I}_{0,k}} \nabla_{\bs{\xi}} f(\bs{\xi}_0,\bs{\theta}_{\alpha})  + \sum _{\ell=1}^k \dfrac{1}{M_{\ell,k}} \sum_{\alpha\in\cl{I}_{\ell,k}} (\nabla_{\bs{\xi}} f(\bs{\xi}_\ell,\bs{\theta}_{\alpha}) - \nabla_{\bs{\xi}} f(\bs{\xi}_{\ell-1},\bs{\theta}_{\alpha})).
\end{equation}
Indeed, we may think that in \texttt{MICE}, the iterations play the same role as the levels of approximation in  \texttt{MLMC}.
However, there are several major differences with  \texttt{MLMC}, namely \emph{i}) \texttt{MICE} exploits sunk cost of previous computations, computing afresh only what is necessary to have enough accuracy on the current iteration \emph{ii}) there is dependence in \texttt{MICE} across iterations and \emph{iii}) in \texttt{MICE}, the sample cost for the gradients is the same in different iterations while in  \texttt{MLMC} one usually has higher cost per sample for deeper, more accurate levels.

Indeed, assuming the availability of a convergent hierarchy of approximations and following the \texttt{MLMC} lines, the work \cite{multilevel_RobMon} proposed and analyzed multilevel stochastic approximation algorithms, essentially recovering the classical error bounds for multilevel Monte Carlo approximations in this more complex context.
In a similar  \texttt{MLMC} hierarchical approximation framework, the work by Yang, Wand, and Fang \cite{yang2019multilevel} proposed a stochastic gradient algorithm for solving optimization problems with nested expectations as objective functions.
Last, the combination of \texttt{MICE} and the \texttt{MLMC} ideas like those in  \cite{multilevel_RobMon} and \cite{yang2019multilevel}  is thus a natural research avenue to pursue.
\end{rmk}

Let us consider now the conditional bias of the \texttt{MICE} gradient estimator introduced in Definition \eqref{df:general.mice.estimator} when we
reuse samples from previous iterations.
\begin{rmk}[Conditional bias in \texttt{MICE}]\label{rem:mice.bias}
To introduce the bias issue, let us consider a regular case given by \eqref{eq:regular}, with only two iterates, $(\bs{\xi}_0,\bs{\xi}_1)$:
\begin{equation}\label{eq:regular_rem}
\nabla_{\bs{\xi}}\cl{F}_1=\underbrace{\dfrac{M_{0,0}}{M_{0,1}}\, \nabla_{\bs{\xi}}\cl{F}_0}_{\text{Samples computed at $k=0$}}
+\underbrace{\Bigg\{\hat{\bs{\mu}}_{1,1}
+\dfrac{1}{M_{0,1}}\sum_{\alpha\in\cl{I}_{0,1}\!\setminus\cl{I}_{0, 0}}
\nabla_{\bs{\xi}} f(\bs{\xi}_0,\theta_\alpha)
\Bigg\}}_{\text{Samples computed at $k=1$}},
\end{equation}
with $\hat{\bs{\mu}}_{1,1}$ as in \eqref{eq:muhat_average_def}.
Let us assume that we iterate from $\bs{\xi}_0$ to $\bs{\xi}_1$ using \texttt{SGD-MICE} with step length $\eta>0$, namely that
\begin{equation}\label{eq:sgdb}
\bs{\xi}_{1} = \bs{\xi}_{0} - \eta \nabla_{\bs{\xi}}\cl{F}_0.
\end{equation}
Then, from \eqref{eq:regular_rem} and \eqref{eq:sgdb} we can see that
\begin{equation}
\bb{E}[ \nabla_{\bs{\xi}}\cl{F}_1 |\bs{\xi}_0,\bs{\xi}_1]=
 \dfrac{M_{0,0}}{M_{0,1}}\, \frac{\bs{\xi}_{0} -\bs{\xi}_{1} }{\eta}
 + \nabla_{\bs{\xi}}{F}(\bs{\xi}_1) - \nabla_{\bs{\xi}}{F}(\bs{\xi}_0)
 +  \dfrac{M_{0,1}-M_{0,0}}{M_{0,1}}\,\nabla_{\bs{\xi}}{F}(\bs{\xi}_0),
\end{equation}
yielding the conditional bias
\begin{equation}\label{eq:BIAS_ex}
\bb{E}[ \nabla_{\bs{\xi}}\cl{F}_1 |\bs{\xi}_0,\bs{\xi}_1] - \nabla_{\bs{\xi}}{F}(\bs{\xi}_1) =
 \dfrac{M_{0,0}}{M_{0,1}}\, \left(
 \frac{\bs{\xi}_{0} -\bs{\xi}_{1} }{\eta} -  \nabla_{\bs{\xi}}{F}(\bs{\xi}_0)
 \right)
\end{equation}
which is not zero in general.
\end{rmk}
The conditional bias in \eqref{eq:BIAS_ex} is subtle: observe that we have
\begin{equation}
\bb{E}\left[ \frac{\bs{\xi}_{0} -\bs{\xi}_{1} }{\eta}
-\nabla_{\bs{\xi}}{F}(\bs{\xi}_0) \bigg| \bs{\xi}_{0}\right] = 0,
\end{equation}
and that
\begin{equation}
\bb{E}\left[ \left\| \frac{\bs{\xi}_{0} -\bs{\xi}_{1} }{\eta}
- \nabla_{\bs{\xi}}{F}(\bs{\xi}_0)\right\|^2
\bigg| \bs{\xi}_{0}\right]
= \dfrac{\sum_{i=1}^{d_{\bs{\xi}}} \,
\bb{V}[\partial_i f(\bs{\xi}_0)|\bs{\xi}_0]}{\eta^2 M_{0,0}}
= \cl{O}(1/M_{0,0}),
\end{equation}
where $\partial_i$ represents the conventional partial derivative with respect to the $i$-th component of $\bs{\xi}$. Fortunately, the effect of the conditional bias on the convergence of the \texttt{MICE}-based stochastic optimizers does not seem to be relevant in practice, as shown in the numerical examples from \S\ref{sc:numerics}. In \S\ref{sc:numerics}, we compare  different versions of \texttt{MICE} with corresponding unbiased ones, observing no relevant deterioration in their convergence.

\subsection{\texttt{MICE} gradient statistical error}\label{sc:mice.statistical.error}

To determine the optimal number of samples per iteration $\ell\in\cl{L}_k$, we begin by defining the square of the statistical error, $\cl{E}_{\text{stat}}$, as the squared $L^2$-distance between \texttt{MICE} estimator \eqref{eq:mice} and its conditional expectation, which leads to
\begin{equation}\label{eq:mse.def}
\left(\cl{E}_{\text{stat}}(k)\right)^2
\coloneqq
\bb{E} \left[  \left\|
\nabla_{\bs{\xi}}\cl{F}_k - \bb{E} \left[\nabla_{\bs{\xi}}\cl{F}_k\vert\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}\right]
\right\|^2
\vert \{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}
        \right].
\end{equation}
\begin{rmk}[Conditional statistical error in \texttt{MICE}]\label{rem:mice.stat_error}
To introduce the statistical error as we did with the bias in Remark \ref{rem:mice.bias}, let us consider a regular case given by \eqref{eq:regular}, with only two iterates, $(\bs{\xi}_0,\bs{\xi}_1)$:
\begin{equation}\label{eq:regular_rem1}
\nabla_{\bs{\xi}}\cl{F}_1=\underbrace{\dfrac{M_{0,0}}{M_{0,1}}\, \nabla_{\bs{\xi}}\cl{F}_0}_{\text{Samples computed at $k=0$}}
+\underbrace{\Bigg\{\hat{\bs{\mu}}_{1,1}
+\dfrac{1}{M_{0,1}}\sum_{\alpha\in\cl{I}_{0,1}\!\setminus\cl{I}_{0, 0}}
\nabla_{\bs{\xi}} f(\bs{\xi}_0,\theta_\alpha)
\Bigg\}}_{\text{Samples computed at $k=1$}}
\end{equation}
with $\hat{\bs{\mu}}_{1,1}$ as in \eqref{eq:muhat_average_def}.
Taking the expected value of expression \eqref{eq:regular_rem1} conditioned on the pair $(\bs{\xi}_0,\bs{\xi}_1)$, we have that
\begin{equation}\label{eq:expectation.regular_rem1}
\bb{E}[ \nabla_{\bs{\xi}}\cl{F}_1 |\bs{\xi}_0,\bs{\xi}_1]=  \dfrac{M_{0,0}}{M_{0,1}}\, \frac{\bs{\xi}_{0} -\bs{\xi}_{1} }{\eta}
 + \nabla_{\bs{\xi}}{F}(\bs{\xi}_1) - \nabla_{\bs{\xi}}{F}(\bs{\xi}_0)
 +  \dfrac{M_{0,1}-M_{0,0}}{M_{0,1}}\,\nabla_{\bs{\xi}}{F}(\bs{\xi}_0).
\end{equation}
Next, subtracting expression \eqref{eq:regular_rem1} from \eqref{eq:expectation.regular_rem1}, we arrive at
\begin{multline}
\nabla_{\bs{\xi}}\cl{F}_1- \bb{E}[ \nabla_{\bs{\xi}}\cl{F}_1 |\bs{\xi}_0,\bs{\xi}_1] = \dfrac{M_{0,0}}{M_{0,1}}\, \left(\nabla_{\bs{\xi}}\cl{F}_0 - \frac{\bs{\xi}_{0} -\bs{\xi}_{1} }{\eta}\right) \\
+ \left(\hat{\bs{\mu}}_{1,1} - \nabla_{\bs{\xi}}{F}(\bs{\xi}_1) +\nabla_{\bs{\xi}}{F}(\bs{\xi}_0) \right) + \dfrac{1}{M_{0,1}}\sum_{\alpha\in\cl{I}_{0,1} \!\setminus \cl{I}_{0, 0}} (\nabla_{\bs{\xi}} f(\bs{\xi}_0,\theta_\alpha)-\nabla_{\bs{\xi}}{F}(\bs{\xi}_0)),
\end{multline}
and then we can conclude, using conditioning, that the statistical error takes the form
\begin{align}
\left(\cl{E}_{\text{stat}}(1)\right)^2 = {}& \frac{M_{0,1}-M_{0,0}}{(M_{0,1})^2}
\bb{E}[\| \nabla_{\bs{\xi}} f(\bs{\xi}_0,\theta)-\nabla_{\bs{\xi}}{F}(\bs{\xi}_0) \|^2|\bs{\xi}_0], \nonumber\\
& + \frac{1}{M_{1,1}} \bb{E}[\| (\nabla_{\bs{\xi}} f(\bs{\xi}_1,\theta)-\nabla_{\bs{\xi}}{F}(\bs{\xi}_1))
-(\nabla_{\bs{\xi}} f(\bs{\xi}_0,\theta)-\nabla_{\bs{\xi}}{F}(\bs{\xi}_0))\|^2|\bs{\xi}_0,\bs{\xi}_1] \nonumber\\
& \le
\frac{1}{M_{0,1}} \bb{E}[\| \nabla_{\bs{\xi}} f(\bs{\xi}_0,\theta)-\nabla_{\bs{\xi}}{F}(\bs{\xi}_0) \|^2|\bs{\xi}_0] \nonumber\\
& + \frac{1}{M_{1,1}} \bb{E}[\| (\nabla_{\bs{\xi}} f(\bs{\xi}_1,\theta)-\nabla_{\bs{\xi}}{F}(\bs{\xi}_1)) -(\nabla_{\bs{\xi}} f(\bs{\xi}_0,\theta)-\nabla_{\bs{\xi}}{F}(\bs{\xi}_0))\|^2|\bs{\xi}_0,\bs{\xi}_1].
\end{align}
\end{rmk}
The previous remark motivates the following approximation for the statistical error.
Letting $\Delta_{\ell,k,\alpha}^{(i)}$ be the $i$-th component of the $d_{\bs{\xi}}$ dimensional vector $\Delta_{\ell,k,\alpha}$, we define
\begin{equation}\label{eq:Vl}
V_{\ell,k}^{(i)} \coloneqq \bb{V}\left[\Delta_{\ell,k}^{(i)}\bigg\vert\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k} \right], \qquad
 \quad 1 \le i \le d_{\bs{\xi}},
\end{equation}
and
\begin{equation}\label{eq:Velldef}
V_{\ell,k}\coloneqq\sum_{i=1}^{d_{\bs{\xi}}} V_{\ell,k}^{(i)}.
\end{equation}
Then, in this work, we use the following approximation
\begin{equation}\label{eq:variance.estimator}
\left(\cl{E}_{\text{stat}}(k)\right)^2 \approx \sum _{\ell\in\cl{L}_k} \frac{V_{\ell,k}}{M_{\ell,k}}.
\end{equation}
Note that, in practice, the terms $V_{\ell,k}$ are computed using sample approximations for each $\ell\in\cl{L}_k$. In the idealized convergence analysis in \S\ref{sc:mice.convergence}, we assume that they are computed exactly.

\subsection{Multi-iteration optimal setting for gradient variability control} \label{sc:mice.optimal.setting}

First, let the gradient sampling cost and the total \texttt{MICE} work be defined as in the following Remark.
\begin{rmk}[\texttt{MICE} gradient sampling cost and work]\label{rem:MICE_cost_work}
The number of gradient evaluations is $1$ for $\Delta_{\ell, k, \alpha}$ when $\ell=\min{\cl{L}_k}$ and $2$ otherwise.
For this reason, we define the auxiliary index function
\begin{equation}
\mathbbm{1}_{\mskip-2mu\underline{\cl{L}_k}}(\ell)\coloneqq
\begin{cases}
0 & \text{if}\quad\ell=\min{\cl{L}_k},\\
1 & \text{otherwise}
\end{cases},
\end{equation}
and define the gradient sampling cost in number of gradient evaluations as
\begin{equation}
  \cl{C}(\{M_{\ell,k}\}_{\ell\in\cl{L}_k})
  \coloneqq
  \sum_{\ell\in\cl{L}_k} (1+\mathbbm{1}_{\mskip-2mu\underline{\cl{L}_k}}(\ell)) M_{\ell,k}.
\end{equation}
The total work of a \texttt{MICE} evaluation is then the sum of the cost of sampling the gradients and the cost of aggregating the gradients as
\begin{equation}
  \cl{W}\left(
  \{M_{\ell,k}\}_{\ell\in\cl{L}_k}
  \right)
  \coloneqq \cl{C}(
  \{M_{\ell,k}\}_{\ell\in\cl{L}_k}
  ) C_\nabla
  + |\cl{L}_k| C_\text{aggr},
\end{equation}
where $C_\nabla$ is the work of sampling $\nabla_{\bs{\xi}} f$ and $C_\text{aggr}$ is the work of averaging the $\Delta_{\ell, k}$ to construct $\cl{F}_k$.
Then, the work done in iteration $k$ to update \texttt{MICE} is
\begin{equation}\label{eq:work_update_MICE}
  \Delta \cl{W}\left(
  \{M_{\ell,k}\}_{\ell\in\cl{L}_k}
  \right)
  \coloneqq
  \sum _{\ell\in\cl{L}_k} (1+\mathbbm{1}_{\mskip-2mu\underline{\cl{L}_k}}(\ell))
  ( M_{\ell,k} - M_{\ell,k-1})
   C_\nabla
  + |\cl{L}_k| C_\text{aggr}.
\end{equation}
\end{rmk}
Motivated by the analysis of the idealized \texttt{MICE} in \S\ref{sc:mice.conv_and_work}, here we choose the number of samples for the hierarchy $\cl{L}_k$
by approximate minimization of the gradient sampling cost \eqref{eq:work_update_MICE} subject to a given tolerance $\epsilon>0$ on the relative statistical error in the mean gradient approximation, cf. \eqref{eq:variance.estimator}, that is
\begin{equation}\label{eq:M_relstatprob}
\begin{aligned}
\{M^*_{\ell,k}\}_{\ell\in\cl{L}_k}
 = & \underset{\{M_{\ell,k}\}_{\ell\in\cl{L}_k}}{\text{arg min}}
\cl{C}\left(
\{M_{\ell,k}\}_{\ell\in\cl{L}_k}
\right) \\
& \text{subject to} \quad \left(\cl{E}_{\text{stat}}(k)\right)^2 \le \epsilon^2\norm{\nabla_{\bs{\xi}}F(\bs{\xi}_k)}{}^2.
\end{aligned}
\end{equation}

Thus, an approximate integer-valued solution based on Lagrangian relaxation to problem \eqref{eq:M_relstatprob} is
\begin{equation}\label{eq:optimal.sampling.1}
M^*_{\ell,k} = \left\lceil
\frac{1}{\epsilon^2 \norm{\nabla_{\bs{\xi}}F(\bs{\xi}_k)}{}^2}
\left( \sum _{\ell^\prime\in\cl{L}_k} \left(V_{\ell^\prime\!,k}(1+\mathbbm{1}_{\mskip-2mu\underline{\cl{L}_k}}(\ell^\prime))\right)^{1/2} \right) \left(\dfrac{V_{\ell,k}}{(1+\mathbbm{1}_{\mskip-2mu\underline{\cl{L}_k}}(\ell))}\right)^{1/2}\right\rceil,\qquad \forall \ell\in\cl{L}_k.
\end{equation}
In general, in considering the cost of computing new gradients at the iteration $k$, the expenditure already carried out up to the iteration $k-1$ is sunk cost and must not be included, as described in Remark \ref{rem:add.iteration.cost}, that is, one should only consider the incremental cost of going from $k-1$ to $k$. However, since the cost is linear with respect to the samples $M_{\ell,k}$, when considering the increment $\Delta M_{\ell,k}$ in \eqref{eq:M_relstatprob} for $\ell\in\cl{L}_k$, rather than $M_{\ell,k}$, the optimal setting $M^*_{\ell,k}$ remains the same, $M^*_{\ell,k}=M^*_{\ell,k-1} + \Delta M^*_{\ell,k}$.
Moreover, in the variance constraint of problem \eqref{eq:M_relstatprob}, since we do not have access to the  norm of the  mean gradient, $\norm{\nabla_{\bs{\xi}}F(\bs{\xi}_k)}{}$, we use a resampling technique combined with the \texttt{MICE} estimator as an approximation; see Remark \ref{rem:grad_resampling}.

\begin{rmk}[Multi-iteration optimal setting in finite populations]\label{rem:mice.optimal.setting.discrete}
When working with a finite amount of data, the variance of the estimator should account for the ratio between the actual number of samples $M_{\ell,k}$ used in the estimator and the total population $N$, see, for instance, the work by Isserlis \cite{isserlis1918value},
\begin{equation}\label{eq:variance.estimator.finite.pop}
\left(\cl{E}_{\text{stat}}(k)\right)^2  \approx \sum _{\ell\in\cl{L}_k} \frac{V_{\ell,k}}{M_{\ell,k}} \left(\frac{N-M_{\ell,k}}{N-1}\right).
\end{equation}
Similarly as in \eqref{eq:M_relstatprob}, we now use the approximation \eqref{eq:variance.estimator.finite.pop} in the statistical constraint, yielding the following optimal samples,
\begin{equation}\label{eq:optimal.sampling.2}
M^*_{\ell,k} = \Bigg\lceil\frac{N}{N-1}
\frac{
 \left(\sum_{\ell^\prime\in\cl{L}_k}(V_{\ell^\prime\!,k}(1+\mathbbm{1}_{\mskip-2mu\underline{\cl{L}_k}}(\ell^\prime)))^{1/2}\right) \left(\frac{V_{\ell,k}}{(1+\mathbbm{1}_{\mskip-2mu\underline{\cl{L}_k}}(\ell))}\right)^{1/2}
 }{
 \epsilon^2 \norm{\nabla_{\bs{\xi}}F(\bs{\xi}_k)}{}^2
+\frac{1}{N-1}\sum_{\ell^{\prime\prime}\in\cl{L}_k} V_{\ell^{\prime\prime}\!,k}
 }
 \Bigg\rceil,\qquad \forall \ell\in\cl{L}_k.
\end{equation}
\end{rmk}

\begin{rmk}[Gradient resampling for calculating sample sizes]\label{rem:grad_resampling}
To calculate the optimal sample sizes according to \eqref{eq:optimal.sampling.1} and \eqref{eq:optimal.sampling.2}, we need the norm of the mean gradient at $\bs{\xi}_k$, $\|\nabla_{\bs{\xi}} F(\bs{\xi}_k)\|$.
Since the mean gradient is not available, we can use the \texttt{MICE} estimate as an approximation.
However, this estimate might not be accurate enough, leading to sample sizes that are lower than necessary. If that happens then the constraint on \eqref{eq:M_relstatprob} may not be satisfied and the stochastic optimizer may loose efficiency.
Properly estimating $\|\nabla_{\bs{\xi}} F(\bs{\xi}_k)\|$ gets harder as optimization proceeds since $\nabla_{\bs{\xi}} F(\bs{\xi}_k)$ converges to zero while the sample gradients, $\nabla_{\bs{\xi}} f(\bs{\xi}_k,\theta)$, do not in general.
To address this issue, we perform a jackknife \cite{efron1982jackknife} resampling of the approximate mean gradient using sample subsets for each iteration $\ell \in \cl{L}_k$ and use the $p_\text{re}$ percentile of the norms resampled when calculating $M_{\ell,k}$.

First, for each iteration $\ell \in \cl{L}_k$, we partition the index set $\cl{I}_{\ell, k}$ in $n_\text{part}$ disjoint sets $\cl{I}^{(1)}_{\ell, k}, \cl{I}^{(2)}_{\ell, k},.., \cl{I}^{(n_\text{part})}_{\ell, k}$ with the same cardinality.
Then, we create, for each of these sets, their complement with respect to $\cl{I}_{\ell, k}$, i.e., $ \overline{\cl{I}}^{(i)}_{\ell, k} = \cl{I}_{\ell, k} \setminus \cl{I}^{(i)}_{\ell, k}$ for all $i=1,2,..,n_\text{part}$.
We use these complements to construct a set of gradient mean estimates $\{\nabla_{\bs{\xi}}\cl{F}_{k, \nu}\}_{\nu=1}^{n_\text{samp}}$ with
\begin{equation}\label{eq:resampling}
\nabla_{\bs{\xi}}\cl{F}_{k, \nu}
\coloneqq
\sum_{\ell\in\cl{L}_k}
\left(
\left| \overline{\cl{I}}^{\left(i_{\ell, \nu}\right)}_{\ell,k} \right|^{-1} \sum_{\alpha\in\overline{\cl{I}}^{\left(i_{\ell, \nu}\right)}_{\ell,k}}\Delta_{\ell,k,\alpha}
\right),
\end{equation}
by independently sampling $i_{\ell, \nu} \sim \cl{U}\{1,n_\text{part}\}$.
The gradient average norm chosen to calculate the sample-sizes is the largest value below the $p_\text{re}$ percentile, i.e., we order $\nabla_{\bs{\xi}}\cl{F}_{k, \nu}$ increasingly by $\nu$ and choose $\nu^* = \lfloor n_\text{samp} \, p_\text{re} /100 \rfloor$, then,
\begin{equation}
  \|\nabla_{\bs{\xi}} \cl{F}_k^\text{re}\| \coloneqq
  \|\nabla_{\bs{\xi}} \cl{F}_{k, \nu^*}\|.
\end{equation}

To control the work of the resampling procedure, we set its cost as a fraction $\delta_\text{re}$ of work to be done at the current iteration, that is
\begin{equation}
  \cl{W}_\text{re} = \delta_\text{re} C_{\nabla}
  \Delta \cl{W}\left(
  \{M_{\ell,k}^*\}_{\ell\in\cl{L}_k}
  \right)
  ,
\end{equation}
for a gradient cost $C_{\nabla}$.
To obtain the number of samples $n_\text{samp}$, we assume that the resampling cost is $\cl{W}_\text{re} = C_\text{samp} |\cl{L}_k| n_\text{samp}$, where $C_\text{samp}$ is the cost of resampling one iteration from the hierarchy, thus arriving at
\begin{equation}
  n_\text{samp} =
  \max\left\{
  \underline{n_{\text{samp}}},
  \left\lfloor
  \frac{
  \delta_\text{re}
  \Delta \cl{W}\left(
  \{M_{\ell,k}^*\}_{\ell\in\cl{L}_k}
  \right)
  }{
  \frac{C_\text{samp}}{C_{\nabla}} |\cl{L}_k|
  }
  \right\rfloor
  \right\},
\end{equation}
where $\underline{n_{\text{samp}}}$ is the minimum resampling size.
The parameters $\delta_\text{re}$, $n_\text{part}$, and $\underline{n_{\text{samp}}}$ need to be set, whereas the ratio $\nicefrac{C_\text{samp}}{C_\nabla}$ can be measured on-the-fly.
Note that the parameter $\delta_\text{re}$ controls the resampling overhead in terms of processing and the parameter $n_\text{part}$ controls the memory overhead due to storing the partitions.
From our numerical tests, we recommend $n_\text{part}$ to be set between $3$ and $10$, $\delta_\text{re}$ between $0.1$ (for expensive gradients) and $1$, $\underline{n_{\text{samp}}}$ of $10$, and $p_\text{re}=5 \%.$
\end{rmk}

\subsection{Dropping iterations of the \texttt{MICE} hierarchy}\label{sc:mice.dropping}

Given our estimator's stochastic nature, at the current iteration $k$, we may wonder if the iteration $k -1$ should be kept or dropped out from the \texttt{MICE} hierarchy since it may not reduce the computational work.
The procedure we follow here draws directly from an idea introduced by Giles \cite{Gil15} for the \texttt{MLMC} method.
Although the numerical approach is the same, we construct the algorithm in a greedy manner. We only check the case of dropping the previous iteration in the current hierarchy. In this approach, we never drop the initial iteration $\ell_0$.
We begin this analysis by noting that the cost of evaluating $\Delta_{\ell, k, \alpha}$ for $\ell>\min \cl{L}_k$ is of $2$ evaluations of $\nabla_{\bs{\xi}} f$.
In the optimal setting \eqref{eq:optimal.sampling.1}, we have that $M^*_{\ell,k}\propto\sqrt{V_{\ell,k}}$ for all $\ell>\min \cl{L}_k$. Next, considering two consecutive iterations $k$ and $k-1$, we arrive at the following relation
\begin{equation}\label{eq:optimal.M.propto}
M^*_{k,k}\sqrt{\dfrac{1}{V_{k,k}}} \approx M^*_{k-1,k}\sqrt{\dfrac{1}{V_{k-1,k}}}.
\end{equation}
The combined variance of the estimator relating these two terms, for $k$ and $k-1$, leads to
\begin{equation}\label{eq:combined.variance}
\dfrac{V_{k,k}}{M^*_{k,k}}+\dfrac{V_{k-1,k}}{M^*_{k-1,k}} \approx \dfrac{V_{k,k}}{M^*_{k,k}}\left(1+\sqrt{\dfrac{V_{k-1,k}}{V_{k,k}}}\right),
\end{equation}
while the combined cost is
\begin{equation}\label{eq:combined.cost}
M^*_{k,k}+M^*_{k-1,k} \approx M^*_{k,k}\left(1+\sqrt{\dfrac{V_{k-1,k}}{V_{k,k}}}\right),
\end{equation}
and the product of the combined cost and variance yields
\begin{equation}\label{eq:combined.variance.cost}
V_{k,k}\left(1+\sqrt{\dfrac{V_{k-1,k}}{V_{k,k}}}\right)^2.
\end{equation}

Next, consider the variance of the difference related to the pair $(k,p_k(k-1))$ as $\bar{V}_{k,k}$. By comparing the variance of the pairs $(k,k-1)$ and $(k,p_k(k-1))$, we arrive at
\begin{equation}\label{eq:inequality.drop}
\bar{V}_{k,k}<V_{k,k}\left(1+\sqrt{\dfrac{V_{k-1,k}}{V_{k,k}}}\right)^2,
\end{equation}
the inequality condition expressing that, if true, using the pair $(k,k-1)$ is computationally more expensive than using the pair $(k,p_k(k-1))$.

Finally, the inequality condition \eqref{eq:inequality.drop}, is written algorithmically as
\begin{equation}\label{eq:drop.condition}
\begin{split}
\texttt{if} &\quad \dfrac{\bar{V}_{k,k}}{(\sqrt{V_{k-1,k}} + \sqrt{V_{k,k}})^2} \le 1+\delta_{\mathrm{drop}}, \\
\texttt{then} &\quad \text{drop iteration $k-1$ of \texttt{MICE}} \colon\quad\cl{L}_k\leftarrow\cl{L}_k\!\setminus\{k-1\},
\end{split}
\end{equation}
where $\delta_{\mathrm{drop}}\ge 0$ is a real parameter that favors dropping.
The larger $\delta_{\mathrm{drop}}$, the smaller the resulting hierarchy $\cl{L}_k$ and the storage memory overhead.
We recommend values between $0.1$ and $0.5$ for $\delta_{\mathrm{drop}}$.

\subsection{Restarting the \texttt{MICE} hierarchy}\label{sc:mice.restart}
As we verified in the previous section on whether or not we should keep the iteration $\ell = k -1$ in the \texttt{MICE} hierarchy, we also may wonder if restarting the estimator may be less expensive than updating it.
Usually, in the literature of control variates techniques for stochastic optimization, the restart step is performed after a fixed number of iterations; see, for instance, \cite{Joh13,defazio2014saga,Nit14,Kon16}.
Here, the restart is instead done to minimize the work required to achieve the desired tolerance $\epsilon$.

Given the optimal work $\cl{W}(\{M_{\ell,k}^*\}_{\ell\in\cl{L}_k})$ (for $M_{\ell,k}^*$ as presented in \eqref{eq:optimal.sampling.1} and \eqref{eq:optimal.sampling.2}), the additional work required to perform iteration $k$ is given by
$\Delta \cl{W}\left(
\{M_{\ell,k}^*\}_{\ell\in\cl{L}_k}
\right)$.

To decide whether to restart or not, we compare it with the cost of restarting at $k$,
\begin{equation}
  \cl{W}_{\text{rest}, k}^* =
  M_{\text{rest},k}^* C_\nabla=
  \left\lceil
  \frac{
  \bb{V}[\norm{\nabla_{\bs{\xi}} f(\bs{\xi}_k, \theta)}{} |\{\bs{\xi_\ell}\}_{\ell\in\cl{L}_k}]
  }{
  \epsilon^2
  \norm{\nabla_{\bs{\xi}} F(\bs{\xi_k})}{}^2
  }
  \right\rceil C_\nabla.
\end{equation}
Thus,
\begin{equation}\label{eq:restart.condition}
\begin{split}
\texttt{if} &\quad \dfrac{
\cl{W}_{\text{rest}, k}^*
}{
\Delta \cl{W}\left(
\{M_{\ell,k}^*\}_{\ell\in\cl{L}_k}
\right)
}
\le 1+\delta_{\mathrm{rest}}, \\
\texttt{then} &\quad \text{restart \texttt{MICE} at } \bs{\xi}_k\colon\quad\cl{L}_k\leftarrow\{k\},
\end{split}
\end{equation}
for a real parameter $\delta_{\mathrm{rest}}\ge 0$ that favors restarting.
The larger $\delta_{\mathrm{rest}}$, the smaller the resulting hierarchy $\cl{L}_k$, which may be desired to decrease the storage memory.
Note that if a restart is performed, then the index set is also restarted with the current iteration, that is, $\cl{L}_k\leftarrow\{k\}$.
In all cases, we use the resampling presented in Remark \ref{rem:grad_resampling} to get an estimate of the norm of the gradient at iteration $k$, i.e., $\norm{\nabla_{\bs{\xi}} F(\bs{\xi_k})}{} \approx \norm{\nabla_{\bs{\xi}}\cl{F}^{\text{re}}_k}{}$.
We recommend values between $0$ and $0.3$ for $\delta_{\mathrm{rest}}$.

\subsection{Clipping the \texttt{MICE} hierarchy} \label{sc:mice.clipping}

In some cases, it may be advantageous to discard only some initial iterations and not the whole hierarchy. We refer to this procedure as clipping the hierarchy.
We propose two different approaches to decide when and where to clip the hierarchy.

\subsubsection{Clipping ``A''} \label{sc:mice.clipping.A}

Here, for a given hierarchy $\cl{L}_k$, we check, for every $\ell \in \cl{L}_k$, what would be the cost of clipping the hierarchy at $\ell$ and choose the case where the cost is minimized.
For that sake, we create a family of clipped hierarchies $\cl{S} = \left\{\cl{L}_k^{\ge \ell}\right\}_{\ell \in \cl{L}_k}$, where

\begin{equation}
  \cl{L}_k^{\ge \ell} \coloneqq \left\{\ell' \in \cl{L}_k: \ell' \ge \ell\right\}.
\end{equation}
For each candidate clipped hierarchy $\cl{L}_k^{\ge \ell} \in \cl{S}$, we calculate the optimal sample sizes using either \eqref{eq:optimal.sampling.1} or \eqref{eq:optimal.sampling.2},
\begin{equation}
  \left\{
  \{M_{\ell,k}^*\}_{\ell\in
  \cl{L}_k^{\ge \ell}}
  \right\}_{
  \cl{L}_k^{\ge \ell}
  \in
  \cl{S}
  }.
\end{equation}
Then, the best clipped hierarchy at iteration $k$ is chosen as the one that minimizes the additional work,
\begin{equation}
  \cl{L}_k^{\ge \ell^*}
   = \underset{\cl{L}_k^{\ge \ell} \in \cl{S}}{\text{arg min}} \,
  \Delta\cl{W}
  \left(
  \cl{L}_k^{\ge \ell}
  \right),
\end{equation}
where we abuse notation by suppressing the dependency of $\Delta\cl{W}$ on $M_{\ell, k}^*$, noting that $M_{\ell, k}^*$ is completely defined by $\cl{L}_k$ using \eqref{eq:optimal.sampling.1} or \eqref{eq:optimal.sampling.2}.
This clipping technique can be applied in both the continuous and discrete cases.

\subsubsection{Clipping ``B''} \label{sc:mice.clipping.B}
Another alternative is, for finite sets of size $N$, to clip the hierarchy whenever $M^*_{\ell, k } = N$ for $\ell \in \cl{L}_k$, thus

\begin{equation}
  \ell^* = \max\left\{\ell \in \cl{L}: M_{\ell, k}^* = N\right\},
\end{equation}
then,
\begin{equation}
  \cl{L}_k = \cl{L}_k^{\ge \ell^*}.
\end{equation}

Clipping ``A'' adds an extra computation overhead when calculating $M_{\ell, k}$ for each $\ell \in \cl{L}_k$ each iteration $k$.
However, clipping shortens the hierarchy, thus possibly reducing the general overhead of \texttt{MICE}.
Moreover, clipping the hierarchy may reduce the frequency of restarts and the bias of the \texttt{MICE} estimator.

\section{Idealized \texttt{SGD-MICE} convergence and gradient sampling cost analysis} \label{sc:mice.conv_and_work}

As shown in Remarks \ref{rem:mice.bias} and \ref{rem:mice.stat_error}, \texttt{MICE} estimator has a bias from sampling gradients at $\bs{\xi}_\ell$ for $\ell<k$ when we condition the \texttt{MICE} estimator on $\bs{\xi}_k$.
In the analysis of the convergence and gradient sampling cost in this section, we ignore this bias, which we refer to as idealized \texttt{MICE}.
Later, in the numerical section, we compare the performance of \texttt{MICE} and idealized \texttt{MICE} to show that the effect of the bias is negligible.
To perform this comparison, we created an algorithm that satisfies all the conditions used to prove the convergence of \texttt{MICE}, which we present in Appendix \ref{sc:ideal_mice}.

\subsection{Idealized \texttt{SGD-MICE} convergence analysis}\label{sc:mice.convergence}
We now aim at establishing the rate of convergence of the idealized version of \texttt{SGD-MICE}. Here, we crucially assume that the relative control of the statistical error in the mean gradient estimation is satisfied; see Assumption \ref{as:grad_rel_err}.
Additionally, we make some standard assumptions on the strong convexity and regularity of the objective function $F$.
\begin{ass}[$L$-Lipschitz gradient]\label{as:lipschitz}
The gradient of $F \colon \bb{R}^{d_{\bs{\xi}}}\mapsto\bb{R}$ is $L$-Lipschitz, for some $L > 0$, then
\begin{equation}\label{eq:lipschitz}
\norm{\nabla_{\bs{\xi}} F(\bs{x}) - \nabla_{\bs{\xi}} F(\bs{y})}{} \le L \norm{\bs{y} - \bs{x}}{}, \qquad \forall \bs{x}, \bs{y} \in \Xi.
\end{equation}
\end{ass}
\begin{ass}[Strongly convex]\label{as:strongly.convex}
If $F$ is $\mu$-strongly convex, then, for some $\mu>0$,
\begin{equation}\label{eq:strongly.convex}
F(\bs{y}) \ge F(\bs{x}) + \langle \nabla_{\bs{\xi}} F(\bs{x}), \bs{y} - \bs{x} \rangle + \dfrac{\mu}{2} \norm{\bs{y} - \bs{x}}{}^2, \qquad \forall\,\bs{x}, \bs{y} \in \Xi.
\end{equation}
\end{ass}
Assumptions \ref{as:lipschitz} and \ref{as:strongly.convex} allow for a simple, idealized stopping criterion at iteration $k^* =k^*(tol)$, namely that $k^*$ is the minimum integer that satisfies
\begin{equation}\label{eq:stopping_ideal}
\norm{\nabla_{\bs{\xi}}F(\bs{\xi}_{k^*})}{} < \mu  \, tol.
\end{equation}
\begin{ass}[Gradient variability control]\label{as:grad_rel_err}
Assume that the update of \texttt{SGD} \eqref{eq:sgd} uses a conditionally unbiased mean gradient estimator $\bs{\upsilon}_k
$, satisfying
\begin{equation}\label{eq:grad_ass_unbiased}
 \bb{E}[\bs{\upsilon}_k |\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}] = \nabla_{\bs{\xi}} F(\bs{\xi}_k)
 \end{equation}
and, for a given value of $\epsilon>0,$
\begin{equation}\label{eq:grad_rel_err}
 \bb{E}[\norm{\nabla_{\bs{\xi}} F(\bs{\xi}_k)-\bs{\upsilon}_k}{}^2|\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}]
 \le\norm{\nabla_{\bs{\xi}} F(\bs{\xi}_k)}{}^2\epsilon^2.
 \end{equation}
\end{ass}
\begin{prop}[Idealized \texttt{SGD-MICE} convergence]\label{prp:convergence.mice}
If Assumptions \ref{as:lipschitz}, \ref{as:strongly.convex} and \ref{as:grad_rel_err} hold, with a constant step-size given by
\begin{equation}\label{eq:step.size.strongly.convex}
\eta = \dfrac{2}{(L+\mu)(1 + \epsilon^2)},
\end{equation}
the \texttt{SGD} method given by \eqref{eq:sgd}
 enjoys linear convergence rate, namely
\begin{align}\label{eq:residual.recursion}
\bb{E}[\norm{\bs{\xi}_{k+1} - \bs{\xi}^*}{}^2] &\le \left[
\frac{\left(\dfrac{\kappa-1}{\kappa+1}\right)^2+\epsilon^2}{1+\epsilon^2}
\right]^k\bb{E}[\norm{\bs{\xi}_0 - \bs{\xi}^*}{}^2].
\end{align}
Here, $\epsilon$ is the relative statistical tolerance, while the Lipschitz gradient constant $L$ is given in \eqref{eq:lipschitz}, the strongly convex constant $\mu$ in \eqref{eq:strongly.convex}, and the condition number $\kappa=L/\mu$.
\end{prop}
\begin{proof}
From the update of \texttt{SGD} \eqref{eq:sgd} with constant step-size and the conditionally unbiased mean gradient estimator $\bs{\upsilon}_k
$, we have that
\begin{align}
\norm{\bs{\xi}_{k+1} - \bs{\xi}^*}{}^2 &= \norm{\bs{\xi}_k - \bs{\xi}^* - \eta\bs{\upsilon}_k}{}^2,\nonumber\\
&= \norm{\bs{\xi}_k - \bs{\xi}^*}{}^2 - 2 \eta \lrang{\bs{\xi}_k - \bs{\xi}^*}{\bs{\upsilon}_k} + \eta^2\norm{\bs{\upsilon}_k}{}^2.
\end{align}
Taking  conditional expectation  on $\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}$ and using \eqref{eq:grad_ass_unbiased}, we arrive at
\begin{align}
\bb{E}[\norm{\bs{\xi}_{k+1} - \bs{\xi}^*}{}^2|\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}] &
= \bb{E}[\norm{\bs{\xi}_k - \bs{\xi}^* - \eta\bs{\upsilon}_k}{}^2|\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}]\nonumber\\
&= \norm{\bs{\xi}_k - \bs{\xi}^*}{}^2 - 2 \eta \lrang{\bs{\xi}_k - \bs{\xi}^*}{\nabla_{\bs{\xi}} F(\bs{\xi}_k)} + \eta^2\bb{E}[\norm{\bs{\upsilon}_k}{}^2|\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}].
\end{align}
Given Assumptions \ref{as:lipschitz} and \ref{as:strongly.convex}, the coercivity of the gradient classical theorem (see the book by Nesterov \cite[Theorem 2.1.12]{nesterov2018lectures}) reads
\begin{equation}\label{eq:lipschitz.strongly.convex}
\lrang{\nabla_{\bs{\xi}} F(\bs{\xi}_k) - \nabla_{\bs{\xi}} F(\bs{\xi}^*)}{\bs{\xi}_k - \bs{\xi}^*} \ge \frac{\mu L}{\mu + L} \norm{\bs{\xi}_k - \bs{\xi}^*}{}^2 + \frac{1}{\mu + L} \norm{\nabla_{\bs{\xi}} F(\bs{\xi}_k) - \nabla_{\bs{\xi}} F(\bs{\xi}^*)}{}^2.
\end{equation}
Combining the two previous expressions, yields
\begin{multline}\label{eq:interm}
\bb{E}[\norm{\bs{\xi}_{k+1} - \bs{\xi}^*}{}^2|\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}] \le \left(1-\dfrac{2 \eta \mu L}{L+\mu}\right)\norm{\bs{\xi}_k - \bs{\xi}^*}{}^2
- \dfrac{2\eta}{L+\mu}\norm{\nabla_{\bs{\xi}} F(\bs{\xi}_k)}{}^2 + \eta^2\bb{E}[\norm{\bs{\upsilon}_k}{}^2|\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}].
\end{multline}
Now, using  the identity
\begin{equation}
\bb{E}[\norm{\nabla_{\bs{\xi}} F(\bs{\xi}_k)-\bs{\upsilon}_k}{}^2|\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}]=-\norm{\nabla_{\bs{\xi}} F(\bs{\xi}_k)}{}^2+\bb{E}[\norm{\bs{\upsilon}_k}{}^2|\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}],
\end{equation}
and the \texttt{MICE} control on the relative variability of the gradient stated in Assumption \ref{as:grad_rel_err}, \eqref{eq:grad_rel_err}
\begin{equation}
\bb{E}[\norm{\nabla_{\bs{\xi}} F(\bs{\xi}_k)-\bs{\upsilon}_k}{}^2|\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}]\le\norm{\nabla_{\bs{\xi}} F(\bs{\xi}_k)}{}^2\epsilon^2,
\end{equation}
together with \eqref{eq:interm} we obtain
\begin{multline} \label{eq:conv_proof_1}
\bb{E}[\norm{\bs{\xi}_{k+1} - \bs{\xi}^*}{}^2|\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}] \le \left(1-\dfrac{2 \eta \mu L}{L+\mu}\right)\norm{\bs{\xi}_k - \bs{\xi}^*}{}^2
+ \eta\left(-\dfrac{2}{L+\mu}+\eta(1+\epsilon^2)\right)\norm{\nabla_{\bs{\xi}} F(\bs{\xi}_k)}{}^2.
\end{multline}
Setting the step-size to be the largest step-size that guarantees linear convergence in the $L^2$ sense, we get
\begin{equation} \label{eq:opt_step}
\eta=\dfrac{2}{(L+\mu)(1+\epsilon^2)}.
\end{equation}
Substituting this step-size in \eqref{eq:conv_proof_1}, we obtain
\begin{equation}
\bb{E}[\norm{\bs{\xi}_{k+1} - \bs{\xi}^*}{}^2|\{\bs{\xi}_\ell\}_{\ell\in\cl{L}_k}] \le \left(\left(\dfrac{\kappa-1}{\kappa+1}\right)^2+\epsilon^2\right)\frac{1}{1+\epsilon^2}\norm{\bs{\xi}_k - \bs{\xi}^*}{}^2.
\end{equation}
Next, with the tower property and by unrolling the recursion, we conclude this proof.
\end{proof}

\begin{rmk}[Canonical convergence for infinite populations of \emph{strongly convex} functions]\label{rem:SGD.comparison}
For \emph{strongly convex} functions, the classical convergence analysis of the Vanilla \texttt{SGD} is sublinear $\propto k^{-1}$, see the book by Nemirovski \cite[Proposition 14.1.2]{nemirovski1994efficient}, whereas the idealized \texttt{SGD-MICE} with constant step-size \eqref{eq:step.size.strongly.convex} attains linear convergence as in \eqref{eq:residual.recursion}.
Moreover, observe that as  $\epsilon\to0$, the idealized \texttt{SGD-MICE} recovers the convergence of the deterministic Gradient Descent (\texttt{GD}) method; see the book by Nesterov \cite[Theorem 2.1.13]{nesterov2018lectures} for deterministic convex functions.
However, as we show in the \S\ref{sc:mice.work} below, it is not advantageous to take $\epsilon$ too small in  \texttt{SGD-MICE}.
\end{rmk}

\subsection{Idealized \texttt{SGD-MICE} gradient sampling cost analysis}\label{sc:mice.work}

\begin{description}
\item[Introduce the gradient sampling cost for \texttt{SGD-MICE}]
For the sake of simplicity and given the cumulative nature of the computational gradient sampling cost in \texttt{MICE}, we analyze the total gradient sampling cost on a set of iterations $\{\bs{\xi}_\ell\}_{\ell=0}^{k^*}$ converging to $\bs{\xi}^\ast$ as per Proposition \ref{prp:convergence.mice}.
Observe that in this simplified setting, the number of iterations required to stop the iteration, $k^*= k^*(tol)$, and both the sequences $(\bs{\xi}_\ell)$ and $(M_{\ell,k})$  are still random. Indeed, we define
$$
k^*= \min\{k\ge0\colon \|\bs{\xi}_k-\bs{\xi}^*\|\le tol\}.
$$
Proposition \ref{prp:convergence.mice} establishes the $L^2$ convergence rate of the distance to the optimal point $\bs{\xi}^\ast$. The following corollary gives a similar rate that holds a.s.
\begin{cor}\label{cor:as_rate_conv}
Under the same assumptions in Proposition \ref{prp:convergence.mice},
define $0<\delta,\alpha<1$ by
\begin{equation}\label{eq:delta_def}
\delta\coloneqq 1-\left(\dfrac{\kappa-1}{\kappa+1}\right)^2
\end{equation}
and
\begin{equation}\label{eq:alpha_def}
\alpha = \frac{\delta}{1+\epsilon^2},
\end{equation}
respectively.
Then, for $0<\beta<\alpha$, we have, a.s. that
\begin{equation}\label{eq:as_rate}
\|\bs{\xi}_k - \bs{\xi}^*\| \le C (1-\beta )^{k/2} \text{ for all }k \text{ with $C$ not depending on }k.
\end{equation}
\end{cor}
\begin{proof}
Observe that from Proposition \ref{prp:convergence.mice} we have
\begin{equation}
\bb{E}[\|\bs{\xi}_k - \bs{\xi}^*\|^2] \le \left(1-\alpha\right)^k \bb{E}[\|\bs{\xi}_0 - \bs{\xi}^*\|^2].
\end{equation}
Define,  the sequence
\begin{equation}
a_k \coloneqq \frac{\|\bs{\xi}_k - \bs{\xi}^*\|}{ (1-\beta )^{k/2}}.
\end{equation}
Then we have that  $\sum_k \bb{E}[a_k^2]<\infty$ and that $\sum_k \bb{E}[a_k]\le \sum_k \bb{E}[a_k^2]^{1/2}<\infty$  so Kolmogorov's two series theorem implies that
\begin{equation}
\sum_k a_k<\infty, \text{ a.s.},
\end{equation}
and from here, we conclude that $a_k$  is bounded a.s., implying \eqref{eq:as_rate}.
\end{proof}
Observe that for $\kappa\gg1$, $\alpha$ defined in Corollary \ref{cor:as_rate_conv} satisfies
\begin{equation}
\alpha \sim \frac{4}{\kappa\left({1+\epsilon^2}\right)} \ll 1.
\end{equation}
Also, we thus have that the idealized \texttt{SGD-MICE} has linear convergence a.s., meaning that the number of required iterations to stop
satisfies $k^*=\cl{O}(\log(tol^{-1}))$ a.s.
We can actually characterize the distribution of $k^*$ as follows:
\begin{cor}[\texttt{SGD-MICE} number of iterations]\label{cor:convergence.seq}
If assumptions in Proposition \ref{prp:convergence.mice} hold then, letting
\begin{equation}
k_1 \coloneqq \frac{\log(tol^{-2}\bb{E}[\|\bs{\xi}_0 - \bs{\xi}^*\|^2])}{\log(1/(1-\alpha))},
\end{equation}
we have
\begin{equation}\label{eq:distrib}
\bb{P}[k^* \ge k ] \le \left\{
\begin{aligned}
1, & \text{ if }k< k_1 \\
(1-\alpha)^{ k-k_1}& \text{ otherwise.}
\end{aligned}
\right.
\end{equation}
Moreover, we have
\begin{equation}\label{eq:expect}
\bb{E}[k^* ] \le \frac{1}{\alpha} + \max\left\{0,\frac{\log(tol^{-2}\bb{E}[\|\bs{\xi}_0 - \bs{\xi}^*\|^2])}{\log(1/(1-\alpha))}\right\}.
\end{equation}
\end{cor}
\begin{proof}
First observe that
\begin{equation}
\bb{P}[k^* \ge k ] \le \bb{P}[\|\bs{\xi}_k - \bs{\xi}^*\| \ge tol ].
\end{equation}
Then apply Chebyshev's inequality and the exponential convergence in $L^2$  in Proposition \ref{prp:convergence.mice}, yielding
\begin{equation}\label{eq:distrib1}
\bb{P}[k^* \ge k ] \le \min\left\{1,
\left(1-\alpha\right)^k \bb{E}[\|\bs{\xi}_0 - \bs{\xi}^*\|^2] tol^{-2} \right\}.
\end{equation}
The result \eqref{eq:distrib} follows then directly.
To show \eqref{eq:expect}, simply use \eqref{eq:distrib1} and that
\begin{equation}\label{eq:distrib2}
\bb{E}[k^*] = \sum_{k\ge 0 }\bb{P}[k^* \ge k ] \le \max\{0,k_1\}+\frac{1}{\alpha}.
\end{equation}

\end{proof}

We now define the total gradient sampling cost, $\cl{C}_k= \cl{C}_k(tol)$, as the sum over all iterations of the number of samples, taking into account that on the base iteration the cost is halved, namely
\begin{equation}\label{eq:cost.estimation}
\cl{C}_k \coloneqq  M_0 + 2 \sum_{\ell=1}^{k} M_\ell.
\end{equation}
Now let us recall that $k^*$ is the iteration when we reach the desired accuracy and we stop. This means that
$\|\bs{\xi}_{k^*} - \bs{\xi}^*\|\le tol$ and that $\|\bs{\xi}_{k^*-1} - \bs{\xi}^*\|> tol$.

Evaluating the gradient sampling cost  \eqref{eq:cost.estimation} at iteration $k^*-1$ with the optimal sample sizes $M^*_\ell$ \eqref{eq:optimal.sampling.1}, and keeping in mind the enforced minimum number of samples per iteration $M_{min}\ge1$, we obtain
\begin{equation}\label{eq:cost}
\cl{C}_{k^*} \le 2\left(
\frac{1}{\epsilon^2 \norm{\nabla_{\bs{\xi}}F(\bs{\xi}_{k^*-1})}{}^2}
\left(\sum_{\ell \in \cl{L}_{k^*}} \sqrt{V_{\ell,k^*} } \right)^2\right) + M_{min} \,|\cl{L}_{k^*}|.
\end{equation}
Here in \eqref{eq:cost}, $V_{\ell,k}$ are the conditional variances defined in \eqref{eq:Velldef}.
Due to strong convexity and that $\|\bs{\xi}_{k^*-1} - \bs{\xi}^*\|\ge tol$ we have
\begin{equation}
\norm{\nabla_{\bs{\xi}}F(\bs{\xi}_{k^*-1})}{} \ge \mu  \, tol
\end{equation}
and combining this last inequality with \eqref{eq:cost} yields
\begin{equation}\label{eq:cost2}
\cl{C}_{k^*} \le 2
\left(
\frac{\sum_{\ell \in \cl{L}_{k^*}} \sqrt{V_{\ell,k^*} } }{\epsilon\,\mu \,tol}
\right)^2
 + M_{min} \,|\cl{L}_{k^*}|.
\end{equation}
To carry out the analysis of the total gradient sampling cost, we need to make an extra assumption that relates to the $L^2$ Lipschitz continuity of the gradients, namely
\begin{ass}[Strong Lipschitz gradient]\label{as:lipschitz2}
There exists $\hat L>0$ such that
\begin{equation}
 \bb{E}[\|\nabla_{\bs{\xi}} f(\bs{x},\bs{\theta}) - \nabla_{\bs{\xi}} f(\bs{y},\bs{\theta}) \|^2|\bs{x},\bs{y}] \le  (\hat L)^2 \|\bs{x} - \bs{y}\|^2, \text{for all }\bs{x},\bs{y}.
\end{equation}
\end{ass}
Observe that we have  $L\le \hat L$, with $L$ from Assumption \ref{as:lipschitz}. Indeed, Assumption \ref{as:lipschitz2} is the natural counterpart of the strong rate of convergence assumption that is standard in MLMC analysis, cf. \cite{Gil15}.

\begin{cor}[\texttt{SGD-MICE} total gradient sampling cost]\label{cor:work}
Consider \texttt{SGD-MICE} without Dropping or Restart.
If assumptions in Proposition \ref{prp:convergence.mice}  and Assumption \ref{as:lipschitz2} hold, we have, letting
\begin{equation}
\hat{\kappa}\coloneqq \hat L/\mu\ge \kappa,
\end{equation}
that
\begin{equation}\label{eq:expectw}
\begin{aligned}
\bb{E}[\cl{C}_{k^*} ] \le &  tol^{-2}\,2\,\left(\frac{2 \hat \kappa}{ \epsilon } \right)^2\frac{\bb{E}[\|\bs{\xi}_0 - \bs{\xi}^*\|^2]}{(1 - \left(1-\alpha\right)^{1/2})^2}
\\
& + M_{min} \left(
\frac{1}{\alpha} +  \max\left\{0,\frac{2\log(tol^{-1}) + \log(\bb{E}[\|\bs{\xi}_0 - \bs{\xi}^*\|^2])}{\log(1/(1-\alpha))}\right\}
\right).
\end{aligned}
\end{equation}
\end{cor}
\begin{proof}
To bound the expected value of the gradient sampling cost on \eqref{eq:cost2}, observe that, since there is neither Dropping nor Restart, $|\cl{L}_{k^*}|\le k^*.$ Then Corollary \ref{cor:convergence.seq} and \eqref{eq:cost} imply the second term in \eqref{eq:expectw}.
To bound the other term in \eqref{eq:cost}, we need to estimate the conditional variances using the pathwise Lipschitz constant of the gradients, i.e.
using Assumption \ref{as:lipschitz2} for some $\hat{L}>0$ we have
\begin{equation}
V_{\ell,k^*} \le \bb{E}[\|\nabla_{\bs{\xi}} f(\bs{\xi}_\ell,\bs{\theta}) - \nabla_{\bs{\xi}} f(\bs{\xi}_{\ell-1},\bs{\theta}) \|^2|\bs{\xi}_{\ell-1},\bs{\xi}_{\ell}] \le (\hat L)^2 \|\bs{\xi}_\ell - \bs{\xi}_{\ell-1}\|^2,
\end{equation}
and therefore
\begin{equation}
\sqrt{V_{\ell,k^*}} \le \hat{L} \|\bs{\xi}_\ell - \bs{\xi}_{\ell-1}\|,
\end{equation}
with $\bs{\xi}_{-1} =0.$
Thus
\begin{equation}
\sum_{\ell \in \cl{L}_{k^*}} \sqrt{V_{\ell,k^*} }  \le \hat{L} \sum_{0 \le \ell \le {k^*}}  \|\bs{\xi}_\ell - \bs{\xi}_{\ell-1}\| \le 2 \hat{L} \sum_{0 \le \ell \le {k^*}}  \|\bs{\xi}_\ell - \bs{\xi}^*\|
\end{equation}
To conclude, we need to estimate, using Proposition \ref{prp:convergence.mice},
\begin{align}
\bb{E}\left[\left(\sum_{0 \le \ell \le {k^*}}  \|\bs{\xi}_\ell - \bs{\xi}^*\|\right)^2\right] \le&{}
\sum_{0 \le \ell_1,\ell_2 < {\infty}}  \bb{E}[\|\bs{\xi}_{\ell_1} - \bs{\xi}^*\|\,\|\bs{\xi}_{\ell_2} - \bs{\xi}^*\|] \nonumber\\
\le &{} \bb{E}[\|\bs{\xi}_0 - \bs{\xi}^*\|^2] \, \left(\sum_{0 \le \ell < \infty} \left(1-\alpha\right)^{k/2} \right)^2 \nonumber\\
\le &{} \frac{\bb{E}[\|\bs{\xi}_0 - \bs{\xi}^*\|^2]}{\left(1 - \left(1-\alpha\right)^{1/2}\right)^2}.
\end{align}
Finally, combining the last two inequalities we arrive at
\begin{equation}
\bb{E}\left[
\left(
\frac{\sum_{\ell \in \cl{L}_{k^*}} \sqrt{V_{\ell,k^*} } }{\epsilon\,\mu \,tol}
\right)^2
\right]
\le
\left(\frac{2 \hat L}{\mu \, \epsilon \, tol} \right)^2\frac{\bb{E}[\|\bs{\xi}_0 - \bs{\xi}^*\|^2]}{(1 - \left(1-\alpha\right)^{1/2})^2},
\end{equation}
which, combined with \eqref{eq:cost2}, finishes the proof.
\end{proof}
Following the lines in {Corollary} \ref{cor:work}, a similar analysis can be carried out to estimate the variance of the gradient sampling cost, yielding that the variance is $\cl{O}(tol^{-4})$ with a proper characterization of the multiplicative constant.
\begin{rmk}[Optimal $\epsilon$]
Here we motivate the existence of an optimal choice for $\epsilon$.
Recall the expected gradient sampling cost in \eqref{eq:expectw} and that, for $\kappa\gg1$,
\begin{equation}
\alpha = \frac{\delta}{1+\epsilon^2}\sim \frac{4}{\kappa(1+\epsilon^2)}.
\end{equation}
Furthermore, with the approximations $\log(1+x)\sim x$ and
$(1 - \left(1-x\right)^{1/2})^2\sim ({x/2})^2$
 for $|x| \ll 1$ we approximate the upper bound in \eqref{eq:expectw} by
\begin{align}\label{eq:expect2}
\bb{E}[\cl{C}_{k^*} ] \le &{} tol^{-2}\,8\,\left(\frac{2 \hat \kappa}{ \epsilon } \right)^2\frac{\bb{E}[\|\bs{\xi}_0 - \bs{\xi}^*\|^2]}{\alpha^2}
\nonumber\\
&{} +\frac{M_{min}}{\alpha}  \left(
 1+ \max\{0,2\log(tol^{-1}) + \log(\bb{E}[\|\bs{\xi}_0 - \bs{\xi}^*\|^2])\}
\right) \nonumber\\
\lesssim &{}
  tol^{-2}\,\frac{1}{2}\,\left(\frac{2 \hat \kappa}{ \epsilon } \right)^2\,{\bb{E}[\|\bs{\xi}_0 - \bs{\xi}^*\|^2]}(\kappa(1+\epsilon^2))^2 \nonumber\\
& + \frac{\kappa(1+\epsilon^2)}{4} {M_{min}} \left(
 1+ \max\{0,2\log(tol^{-1}) + \log(\bb{E}[\|\bs{\xi}_0 - \bs{\xi}^*\|^2])\}
\right).
\end{align}
It can now be argued that the first term in \eqref{eq:expect2} always dominates the second term for all values of $\epsilon$. Thus, the optimal gradient sampling cost corresponds to the minimization of the first term with respect to $\epsilon$, yielding $\epsilon=1$ and the corresponding complexity for the idealized \texttt{SGD-MICE} to be $tol^{-2}.$
\end{rmk}
Indeed, substituting the optimal value of $\epsilon$ yields
\begin{equation}\label{eq:expect3}
\begin{aligned}
\bb{E}[\cl{C}_{k^*} ] \lesssim  & \, tol^{-2}\,8\left(\kappa{\hat \kappa}\right)^2\,{\bb{E}[\|\bs{\xi}_0 - \bs{\xi}^*\|^2]}
+ \frac{\kappa }{2} {M_{min}}  \left( 1 + \max\{0,2\log(tol^{-1}) + \log(\bb{E}[\|\bs{\xi}_0 - \bs{\xi}^*\|^2])\}
\right).\\
\end{aligned}
\end{equation}
\begin{rmk}[Stopping criterion]
In practice, applying the stopping criterion \eqref{eq:stopping_ideal} requires an approximation of the mean gradient norm at each iteration.
A natural approach is to use the  \texttt{MICE} estimator as such an approximation, yielding
\begin{equation}\label{eq:stopping_real}
\norm{\nabla_{\bs{\xi}}\cl{F}_{k^*}}{} < \mu  \, tol,
\end{equation}
provided that the statistical error in the mean gradient is controlled in a relative sense.
This quality assurance requires a certain number of gradient samples.
For example, let us consider the ideal case of stopping when we start inside the stopping region, near the optimal point $\bs{\xi}^\ast$.
To this end, suppose that the initial iteration point, $\bs{\xi}_0$, is such that $\|\bs{\xi}_0-\bs{\xi}^*\| \le tol$. What is the cost needed to stop by sampling gradients at $\bs{\xi}_0$ without iterating at all? Observing that we need a statistical tolerance $ \mu \, tol$, we thus need a number of samples $M$ that  satisfies
\begin{equation}
\frac{\bb{E}[\|\nabla_{\bs{\xi}} f(\bs{\xi}_0,\theta)\|^2]}{(\mu \, tol)^2} \le  M.
\end{equation}
Compare the last estimate with \eqref{eq:expect3}.
\end{rmk}
\end{description}

\section{\texttt{MICE} algorithm}\label{sc:mice_algorithm}

In this section, we describe the \texttt{MICE} algorithm.
In Algorithm \ref{alg:mice}, we present the pseudocode for the \texttt{MICE} estimator, with the verification for dropping iterations from the \texttt{MICE} hierarchy and with the restart technique.
Two coupling algorithms for the multi-iteration stochastic optimizers are presented in Appendix \ref{sc:algorithms}, these are \texttt{SGD-MICE} and \texttt{Adam-MICE}.
\begin{algorithm}[h]
 \setstretch{1.5}
 \caption{Pseudocode for \texttt{MICE} with dropping and restart techniques.
  \texttt{MICE} has as input the current iteration, $k$; the gradient of $f$, $\nabla_{\bs{\xi}} f$; the current point $\bs{\xi}_k$ where the expectation of $\nabla_{\bs{\xi}} f$ is to be estimated; some distribution $\pi$ from which $\bs{\theta}$ can be sampled; the statistical tolerance for the error in \texttt{MICE} estimates, $\epsilon$; and an initial sample size, $M_{min}$. Moreover, the history of the previous evaluations of \texttt{MICE}, the hierarchy of batches, is supposed to be known.
 }\label{alg:mice}
 \begin{algorithmic}[1]
  \Procedure{\texttt{MICE}}{$k$, $\nabla_{\bs{\xi}} f$, $\bs{\xi}_k$, $\pi$, $\epsilon$, $M_{min}$}
  \State $\cl{L}_k \gets \cl{L}_{k-1} \cup \{k\}$
  \State $\cl{I}_k \gets \{\alpha\}_{\alpha=1}^{M_{min}}$
  \State Sample $\boldsymbol{\theta}_\alpha \sim \pi \qquad\forall \alpha \in \cl{I}_k$
  \State Obtain $\Delta_{k \alpha}$ from \eqref{eq:Delta.sampled} for each $\alpha$
  \State Calculate $V_k$ from \eqref{eq:Vl}
  \State Evaluate $\bar{V}_k$ between the pairs $(k,k-2)$ using \eqref{eq:Vl}
  \If{condition \eqref{eq:drop.condition} is satisfied}

  \State $\cl{L}_k \gets \cl{L}_k \setminus \{k-1\}$
  \Comment{Drop last iteration}
  \EndIf
  \If{condition \eqref{eq:restart.condition} is satisfied}

  \State $\cl{L}_k \gets \{k\}$
  \Comment{Restart hierarchy}
  \EndIf
  \State{Clip the hierarchy using either ``A'' as in \S\ref{sc:mice.clipping.A} or ``B'' as in \S\ref{sc:mice.clipping.B}}
  \State Obtain $\norm{\nabla_{\bs{\xi}}\cl{F}_k^{\text{re}}}{}$ using the resampling as in Remark \ref{rem:grad_resampling}
  \While{$\sum_{\ell \in \cl{L}_k} \nicefrac{V_{\ell,k}}{M_\ell} \ge \epsilon^2\norm{\nabla_{\bs{\xi}}\cl{F}_k^{\text{re}}}{}^2$}
  \State Calculate $\left\{M^*_{\ell, k}\right\}_{\ell \in \cl{L}_l}$ from \eqref{eq:optimal.sampling.1} or \eqref{eq:optimal.sampling.2} using $\norm{\nabla_{\bs{\xi}}\cl{F}_k^{\text{re}}}{}^2$
  \For{$\ell \in \cl{L}_k$}
  \State $\Delta M_\ell = M^*_\ell - M_\ell$
  \State $\cl{I}'_\ell \gets \{\alpha\}_{\alpha=M_\ell+1}^{M_\ell + \Delta M_\ell}$
  \State Sample $\boldsymbol{\theta}_\alpha \sim \pi \qquad\forall \alpha \in \cl{I}'_\ell$
  \State Obtain $\Delta_{\ell, k, \alpha}$ from \eqref{eq:Delta.sampled} for each $\alpha \in \cl{I}'_\ell$
  \State Calculate $V_{\ell,k}$ from \eqref{eq:Vl}
  \EndFor
  \State $M_\ell \gets M^*_\ell$
  \EndWhile
  \State \textbf{return} $\nabla_{\bs{\xi}}\cl{F}_k = \sum_{\ell\in\cl{L}_k}\dfrac{1}{M_{\ell,k}^*}
  \sum_{\alpha\in\cl{I}_{\ell,k}}\Delta_{\ell,k,\alpha}$ from \eqref{eq:mice}
  \EndProcedure
 \end{algorithmic}
\end{algorithm}

In general, keeping all gradient realizations for all iterations in memory may be computationally inefficient, especially for large-dimensional problems.
To avoid this unnecessary memory overhead, we use Welford's online algorithm to estimate the variances $V_{\ell,k}$ online. We keep in memory only the samples mean and second-centered moments and update them in an online fashion \cite{welford1962note}.
Thus, with the Welford's online algorithm, the memory overhead is reduced to the storage of a single gradient, namely the mean, and a scalar, namely $M_\ell V_{\ell,k}$, for each iteration. This procedure makes the memory overhead much smaller than naively storing all gradients and evaluating variances when needed.
For first-order methods such as \texttt{Adam-MICE} and \texttt{SGD-MICE}, with a hierarchy $\cl{L}_k$ at iteration $k$, the memory overhead of \texttt{MICE} is of $| \cl{L}_k| (d_{\xi}+1)$ floating-point numbers.
Thus, for large-scale problems, dropping iterations and restarting the hierarchy are very important to reduce memory allocation.

\begin{rmk}[More general $\bs{\theta}$ probability distributions]\label{rem:theta_distribution}
Although in Assumption \ref{as:theta_distribution} we restricted our attention to the case where the probability distribution of $\bs{\theta}$, $\pi$, does not depend on $\bs{\xi}$, it is possible to use mappings to address more general cases. Indeed, let us consider the case where
\begin{equation}
\bs{\theta} = h(\tilde{\bs{\theta}}, \bs{\xi}),
\end{equation}
for some given smooth function $h$ and such that the distribution of $\tilde{\bs{\theta}}$, $\tilde \pi$, does not depend on $\bs{\xi}$.
Then we can simply write, letting $\tilde f(\bs{\xi},\tilde{\bs{\theta}}) = f(\bs{\xi},h(\tilde{\bs{\theta}},\bs{\xi}))$,
\begin{equation}
F(\bs{\xi}) = E[f(\bs{\xi},\bs{\theta})|\bs{\xi} ] = E[\tilde f(\bs{\xi},\tilde{\bs{\theta}})|\bs{\xi} ]
\end{equation}
and, by sampling $\tilde{\bs{\theta}}$ instead of $\bs{\theta}$, we are back in the setup of Assumption \ref{as:theta_distribution}.
\end{rmk}

\section{Numerical examples}\label{sc:numerics}

In this section, we present some numerical examples to assess the efficiency of \emph{Multi-Iteration Stochastic Optimizers}.
We focus on \texttt{SGD-MICE}, \texttt{Adam-MICE} and compare their performances with \texttt{SGD}, \texttt{Adam}, \texttt{SAG}, \texttt{SAGA}, \texttt{SRVG}, and \texttt{SARAH} methods in stochastic optimization.
When using \texttt{SGD}, with or without \texttt{MICE}, we assume the constants $L$ and $\mu$ to be known and use them to compute the step-size as in \eqref{eq:opt_step}.
As a measure of the performance of the algorithms, we use the optimality gap, which is the difference between the approximate optimal value at iteration $k$ and the exact optimal value,
\begin{equation}\label{eq:opt_gap}
F(\bs{\xi}_k)-F(\bs{\xi}^*).
\end{equation}
In some examples, we know the optimal value and optimal point analytically; otherwise, we estimate numerically by letting optimization algorithms run for many iterations.

As for \texttt{MICE} parameters, in each case we use different values for $\epsilon$, ranging from $0.5$ to $1.0$.
The other parameters are fixed for all problems, showing the robustness of \texttt{MICE} with respect to the tuning: $\delta_{\mathrm{drop}}=0.5$, $\delta_{\mathrm{rest}}=0$, $M_{min}$ is set to $5$ for general iterations and $50$ for restarts, and the maximum hierarchy size is set to $100$.
For the continuous cases, we use the clipping ``A'', whereas, for the finite case, we opted to use clipping ``B''.
As for the resampling parameters, we use $n_\text{part}=5$, $\delta_\text{re}=1.0$, and $p_\text{re}=5\%$ with a minimum resampling size of $10$.
We recommend using these parameters for general problems; however, $\delta_\text{re}$ can be reduced for problems with more computationally expensive gradients. In our numerical examples, we chose to report the run time taken by the different algorithms we tested in addition to the usual number of gradient evaluations. Although we believe that our \texttt{MICE} implementations, namely \texttt{SGD-MICE} and \texttt{Adam-MICE}, can be improved, we hope to show, that the different \texttt{MICE} versions are already robust and competitive.

\subsection{Stochastic quadratic function}\label{ex:1}

\FloatBarrier
This problem is a simple numerical example devised to test the performance of \texttt{SGD-MICE} on the minimization of a strongly convex function.
The function whose expected value we want to minimize is
\begin{equation}
  f(\bs{\xi}, \theta) = \frac{1}{2} \bs{\xi} \cdot \bs{H}(\theta) \cdot \bs{\xi}
  - \bs{b} \cdot \bs{\xi},
\end{equation}
where
\begin{equation}
 \bs{H}(\theta) \coloneqq
 \bs{I}_2(1-\theta) +
 \begin{bmatrix}
  2 \kappa & 0.5 \\
  0.5      & 1
 \end{bmatrix}
 \theta,
\end{equation}
$\bs{I}_2$ is the identity matrix of size $2$, $\bs{b}$ is a vector of ones, and $\theta \sim \cl{U}(0, 1)$.
We use $\kappa=100$ and initial guess
$\bs{\xi}_0=(20, 50)$.
The objective function to be minimized is
\begin{equation}\label{eq:quad.example}
  F(\bs{\xi}) = \frac{1}{2} \bs{\xi} \cdot
  \bb{E}[\bs{H}(\bs{\theta})] \cdot \bs{\xi}
  - \bs{b} \cdot \bs{\xi},
\end{equation}
where
\begin{equation}
 \bb{E}[\bs{H}(\theta)] =
 \begin{bmatrix}
  \kappa+0.5 & 0.25 \\
  0.25       & 1
 \end{bmatrix}.
\end{equation}
The optimal point of this problem is $\bs{\xi}^* = \bb{E}[\bs{H}(\theta)]^{-1} \bs{b}$.
To perform optimization using \texttt{SGD-MICE} and \texttt{SGD}, we use the unbiased gradient estimator
\begin{equation}
  \nabla_{\bs{\xi}} f(\bs{\xi}, \theta) = \bs{H}(\theta) \cdot \bs{\xi} - \bs{b}.
\end{equation}
We use the eigenvalues of the Hessian of the objective function, $ \bb{E}[\bs{H}(\theta)]$, to calculate $L$ and $\mu$ and thus define the step-size as in \eqref{eq:step.size.strongly.convex}.

In Figures \ref{fig:quad_per_k} and \ref{fig:quad_per_grad}, we present the optimality gap \eqref{eq:opt_gap}, the distance to the optimal point and the norm of the gradient estimate versus iteration and number of gradient evaluations, respectively.
It Figure \ref{fig:quad_per_grad} we also plot the iteration per number of gradient evaluations.
We mark the starting point, restarts, and ending point with blue, red, and purple squares, respectively; the dropped points with black $\times$, and the remaining iterations in the \texttt{MICE} hierarchy with green dots.
In Figure \ref{fig:quad_per_k}, we depict that \texttt{SGD-MICE} attains the $L^2$ convergence with a constant step-size, as predicted in Proposition \ref{prp:convergence.mice}.
In Figure \ref{fig:quad_per_grad}, we present the convergence plots versus gradient evaluations, exhibiting numerical rates of $\cl{O}(M^{-\nicefrac{1}{2}})$ for the distance to the optimal point and for the norm of the gradient and of $\cl{O}(M^{-1})$ for the optimality gap.
These rates are expected as the distance to the optimal point converges linearly (see Proposition \ref{prp:convergence.mice}) and the cost of sampling new gradients per iteration grows as $\cl{C} = \cl{O}(\| \nabla F(\bs{\xi}_k)\|^{-2})$, as shown in \eqref{eq:optimal.sampling.1}.
Consequently, since the objective function is quadratic, the convergence of the optimality gap is $\cl{O}(M^{-1})$.
Note that the convergence is exponential as in the deterministic case until around $2 \times 10^4$ gradient evaluations.
After this point, there is a change of regime in which we achieve the asymptotic rates.
We note that, after $2 \times 10^4$, the cost of performing each iteration grows exponentially.

Figure \ref{fig:quad_err} presents the relative error in gradient evaluations, the $V_{k,k}$, and the length of the hierarchy versus the number of iterations.
Together with the relative error of the gradients, we plot the given tolerance on the statistical error, $\epsilon$.
It can be observed that \texttt{SGD-MICE} is successful in keeping the relative error below $\epsilon$
and that most of the iterations are dropped with the hierarchy achieving a length of 11 after $1170$ iterations.
Also, at around iterations $157$ and $641$, clipping can be observed, i.e., the hierarchy length is reduced without restarting it.
Note that, when the first clipping is done, there is a significant reduction on the relative error.
Once the relative statistical error becomes of size $\epsilon$, around iteration $k=500$ and $2 \times 10^4$ gradient evaluations, the exponential convergence in Figure \ref{fig:quad_per_grad} is lost.
The \texttt{SGD-MICE} algorithm drops all iterations from $180$ until $550$.
When the relative error achieves $\epsilon$, it becomes advantageous to keep some iterations in the hierarchy, consequently progressively reducing the $V_{k, k}$.
In the asymptotic regime, after iteration $1000$, two jumps can be seen in the $V_{k,k}$ exactly at the iterations where iterations are not dropped.

To verify if \texttt{SGD-MICE} indeed achieves the $L^2$ convergence predicted in \S\ref{sc:mice.convergence}, we perform a series of $1000$ independent runs for a few iterations arbitrarily close to the optimal point for \texttt{SGD-MICE}, idealized \texttt{SGD-MICE} (implemented as described in Appendix \ref{sc:ideal_mice}), and the deterministic Gradient Descent (\texttt{GD}). In all cases we use the \texttt{MICE} step-size for strongly convex functions as presented in \eqref{eq:step.size.strongly.convex}.
Idealized \texttt{SGD-MICE} is devised to validate the theory presented in \S\ref{sc:mice.convergence} when Assumptions \ref{as:lipschitz}, \ref{as:strongly.convex}, and \ref{as:grad_rel_err} are satisfied.
Specifically, we want to test if the bias in \texttt{MICE} gradient estimates affects the $L^2$ convergence of \texttt{SGD-MICE}.
In Figure \ref{fig:quad_L2}, we present the mean-squared distance to the optimal point convergence for \texttt{SGD-MICE}, idealized \texttt{SGD-MICE}, and \texttt{GD} with confidence intervals of $99\%$ for each method using bootstrapping.
For reference, we additionally plot the theoretical $L^2$ convergence for \texttt{SGD-MICE}.
Two cases of $\kappa$ are presented, $\kappa=10$ and $100$; and two cases of $\epsilon$ are presented, $\epsilon=0.1$ and $1.0$.
\texttt{SGD-MICE} performs as well as its idealized version, sometimes better.
The better performance of \texttt{SGD-MICE} might be a consequence of the measures taken to make the algorithm robust, for instance, the resampling and $M_{min}$.

To validate the robustness and performance increase of \texttt{SGD-MICE}, we compare it with \texttt{SGD} using Monte Carlo sampling, at every iteration, to estimate the mean gradient, which we call \texttt{SGD-A}.
Likewise \texttt{SGD-MICE}, \texttt{SGD-A} controls the relative statistical error in the mean gradient approximation and uses the resampling technique to estimate its norm.

In Figure \ref{fig:quad_mice_mc}, we present the optimality gap per iteration and number of gradient evaluations for \texttt{SGD-A}, \texttt{SGD-MICE}, and \emph{vanilla} \texttt{SGD} (Robbins--Monro algorithm \cite{robbins1951stochastic}).
\texttt{SGD-MICE} appears to be more robust than \texttt{SGD-A}: the latter starts to oscillate once it gets close to the optimal point and thus is unable to maintain the linear convergence per iteration.
Indeed, we performed $100$ independent runs of \texttt{SGD-A} and verified that, in all the runs, it oscillates in the same way as it gets close to the optimal point.
Also, \texttt{SGD-MICE} requires less gradient evaluations to reduce the optimality gap than \texttt{SGD-A}.
For \emph{vanilla} \texttt{SGD}, we used a decreasing step-size of $\eta_k = (L(1+k/50))^{-1}$, tuned for the best performance of the method.
When it comes to \texttt{SGD}, \texttt{SGD-MICE} performed better, achieving a lower optimality gap by several orders of magnitude for the same number of gradient evaluations.

Although this example is very simple, it illustrates the performance of \texttt{SGD-MICE} in an ideal situation where both $L$ and $\mu$ are known.
Most importantly, we observed that the simplified theory presented in \S\ref{sc:mice} and \S\ref{sc:mice.conv_and_work} still describes  the behavior of \texttt{SGD-MICE} well; the linear $L^2$ convergence is achieved for the optimal step-size even with an $\epsilon=1.0$.
Finally, \texttt{SGD-MICE} was able to automatically decide whether to drop iterations, restart, or clip the hierarchy to minimize the overall work required to attain the linear convergence per iteration.
\begin{rmk}[Alternative quadratic example]
Here, we aim instead at minimizing the expected value of
\begin{equation}
f(\bs{\xi},\bs{\theta})\coloneqq(\bs{\xi}+\theta_0\mbox{\tenbfsl 1\/})\cdot\bs{H}\cdot(\bs{\xi}+\theta_1\mbox{\tenbfsl 1\/})\qquad\text{with}\qquad
\bs{H} =
\begin{bmatrix}
100 & 3\\
3 & 8
\end{bmatrix}.
\end{equation}
Here, $\mbox{\tenbfsl 1\/}$ is a vector fully populated with ones, $\theta_0$ and $\theta_1$ are two independent, Gaussian $N(0,\sigma^2)$ distributed, random variables. In this example, we performed numerical experiments, with $\sigma$ ranging from $10^{-4}$ to $10^{-1}$, to assess the efficiency of the couplings \texttt{SGD-MICE} and \texttt{Adam-MICE}. First, note that for two different values of $\bs{\xi}$, there is perfect correlation between the corresponding gradients, meaning that the variance of the gradient differences is zero. As a consequence, dropping occurs for all intermediate iterations and we only need one gradient difference sample at the current iteration point, $\bs{\xi}_k$. As we approach the optimal point $\bs{\xi}^*,$ only the batch size at the initial point $\bs{\xi}_0$, that is, $M_{0,k}$, increases along the stochastic optimization process.
In this example,  \texttt{MICE} performs optimally and outperforms other alternatives
\end{rmk}

\begin{figure}
  \includegraphics[width=.9\linewidth]{SGD_MICE_all_per_iter.pdf}
  \caption{
  Single run, stochastic quadratic example, Equation \eqref{eq:quad.example} with $\kappa=100$ and $\epsilon=1$.
  Optimality gap (top), distance to the optimal point (center), and norm of gradient estimate (bottom) per iteration for \texttt{SGD-MICE}.
  The starting point, the restarts, and the end are marked respectively as blue, red, and purple squares, iterations dropped with black $\times$, and the remaining \texttt{MICE} points with green circles.
  \texttt{SGD-MICE} is able to achieve linear the $L^2$ convergence as predicted in Proposition \ref{prp:convergence.mice}.
  }
  \label{fig:quad_per_k}
\end{figure}

\begin{figure}
  \includegraphics[width=.9\linewidth]{SGD_MICE_all_per_grads.pdf}
  \caption{Single run, stochastic quadratic example, Equation \eqref{eq:quad.example} with $\kappa=100$ and $\epsilon=1$.
  Optimality gap (top), distance to the optimal point (center top), norm of gradient estimate (center bottom), and number of iterations (bottom) per number of gradient evaluations for \texttt{SGD-MICE}.
  The starting point, the restarts, and the end are marked respectively as blue, red, and purple squares, iterations dropped with black $\times$, and the remaining \texttt{MICE} points with green circles.
  The distance to the optimal point and the norm of the gradient achieve the Monte Carlo rate, $\cl{O}(M^{-\nicefrac{1}{2}})$ , whereas the optimality gap converges at the square of this rate, at $\cl{O}(M^{-1})$.
  }
  \label{fig:quad_per_grad}
\end{figure}

\begin{figure}
  \includegraphics[width=.9\linewidth]{SGD_MICE_err_vkk_chain_size.pdf}
  \caption{Single run, stochastic quadratic example, Equation \eqref{eq:quad.example} with $\kappa=100$ and $\epsilon=1$.
  Relative error in gradient estimates (top), estimated variances $V_{k,k}$ (center), and length of \texttt{MICE} hierarchy per iteration for \texttt{SGD-MICE} (bottom).
  The starting point, the restarts, and the end are marked respectively as blue, red, and purple squares, iterations dropped with black $\times$, and the remaining \texttt{MICE} points with green circles.
  The iterations where clipping is performed are marked with a red line.
  The tolerance on the statistical error is presented as a dashed line over the relative error in gradient estimates.
  }
  \label{fig:quad_err}
\end{figure}

\begin{figure}[h]
  \centering
   \includegraphics[width=\linewidth]{QuadC_L2.pdf}
   \caption{{Estimates based on $1000$ independent runs.} Stochastic quadratic example, Equation \eqref{eq:quad.example}.
   Here, we plot estimates for the $L^2$ error studied in Proposition \ref{prp:convergence.mice}. We carry out this for \texttt{SGD-MICE}, \texttt{SGD} with the idealized \texttt{MICE} (see Appendix \ref{sc:ideal_mice}), and gradient descent (\texttt{GD}) using the optimal \texttt{SGD-MICE} step-size for strongly convex functions, $\eta=\nicefrac{2}{(L+\mu)(1+\epsilon^2)}$.
   The idealized \texttt{MICE} is implemented satisfying the conditions as presented in Proposition \ref{prp:convergence.mice}.
   In this case, we use the constants $\kappa$ as $10$ (top) and $100$ (bottom) and $\epsilon$ as $0.1$ (left) and $1.0$ (right).
   \texttt{SGD-MICE} achieves the theoretically predicted $L^2$ convergence rate, as can be observed by comparing it to the idealized \texttt{SGD-MICE}, showing numerically that, for this example, the bias in \texttt{MICE} estimates is negligible and does not affect the $L^2$ linear convergence rate.
   In fact, due to the minimum sample size and the resampling, for $\epsilon=1.0$, \texttt{SGD-MICE} performed better than theoretically predicted.
   }
   \label{fig:quad_L2}
\end{figure}

\begin{figure}[h]
  \centering
   \includegraphics[width=.8\linewidth]{sgd.pdf}
   \caption{Single run, stochastic quadratic example, Equation \eqref{eq:quad.example} with $\kappa=100$ and $\epsilon=1$.
   For the sake of comparison, we solve this problem using \texttt{SGD-MICE}, \emph{vanilla} \texttt{SGD} (Robbins--Monro algorithm \cite{robbins1951stochastic}), and a version of \texttt{SGD} with fixed step-size and controlled variance, which we refer to as \texttt{SGD-A}.
   For \texttt{SGD-A}, we adaptively control the variance in gradient estimates by increasing the sample-size to keep the relative statistical error below a threshold, just as with \texttt{SGD-MICE}, also using the resampling technique to estimate the gradient norm.
   As for the tolerance on the statistical error and the resampling parameters, we used the same setup as for \texttt{SGD-MICE}.
   For \texttt{SGD}, we use a decreasing step-size of $\eta_k = 1/(L+Lk/50)$, tuned to achieve the best performance possible.
   We present the optimality gap versus iterations (top) and the number of gradient evaluations (bottom) for the three versions of \texttt{SGD}.
   In the top figure, it can be observed that \texttt{SGD-MICE} is able to maintain stability in the asymptotic regime, whereas \texttt{SGD-A} oscillates.
   Moreover, \texttt{SGD-MICE} achieves better results for the same cost in terms of gradient evaluations than \texttt{SGD-A} during the optimization process.
   Thus, we can conclude that \texttt{MICE} not only improves the performance of stochastic gradient optimizers in terms of gradient evaluations as it also improves their stability in the asymptotic regime.
   }
   \label{fig:quad_mice_mc}
 \end{figure}

\FloatBarrier
\subsection{Stochastic Rosenbrock function}\label{ex:2}

The goal of this example is to test the performance of {\texttt{Adam-MICE}, that is, \texttt{Adam} coupled with our gradient estimator \texttt{MICE}}, in minimizing the expected value of the stochastic Rosenbrock function in \eqref{eq:rosenbrock}, showing that \texttt{MICE} can be coupled with {different} first-order optimization methods in a non-intrusive manner.
Here we adapt the deterministic Rosenbrock function to the stochastic setting, specializing our optimization problem \eqref{eq:underlying.problem} with
\begin{equation}
 f(\bs{\xi},\bs{\theta}) = \left(a - \xi_0 + \theta_{0}\right)^{2} + b \left(- \xi_0^{2} + \xi_1 + \theta_{0}^{2} - \theta_{1}^2\right)^{2},
\end{equation}
where $a=1$, $b=100$, $\theta_0,\theta_1\sim\cl{N}(0,\sigma_\theta^2)$.
The objective function to be minimized is thus
\begin{equation}\label{eq:rosenbrock}
F(\bs{\xi}) = (a - \xi_0)^2  + \sigma_\theta^{2} + b\left(4 \sigma_\theta^{4} + (\xi_1 - \xi_0^{2})^2\right),
\end{equation}
{and its} gradient {is given by}
\begin{equation}
 \nabla_{\bs{\xi}}F(\bs{\xi})=
 \begin{bmatrix}
  - 2 a + 4 b \xi_0^{3} - 4 b \xi_0 \xi_1 + 2 \xi_0 \\
  - 2 b \xi_0^{2} + 2 b \xi_1
 \end{bmatrix},
\end{equation}
which coincides with the gradient of the deterministic Rosenbrock function.
Therefore, the optimal point of the stochastic Rosenbrock is the same as the one of the deterministic: $\bs{\xi}^*=(a,a^2)$.
To perform the optimization, we sample the stochastic gradient
\begin{equation}
 \nabla_{\bs{\xi}}f(\bs{\xi}, \bs{\theta})=
 \begin{bmatrix}
  - 2 a + 4 b \xi_0\left(\xi_0^{2} - \xi_1 - \theta_0^2 + \theta_1^2 \right) + 2 \xi_0 - 2 \theta_0 \\
  2 b \left(-\xi_0^{2} + \xi_1 + \theta_0^2 - \theta_1^2 \right)
 \end{bmatrix}.
\end{equation}

Although this is still a low dimensional example, minimizing the Rosenbrock function poses a difficult optimization problem for first-order methods; these tend to advance slowly in the region where the gradient has near-zero norm.
Moreover, when noise is introduced in gradient estimates, their relative error can become large, affecting the optimization convergence.

We compare the convergence of the {classical} \texttt{Adam} algorithm against \texttt{Adam-MICE}.
To illustrate the effect of the dispersion of the random variable $\bs{\theta}$, two distinct standard deviations are analyzed, namely $\sigma_\theta=10^{-4}$ and $\sigma_\theta=10^{-1}$.
As for the optimization setup, we set \texttt{Adam-MICE} with fixed step-size $0.2$ and \texttt{Adam} with a decreasing step-size $\eta_k = \nicefrac{0.02}{\sqrt{k}}$, which we observed to be the best step-sizes for each method.
The stopping criterion for both algorithms is set as $10^7$ gradient evaluations.
For \texttt{Adam-MICE}, we use $\epsilon = 0.7$, whereas for \texttt{Adam} we use a fixed batch size of $100$.

In Figures \ref{fig:rosen_1} and \ref{fig:rosen_2}, we present, for $\sigma_\theta$ of $10^{-1}$ and $10^{-4}$, respectively, the optimality gap for both \texttt{Adam} and \texttt{Adam-MICE} versus {the} number of gradients, iterations, and runtime in seconds.
It is clear that \texttt{Adam-MICE} is more stable than \texttt{Adam}, the latter oscillating as it approximates the optimal point in both cases.
The efficient control of the error in gradient estimates allows \texttt{Adam-MICE} to converge monotonically in the asymptotic phase.
Moreover, the number of iterations and the runtime are much smaller for \texttt{Adam-MICE} than for \texttt{Adam}.

\begin{figure}
 \includegraphics[width=.7\linewidth]{optimality_gap_Rosenbrock_1}
 \caption{Single run, stochastic Rosenbrock function example, Equation \eqref{eq:rosenbrock} with $\sigma_{\theta}=10^{-4}$.
 Optimality gap for \texttt{Adam} and \texttt{Adam-MICE} versus number of gradient evaluations (top), iterations (center), and runtime in seconds (bottom).}
 \label{fig:rosen_1}
\end{figure}

\begin{figure}
 \includegraphics[width=.7\linewidth]{optimality_gap_Rosenbrock_2}
 \caption{Single run, stochastic Rosenbrock function example, Equation \eqref{eq:rosenbrock} with $\sigma_{\theta}=10^{-1}$.
 Optimality gap for \texttt{Adam} and \texttt{Adam-MICE} versus number of gradient evaluations (top), iterations (center), and runtime in seconds (bottom).}
 \label{fig:rosen_2}
\end{figure}

\FloatBarrier

\subsection{Logistic regression} \label{ex:3}
In this example, we train logistic regression models using \texttt{SGD-MICE}, \texttt{SAG} \cite{schmidt2017minimizing}, \texttt{SAGA} \cite{defazio2014saga}, \texttt{SARAH} \cite{nguyen2017sarah}, and \texttt{SVRG} \cite{Joh13} to compare their performances.
Here, we present a more practical application of \texttt{MICE}, where we can test its performance on high-dimensional settings with finite populations.
Therefore, we calculate the statistical error as in Remark \ref{rem:mice.optimal.setting.discrete} and use \eqref{eq:optimal.sampling.2} to obtain the optimal sample sizes.
To train the logistic regression model for binary classification, we use the $\ell_2$-regularized log-loss function
\begin{equation}\label{eq:logistic}
 F(\bs{\xi}) = \frac{1}{N} \sum_{i=1}^N f\left(\bs{\xi}, \bs{\theta}_i= (\bs{x}_i, y_i)\right) = \frac{1}{N} \sum_{i=1}^N \log(1 + \exp(-y_i \bs{\xi} \cdot \bs{x}_i))
 + \frac{\lambda}{2} \|\bs{\xi}\|^2,
\end{equation}
where each data point $(\bs{x}_i, y_i)$ is such that $\bs{x}_i \in \bb{R}^{d_{\bs{\xi}}}$ and $y_i \in \{-1,1\}$.
We use the datasets \emph{mushrooms}, \emph{gisette}, and \emph{Higgs}, obtained from LibSVM\footnote{
 \url{https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary.html}
}.
The size of the datasets $N$, number of features $d_{\bs{\xi}}$, and regularization parameters $\lambda$ are presented in Table \ref{tab:datasets}.
\begin{table}[h]
 \caption{Size, number of features, and regularization parameters for the datasets used in the logistic regression example.}
 \label{tab:datasets}
 \begin{tabular}{cccc}
  \toprule
  Dataset          & Size     & Number of features & $\lambda$ \\ \midrule
  \emph{mushrooms} & 8124     & 112                & $10^{-5}$ \\
  \emph{gisette}   & 6000     & 5000               & $10^{-4}$ \\
  \emph{HIGGS}     & 11000000 & 28                 & $10^{-4}$ \\ \bottomrule
 \end{tabular}
\end{table}

When using \texttt{SGD-MICE} for training the logistic regression model, we use $\epsilon = 0.5$.
For the other methods, we use batchsizes of size $10$.
Since we have finite populations, we use \eqref{eq:optimal.sampling.2} to calculate the sample-sizes.
\texttt{SGD-MICE} step is based on the Lipschitz smoothness and strong-convexity constant of the true objective function as presented in Proposition \ref{prp:convergence.mice}. {Conversely}, the other methods rely on a Lipschitz constant that must hold a.s., which we refer to as $L_\text{a.s.}$.
The step-sizes for \texttt{SAG}, \texttt{SAGA}, \texttt{SARAH}, and \texttt{SVRG} are presented in Table \ref{tab:steps}.
These steps were chosen as the best performing for each case based on the recommendations of their original papers.
\begin{table}[h]
 \caption{Step-sizes chosen for each method for the logistic regression example.}
 \label{tab:steps}
 \begin{tabular}{lcccc}
  \toprule
  Method    & \texttt{SAG}                          & \texttt{SAGA}                        & \texttt{SARAH}              & \texttt{SVRG}               \\ \midrule
  Step-size & $\frac{1}{16(L_\text{a.s.} + \mu N)}$ & $\frac{1}{2(L_\text{a.s.} + \mu N)}$ & $\frac{1}{2 L_\text{a.s.}}$ & $\frac{1}{2 L_\text{a.s.}}$ \\ \bottomrule
 \end{tabular}
\end{table}

Figures \ref{fig:logreg.mushrooms}, \ref{fig:logreg.gisette}, and \ref{fig:logreg.higgs} present the relative optimality gap (the optimality gap normalized by its starting value) for, respectively, the \emph{mushrooms}, \emph{gisette}, and \emph{HIGGS} datasets versus {the} number of gradient evaluations, iterations and runtime in seconds.
In the \emph{mushrooms} dataset, \texttt{SGD-MICE} decreases the optimality gap more than the other methods for the same number of gradient samples during the whole optimization process.
Moreover, the total number of iterations is much smaller than for the other methods.
However, the overhead of our \texttt{SGD-MICE} implementation becomes plain to see when we compare the runtimes; the current implementation of \texttt{SGD-MICE} is more costly than the other methods.
In the \emph{gisette} dataset, the overall behavior of the methods is the same as in the first data set.
Finally, for the \emph{HIGGS} dataset, which is a much larger dataset, \texttt{SGD-MICE} performs better than the other methods both in the number of gradient evaluations and in runtime.
Moreover, \texttt{SGD-MICE} performed around $2 \times 10^3$ iterations while the other methods required more than $10^6$ iterations.
\begin{figure}[h]
  \includegraphics[width=.7\linewidth]{relative_loss_gap_mushrooms}
  \caption{{Single run}, logistic regression example, Equation \eqref{eq:logistic} for the \emph{mushrooms} dataset. Relative optimality gap versus number of gradient evaluations (top), iterations (center), and runtime in seconds (bottom) for \texttt{SGD-MICE}, \texttt{SAG}, \texttt{SAGA}, \texttt{SARAH}, and \texttt{SVRG}.}
  \label{fig:logreg.mushrooms}
\end{figure}
\begin{figure}[h]
  \includegraphics[width=.7\linewidth]{relative_loss_gap_gisette}
  \caption{{Single run}, logistic regression example, Equation \eqref{eq:logistic} for the \emph{gisette} dataset. Relative optimality gap versus number of gradient evaluations (top), iterations (center), and runtime in seconds (bottom) for \texttt{SGD-MICE}, \texttt{SAG}, \texttt{SAGA}, \texttt{SARAH}, and \texttt{SVRG}.}
  \label{fig:logreg.gisette}
\end{figure}
\begin{figure}[h]
  \includegraphics[width=.7\linewidth]{relative_loss_gap_HIGGS}
  \caption{{Single run}, logistic regression example, Equation \eqref{eq:logistic} for the \emph{HIGGS} dataset. Relative optimality gap versus number of gradient evaluations (top), iterations (center), and runtime in seconds (bottom) for \texttt{SGD-MICE}, \texttt{SAG}, \texttt{SAGA}, \texttt{SARAH}, and \texttt{SVRG}.}
  \label{fig:logreg.higgs}
\end{figure}

From the results of this example, we observe that \texttt{MICE} performs well in problems with a reasonably large number of parameters, for instance, $5000$ in the \emph{gisette}, and finite dataset populations ranging from the thousands to the millions.
Moreover, it seems that \texttt{SGD-MICE}'s performance compared to the other methods increases as the population size grows.
Note that both \texttt{SAG} and \texttt{SAGA} need to decrease their step-sizes as the sample-size increases, and that \texttt{SARAH} and \texttt{SVRG} need to reevaluate the full-gradient after a few epochs to keep their convergence.