[name=Corollary]corollary [name=Exercise, style=definition]exercise
[name=Lemma]lemma [name=Proposition]proposition [name=Remark,
style=remark]remark [name=Theorem]theorem [name=Definition]definition
[name=Example]example [name=Question]question [name=Assumption]assumption
[name=Claim]claim

# Reproducibility in Optimization:
Theoretical Framework and Limits

Kwangjun Ahn
MIT EECS
<EMAIL_ADDRESS>&Prateek Jain
Google Research
<EMAIL_ADDRESS>&Ziwei Ji
Google Research
<EMAIL_ADDRESS>&Satyen Kale
Google Research
<EMAIL_ADDRESS>&Praneeth Netrapalli
Google Research
<EMAIL_ADDRESS>&Gil I. Shamir
Google Research, Brain Team
<EMAIL_ADDRESS>Part of this work was done when Kwangjun Ahn was an intern
at Google Research.

###### Abstract

We initiate a formal study of reproducibility in optimization. We define a
quantitative measure of reproducibility of optimization procedures in the face
of noisy or error-prone operations such as inexact or stochastic gradient
computations or inexact initialization. We then analyze several convex
optimization settings of interest such as smooth, non-smooth, and strongly-
convex objective functions and establish tight bounds on the limits of
reproducibility in each setting. Our analysis reveals a fundamental trade-off
between computation and reproducibility: more computation is necessary (and
sufficient) for better reproducibility.

###  1 Introduction

Machine learned models are increasingly entering wider ranges of domains in
our lives, driving a constantly increasing number of important systems. Large
scale systems can be trained in highly parallel and distributed training
environments, with a large amount of randomness in training the models. While
some systems may tolerate such randomness leading to models that differ from
one another every time a model retrains, for many applications, reproducible
models are required, where slight changes in training do not lead to drastic
differences in the model learned.

Beyond practical deployments of machine learned models, the reproducibility
crisis in the machine learning academic world has also been well-documented:
see (Pineau et al., 2021) and the references therein for an excellent
discussion of the reasons for irreproducibility (insufficient exploration of
hyperparameters and experimental setups, lack of sufficient documentation,
inaccessible code, and different computational hardware) and for mitigation
recommendations. Recent papers (Chen et al., 2020; D’Amour et al., 2020;
Dusenberry et al., 2020; Snapp and Shamir, 2021; Summers and Dinneen, 2021; Yu
et al., 2021) have also demonstrated that even when models are trained on
identical datasets with identical optimization algorithms, architectures, and
hyperparameters, they can produce significantly different predictions on the
same example. This type of irreproducibility may be caused by multiple factors
(D’Amour et al., 2020; Fort et al., 2020; Frankle et al., 2020; Shallue et
al., 2018; Snapp and Shamir, 2021; Summers and Dinneen, 2021), such as non-
convexity of the objective, random initialization, nondeterminism in training
such as data shuffling, parallelism, random schedules, hardware used, and
round off quantization errors. Perhaps surprisingly, even if we control for
the randomness by using the same “seed" for model initialization, other
factors such as numerical errors introduced due to nondeterminism of modern
GPUs (see, e.g., (Zhuang et al., 2021)) may still lead to significant
differences. It was empirically shown (see, e.g., Achille et al. (2017)) that
slight deviations early in training can lead to different optima, with
substantial differences in resulting models. Thus we are forced to accept some
fundamental level of irreproducibility that persists even after fixing all
aspects of the training process under our control.

The goal of this paper is to initiate a formal study of the fundamental limits
of irreproducibility. In particular, we focus on the most basic training
process: convex optimization. At first glance, it might seem surprising that
convex optimization procedures can exhibit irreproducibility since they’re
guaranteed to converge to an optimal solution. However, in practice, the
default convex optimization algorithms are iterative first-order methods;
methods that only use a first-order oracle to provide an approximate gradient
of the function at a given point, and converge to an approximately optimal
solution. The first-order oracle is a source of irreproducibility. In
stochastic gradient descent, it returns a random vector whose expectation is
the true gradient. The randomness in the stochastic gradients can lead to
different outcomes of the optimization process. Similarly, there are numerical
errors that can arise in the computation of the gradients due to inherent
nondetermism in modern GPUs. Beyond the first-order oracle, irreproducibility
may also arise in convex optimization procedures because the initial point is
chosen randomly. Thus, we attempt to answer the following questions for convex
optimization procedures operating with the above sources of irreproducibility:

  * $\blacksquare$

What are the fundamental limits of reproducibility for any convex optimization
procedure?

  * $\blacksquare$

Can we design practical and efficient first-order methods that achieve these
limits?

We study these questions in a variety of settings; including general non-
smooth convex functions, smooth convex functions, strongly-convex functions,
finite-sum functions, and stochastic convex optimization, under the different
sources of irreproducibility mentioned above. To the best of our knowledge, no
prior theoretical work considered such questions.

The primary contribution of this paper is conceptual: the development of a
rigorous theoretical framework to study the fundamental limits of
reproducibility in convex optimization. The concepts developed in this
framework can be extended easily to other settings of interest such as non-
convex optimization. The technical contribution of this paper is the
development of lower bounds on the amount of reproducibility, and matching
upper bounds via analysis of specific first-order algorithms in all the
different settings of convex optimization described above. Detailed technical
descriptions of the results appear below in Subsection 1.1.

At a high level, our study provides the same message for all the different
optimization settings we consider. On the lower bound side, we find that any
first-order method would need to trade-off convergence rate (computational
complexity) for more reproduciblity. On the upper bound side, we find that
various forms of gradient descent, when run with lower step-size (and
correspondingly, more iterations) already achieve the fundamental limits of
reproducibility. One (somewhat surprising) consequence of our results is that
advanced techniques like regularization, variance reduction, and acceleration
do not improve reproducibility over standard gradient descent methods.

We show, for example, that when optimizing a Lipschitz non-smooth convex
function $f$ on a bounded domain using a first-order oracle that computes
gradients with even vanishingly small error, two runs of any first-order
method that obtains an $\varepsilon$ suboptimal solution of $f$ after $T$
iterations can generate solutions that are
$\Omega(\nicefrac{{1}}{{(\sqrt{T}\varepsilon)}})$ apart in $\ell_{2}$
distance. Thus, if we run the method for the standard
$T=O(\nicefrac{{1}}{{\varepsilon^{2}}})$ iterations required to obtain
$\varepsilon$-approximate solution for a non-smooth function, then obtained
solutions can deviate by $\Omega(1)$ distance; i.e., the method is _maximally_
irreproducible. To ensure that irreproducibility is small, i.e. that the
solutions are within a small distance $\gamma$ of each other, we will need to
run at least $\Omega(\nicefrac{{1}}{{(\varepsilon^{2}\gamma^{2})}})$
iterations. Interestingly, standard gradient descent with appropriately chosen
learning rate and number of iterations already achieves this trade-off.

Our results demonstrate the challenge of reproducibility even for standard
convex optimization. While we provide matching lower and upper bounds in
certain general settings, in Section 7, we outline several important open
directions. Solutions to these problems should enable better understanding of
reproducibility even for deep learning.

####  1.1 Summary of results

Table 1: Summary of $(\varepsilon,\delta)$-deviation bounds for various convex optimization settings.  | Stochastic Inexact | Non-stochastic Inexact | Inexact Initialization
---|---|---|---
Gradient Oracle | Gradient Oracle | Oracle
Section 3 | Section 4 | Section 5
Smooth | $\Theta(\nicefrac{{\delta^{2}}}{{(T\varepsilon^{2})}})$ | $\Theta(\nicefrac{{\delta^{2}}}{{\varepsilon^{2}}})$ | $\Theta(\delta^{2})$
Smooth Strongly-Cvx. | $\Theta(\nicefrac{{\delta^{2}}}{{T}}\wedge\varepsilon)$ | $\Theta(\delta^{2}\wedge\varepsilon)$ | $\Theta(e^{-\Omega(T)}\delta^{2}\wedge\varepsilon)$
Nonsmooth | $\Theta(\nicefrac{{1}}{{(T\varepsilon^{2})}})$ | $\Theta(\nicefrac{{1}}{{(T\varepsilon^{2})}}+\nicefrac{{\delta^{2}}}{{\varepsilon^{2}}})$ | $\Theta(\nicefrac{{1}}{{(T\varepsilon^{2})}}+\delta^{2})$
Nonsmooth Strongly-Cvx. | $\Theta(\nicefrac{{1}}{{T}}\wedge{\varepsilon})$ | $\Theta((\nicefrac{{1}}{{T}}+{\delta^{2}})\wedge{\varepsilon})$ | $\Theta(\nicefrac{{1}}{{T}}\wedge{\varepsilon})$

Table 1 summarizes our key results for our measure of irreproducibility,
$(\varepsilon,\delta)$-deviation (see Section 2). The
$(\varepsilon,\delta)$-deviation measures the amount of change between the
outputs of two independent runs of an optimization algorithm, that is
guaranteed to achieve $\varepsilon$-suboptimality after $T$ iterations, when
the computations of the algorithm incur errors of magnitude up to $\delta$. We
specifically focus on three different sources of errors: i) stochastic
gradient oracles, ii) gradient oracles with non-deterministic numerical errors
(Section 2), and iii) inexact initialization for the optimizer (Section 2). We
analyze the deviation under these sources of errors for four types of function
classes: smooth convex functions, non-smooth but Lipschitz convex functions
and strongly-convex restrictions of the two. Throughout the paper, $a\wedge b$
denotes the minimum between $a$ and $b$.

All lower bounds are for first-order iterative algorithms (_à la_ Nesterov
(2018)) that we formally define in (FOI). This is a large class of iterative
optimization methods, including Stochastic Gradient Descent (SGD), which
construct successive iterates adaptively in the linear span of previous
iterates. Additionally, for smooth costs and stochastic inexact gradient
oracle, we have an information theoretic lower bound of
$\Omega(\nicefrac{{\delta^{2}}}{{(T\varepsilon^{2})}})$ when
$\varepsilon\lesssim\delta^{2}$ (Subsection B.1). We believe such
informtation-theoretic lower bounds can be shown for all the settings in this
paper. As for the upper bounds, they are all obtained using slowed-down SGD:
i.e. SGD using smaller learning rates and more iterations.

For the non-strongly convex cases, one may expect to have high
irreproducibility if the minima form a large flat region; however,
surprisingly, our upper bounds show that we can always bound the extent of
irreproducibility via slowed-down SGD. In the non-smooth cases, the main
observation is that the deviation does not depend on scale of perturbation by
the gradient oracle, i.e., any $\delta>0$ can lead to fairly irreproducible
solutions. The non-stochastic gradient oracle setting is strictly harder than
the stochastic setting. Naturally, the lower and upper bounds on
reproducibility are worse. Interestingly, even though strong convexity implies
uniqueness of the global optimum, which intuitively should lead to highly
reproducible solutions, we show that when faced with sources of error in
computations, the deviation can still be significantly large for any
algorithm.

Finally, we study reproducibility of optimization in machine learning
settings. Here we have additional structure such as finite-sum minimization
(for optimizing training loss) and stochastic convex optimization (for
optimizing population loss). We define appropriate notions of errors for these
problems and analyze two settings of particular interest. Our main results
(Subsection 6.1 and Subsection 6.2) show that despite the additional structure
in these problems, the bounds given by Table 1 for the specific settings are
nonetheless tight. One consequence is that more sophisticated techniques for
these problems such as variance reduction don’t improve reproducibility.

####  1.2 Related work

##### Related notions.

In the scientific world, the terms reproducibility and replicability are often
used interchangeably, but here we distinguish the two, following Pineau et al.
(2021). Reproducibility refers to the requirement that results obtained by a
computational procedure (e.g. an experiment or a statistical analysis of a
data set) should be the same (or largely similar) when the procedure is
repeated using the same code on the same data, whereas replicability is a
different notion that requires that results be reliably the same or similar
when the data are changed. The field of statistical hypothesis testing
(Lehmann and Romano, 2005) provides rigorous and principled techniques to
minimize false discoveries and thereby promote replicability. The notion of
algorithmic stability can also be seen as quantifying the amount of change in
the output when a single data sample is changed. This notion has been
extensively studied in the context of providing algorithm-dependent
generalization bounds (Bousquet and Elisseeff, 2001; Kutin and Niyogi, 2002)
and in developing differentially private algorithms (Dwork et al., 2006;
McSherry and Talwar, 2007). In very recent concurrent work, Impagliazzo et al.
(2022) define a notion of replicability in learning that is quite different
from ours: they aim to develop algorithms that generate the exact same output
with reasonable probability given a fresh sample. Note that despite the title
of their paper, technically the notion studied is replicability, not
reproducibility, since they study the output of algorithms when the input data
are changed.

In this paper, we specifically focus on reproducibility: how much can the
results of a computation differ when it is re-run on the same data with the
same code? Hence, both hypothesis testing and algorithmic stability are
orthogonal to the study in this paper, although some of our upper bounds use
similar analysis techniques as algorithmic stability. On a different note,
similar to replicability, the boundary between the notions of
irreproducibility and uncertainty in deep models is rather blurred. Several
papers considered different aspects of uncertainty (see, e.g.,
(Lakshminarayanan et al., 2017) and references therein), but this line of work
has been empirical in nature.

##### Inexact oracles in optimization.

The optimization community has studied the consequences of using inexact or
error-prone gradient oracles in optimization. Several papers (e.g. (Aybat et
al., 2020; Devolder et al., 2014; d’Aspremont, 2008; Cohen et al., 2018) have
developed bounds on the optimization error incurred due to the use of inexact
oracles. While the sources of errors are similar to the ones studied in this
paper, the quantities of interest in these papers are convergence rate and
optimization error rather than reproducibility. Interestingly, despite the
different objective, some of the high-level conclusions are similar to our
paper: for example, accelerated gradient methods do not outperform standard
classical methods when used with inexact gradient oracles.

##### Techniques to improve reproducibility in practice.

Several recent empirical papers considered methods that can reduce levels of
irreproducibility in deep models despite nondeterminism in training. Smooth
activations (Du, 2019; Mhaskar, 1997) have been shown (Shamir et al., 2020) to
improve reproducibilty over popular activations, as the Rectified Linear Unit
(Nair and Hinton, 2010). Ensembles (Dietterich, 2000) leverage diversity of
multiple different solutions to produce an average more reproducible one
(Allen-Zhu and Li, 2020). Co-distillation (Anil et al., 2018) and Anti-
distillation (Shamir and Coviello, 2020) leverage ensembles to further push
deployed models to be more reproducible. Imposing constraints (Bhojanapalli et
al., 2021; Shamir, 2018) forces models to prefer some solutions over others,
but may come at the cost of reducing model accuracy performance.

##### Robustness of dynamical systems.

The upper bound results in our paper can be interpreted as robustness results
of the (sub)gradient descent dynamics against disturbances. In particular, our
upper bounds can be viewed as some variants of the input-to-state stability
(Sontag and Wang, 1995) results for the dynamics (see, e.g., (Tu et al., 2022,
Definition 3.2)).

###  2 Problem Formulation

In this section, we define a quantitative measure of irreproducibility
amenable to a theoretical analysis. Intuitively, a computation is reproducible
if it generates the exact same output given the same inputs on two different
runs. Irreproducibility arises because low-level operations of a computation
produce different answers on two runs due to either randomness or non-
determinism.

Our computation of interest is convex optimization via first-order methods,
where initialization and gradient computations are the primary operations that
constitute the computation and are subject to errors leading to inexact
outputs. A natural measure of irreproducibility is the amount of change in the
computed solution to the convex optimization problem under inexact gradient
computations or inexact initialization. However, there are two nuances that
must be carefully handled here. First, a trivial procedure which ignores its
input and outputs a constant solution is perfectly reproducible!
Unfortunately, it is perfectly useless as a convex optimization procedure as
well. Thus, in order to compare different procedures by their reproducibility
metrics, we must assume that the procedures are guaranteed to converge to an
optimal solution. The second nuance is that we need to assume that the errors
in the gradient or initialization computations are bounded in some manner.
Evidently, without such an assumption, any non-trivial convex optimization
procedure will be extremely irreproducible. We now use the above
considerations to develop a precise definition of a measure of
irreproducibility.

##### Convex function classes.

We assume that the function to be optimized is chosen from a certain class,
$\mathcal{F}$, of convex functions, along with their domains, satisfying
suitable regularity conditions (e.g. Lipschitzness, smoothness, strong-
convexity, etc.) to develop convergence rates. For clarity, we will suppress
exact dependence on smoothness, Lipschitzness, and strong-convexity
parameters. In particular, “smooth” will denote a convex function whose
gradients are $O(1)$-Lipschitz continuous, “non-smooth” a convex function that
is $O(1)$-Lipschitz continuous, and “strongly-convex” an $\Omega(1)$-strongly-
convex function. Here, $O(1)$ and $\Omega(1)$ denote universal constants
independent of the dimension or other problem dependent quantities, which we
leave unspecified to ease the exposition.

##### Convex optimization procedures.

A first-order convex optimization procedure for $\mathcal{F}$ is an algorithm
that, given any function $f\in\mathcal{F}$, and access to two (potentially
noisy) oracles – an initialization oracle, which generates the initial point,
and a gradient oracle, which computes gradients for $f$ at any given query
point – generates a candidate solution ${\bm{x}}_{\mathrm{out}}$ for the
problem of minimizing $f$ over its domain. Note that the algorithm can only
access $f$ via the oracles provided. We call such an algorithm
$\varepsilon$-accurate if it guarantees that
$\operatorname*{\mathbb{E}}f({\bm{x}}_{\mathrm{out}})-\inf_{{\bm{x}}\in\text{dom}f}f({\bm{x}})\leq\varepsilon$,
where the expectation is over any randomness in the computation of
${\bm{x}}_{\mathrm{out}}$. Several of our lower bounds require more structure
for the algorithm: specifically, a first-order iterative (FOI) algorithm (_à
la_ Nesterov (2018)) is one that starting from the point ${\bm{x}}_{0}$
generated by the initialization oracle, constructs successive iterates

$\displaystyle{\bm{x}}_{t}={\bm{x}}_{0}-\textstyle\sum_{i=0}^{t-1}\lambda^{(t)}_{i}g({\bm{x}}_{i})\quad\text{for
some }\lambda_{i}^{(t)},~{}i=0,\dots,t-1,$ (FOI)

where $g({\bm{x}}_{i})$ is the output of the gradient oracle query at
${\bm{x}}_{i}$, and outputs ${\bm{x}}_{T}$ for some integer $T>0$. We
emphasize that for all $t$, the coefficients $\lambda^{(t)}_{i}$ can be chosen
adaptively based on all the previous computations. The above class of
algorithms is a canonical one to consider when proving lower bounds against
gradient oracle based optimization algorithms. We refer readers to (Nesterov,
2018, §2.1.2) for more background. For the case of nonsmooth costs, we
additionally assume that the coefficient of the latest gradient is nonzero,
i.e., $\lambda^{(t)}_{t-1}\neq 0$ for all $t$. We also note that one of our
lower bound results (Subsection B.1) is _information-theoretic_ (in the sense
of Nemirovski and Yudin (1983)).

##### Sources of errors in computation.

Errors arise due to inexactness in the outputs of the initialization or
gradient oracles. Queries to these oracles on two different runs of the same
algorithm might yield different outputs, but we will control the errors by
assuming that the outputs are close to some reference point (that remains
fixed over different runs) in a suitable metric. [$\delta$-bounded inexact
initialization oracle] Given a function $f\in\mathcal{F}$ and a reference
initialization point ${\bm{x}}_{0}^{\text{ref}}\in\text{dom}f$, a
$\delta$-bounded inexact initialization oracle for $f$ is one that generates
an initial point ${\bm{x}}_{0}\in\text{dom}f$ such that
$\|{\bm{x}}_{0}-{\bm{x}}_{0}^{\text{ref}}\|\leq\delta$. The gradient
computation oracle is said to be $\delta$-bounded if for any $f\in\mathcal{F}$
and any point ${\bm{x}}\in\text{dom}f$, it outputs a vector $g({\bm{x}})$ such
that $\operatorname*{\mathbb{E}}\left\lVert g({\bm{x}})-\nabla
f({\bm{x}})\right\rVert^{2}\leq\delta^{2}$ for some $\nabla
f({\bm{x}})\in\partial f({\bm{x}})$, where the expectation is over any
randomness in the computation of $g({\bm{x}})$. We consider both stochastic
and non-stochastic inexact $\delta$-bounded gradient oracles. A stochastic
gradient oracle has the additional property that its output $g({\bm{x}})$ is a
random vector such that $\operatorname*{\mathbb{E}}g({\bm{x}})=\nabla
f({\bm{x}})$, with different queries being independent of each other.
Stochastic inexact gradient oracles arise naturally in machine learning
applications due to randomness in minibatching. Non-stochastic inexact
gradient oracles model non-deterministic numerical errors due to the
accumulation of floating point errors; for giant machine learning models with
billions of parameters, individual floating point errors could add up to a
noticeable large error. We formally define the two types of inexact oracles
below. [$\delta$-bounded inexact gradient oracle] Given a function
$f\in\mathcal{F}$, and ${\bm{x}}\in\text{dom}f$, let $\partial f({\bm{x}})$
denote the sub-differential of $f$ at ${\bm{x}}$.

  * (a)

A stochastic inexact $\delta$-bounded gradient oracle outputs a random vector
$g({\bm{x}})$ such that $\operatorname*{\mathbb{E}}{g({\bm{x}})}=\nabla
f({\bm{x}})$ and $\operatorname*{\mathbb{E}}\left\lVert g({\bm{x}})-\nabla
f({\bm{x}})\right\rVert^{2}\leq\delta^{2}$ for some $\nabla
f({\bm{x}})\in\partial f({\bm{x}})$. The expectation is over the randomness in
$g({\bm{x}})$ which is also assumed to be independent for each oracle call.

  * (b)

A non-stochastic inexact $\delta$-bounded gradient oracle outputs a non-
deterministic vector $g({\bm{x}})$ such that $\left\lVert g({\bm{x}})-\nabla
f({\bm{x}})\right\rVert^{2}\leq\delta^{2}$ for some $\nabla
f({\bm{x}})\in\partial f({\bm{x}})$.

##### Measure of irreproducibility.

Let $\mathcal{A}$ be a first-order, $\varepsilon$-accurate convex optimization
procedure for $\mathcal{F}$ with access to either a $\delta$-bounded
initialization oracle (and an exact gradient oracle), or a $\delta$-bounded
gradient oracle (and an exact initialization oracle). The
$(\varepsilon,\delta)$-deviation is [$(\varepsilon,\delta)$-deviation] Given a
function $f\in\mathcal{F}$, let ${\bm{x}}_{f}$ and ${\bm{x}}^{\prime}_{f}$
denote the outputs of $\mathcal{A}$ on two independent runs of $\mathcal{A}$
that result in $\varepsilon$-accurate solutions.

  * (a)

If a stochastic inexact $\delta$-bounded gradient oracle is used, the
$(\varepsilon,\delta)$-deviation of $\mathcal{A}$ is defined as
$\sup_{f\in\mathcal{F}}\operatorname*{\mathbb{E}}\|{\bm{x}}_{f}-{\bm{x}}^{\prime}_{f}\|^{2}$
where the randomness is over the stochastic oracle in the two runs.

  * (b)

If either an inexact $\delta$-bounded initialization oracle or a
$\delta$-bounded non-stochastic inexact gradient oracle is used, the
$(\varepsilon,\delta)$-deviation of $\mathcal{A}$ is defined as
$\sup_{f\in\mathcal{F}}\sup\|{\bm{x}}_{f}-{\bm{x}}^{\prime}_{f}\|^{2}$, where
the inner supremum is over the two runs.

Here we note that ${\bm{x}}_{f}$ is not necessarily the last iterate of a
first-order algorithm. Depending on cases, sometimes we choose ${\bm{x}}_{f}$
to be the (weighted) average iterate. We consider
$\|x_{f}-x^{\prime}_{f}\|^{2}$ for the notion of deviation instead of
$\|x_{f}-x^{*}\|^{2}$; note that $x^{*}$ may not even be unique without strong
convexity.

##### Other notions?

Alternate definitions are certainly plausible. For example, in machine
learning applications, the computed solution ${\bm{x}}$ may be used as a
parameter vector for a predictor function $g_{{\bm{x}}}$ which maps inputs $z$
to real-valued predictions $\hat{y}$ whose quality is measured by a loss
function $\ell(\hat{y},y)$ where $y$ is the true label of $z$. Let ${\bm{x}}$
and ${\bm{x}}^{\prime}$ be outputs of two different runs of the optimization
algorithm. Then one can also define $(\varepsilon,\delta)$-deviation based on
_prediction reproducibility_ : e.g.
$\sup_{z}|g_{{\bm{x}}}(z)-g_{{\bm{x}}^{\prime}}(z)|$ or
$\operatorname*{\mathbb{E}}_{z}|g_{{\bm{x}}}(z)-g_{{\bm{x}}^{\prime}}(z)|$, or
_loss reproducibility_ : e.g.
$\sup_{(z,y)}|\ell(g_{{\bm{x}}}(z),y)-\ell(g_{{\bm{x}}^{\prime}}(z),y)|$ or
$\operatorname*{\mathbb{E}}_{(z,y)}|\ell(g_{{\bm{x}}}(z),y)-\ell(g_{{\bm{x}}^{\prime}}(z),y)|$,
where the expectation is over the distribution of the examples. Nevertheless,
we adopt Section 2, which is based on _parameter reproducibility_ , in this
paper for the following reasons:

  * •

The convex optimization problems studied in this paper are more
basic/fundamental than more structured ML optimization. To the best of our
knowledge, there is no pre-existing theory even for this basic setting.
Parameter reproducibility is a more natural definition here since there is no
notion of prediction or loss.

  * •

In many ML applications, the predictor function $g_{{\bm{x}}}$ is Lipschitz in
${\bm{x}}$ for any input $z$. In such cases, $(\varepsilon,\delta)$-deviation
bounds for parameter reproducibility immediately transform into
$(\varepsilon,\delta)$-deviation bounds for prediction reproducibility.
Similar transformations are also generally possible
$(\varepsilon,\delta)$-deviation bounds for loss reproducibility.

  * •

Without knowledge of how a learned parametric function is deployed for making
predictions, it is difficult to analyze prediction reproducibility even in ML
settings. Hence, parameter deviations provide a reasonable first
approximation. Furthermore, in real systems, we often optimize multiple
metrics (e.g. performance on sub-segments of populations, for fairness).
Parameter reproducibility gives greater assurance on all metrics. A similar
argument applies when the test distribution is different from the training
distribution.

  * •

Finally, several recent works (Shamir et al., 2020; D’Amour et al., 2020) have
empirically observed that even if two different runs of the algorithm resulted
in parameters that have nearly the same loss, the predictions on test examples
could be very different. One reason this happens is because surrogate losses
are used in place of the true metric for optimization. So simply achieving
loss reproducibility may not be sufficient for practical applications.

###  3 Reproducibility with Stochastic Inexact Gradient Oracles

In this section we consider reproducibility of optimizing a convex function
$f\in\mathcal{F}$ where we can access $f$ only via a stochastic gradient
oracle (see Section 2). This setting covers several important ML optimization
scenarios, e.g., when the training data is randomly sampled from a population
or when the selection of mini-batches is randomized. Our main result is the
following theorem.

For any $\varepsilon,\delta>0$, and number of iterations $T$, the
$(\varepsilon,\delta)$-deviation for optimizing convex functions with a
stochastic inexact gradient oracle is as follows. Unless indicated otherwise,
the lower bounds hold for any FOI algorithm, and the upper bound is achieved
by stochastic gradient descent for
$T=\Omega(\nicefrac{{1}}{{\varepsilon^{2}}})$ in the non-strongly-convex
settings and $T=\Omega(\nicefrac{{1}}{{\varepsilon}})$ in the strongly-convex
settings.

  * $\blacksquare$

Smooth functions: $(\varepsilon,\delta)$-deviation is
$\Theta(\nicefrac{{\delta^{2}}}{{(T\varepsilon^{2})}})$. Furthermore, for the
case $\varepsilon\lesssim\delta^{2}$ there is a matching information theoretic
lower bound.

  * $\blacksquare$

Smooth and strongly convex functions: $(\varepsilon,\delta)$-deviation is
$\Theta(\nicefrac{{\delta^{2}}}{{T}}\wedge\varepsilon)$.

  * $\blacksquare$

Lipschitz (non-smooth) functions: $(\varepsilon,\delta)$-deviation is
$\Theta(\nicefrac{{1}}{{(T\varepsilon^{2})}})$.

  * $\blacksquare$

Lipschitz (non-smooth) and strongly convex functions:
$(\varepsilon,\delta)$-deviation is
$\Theta(\nicefrac{{1}}{{T}}\wedge\varepsilon)$.

This theorem is proved in several pieces, one corresponding to each setting
(i.e., smooth, smooth & strongly convex etc.) and type of bound (i.e., upper
or lower). The precise details are given in Appendix A. A few remarks are in
order.

Smoothness: Intuitively, smoothness would ensure that a slight error in
gradient computation need not imply catastrophic deviation in the iterates.
Our matching upper and lower bounds confirm this intuition.

Non-smoothness: Our results show that even a slight amount of noise in the
gradients can lead to drastic irreproducibility in the non-smooth case.
Intuitively, the reason for this phenomenon is that non-smooth functions have
non-differentiable points where the slightest amount of noise in the gradients
can lead to drastically different behavior. This is in line with the empirical
observations of Shamir et al. (2020).

Strong convexity: Note that the deviation is smaller than other cases due to
the existence of a unique minimizer; the $\varepsilon$-accuracy in the cost
already implies an $O(\varepsilon)$ upper bound on the deviation.

Furthermore, in all the settings, our results show that using an algorithm
with a larger number of iterations is more helpful; intuitively that is
because more gradient samples from the oracle can help reduce the sample noise
as an averaging. For example, for smooth functions, if a gradient descent type
method is run for the standard $T=O(\nicefrac{{1}}{{\varepsilon^{2}}})$
iterations, then the solution obtained might still suffer deviation of
$\delta^{2}$, which is independent of $\varepsilon$. Thus, in order to obtain
low deviation, we are forced to run the method for
$\omega(\nicefrac{{1}}{{\varepsilon^{2}}})$ iterations. We also note that by
relying on the standard convergence rate lower bounds (Nemirovski and Yudin,
1983), the requirement on $T$ e.g., $T=\Omega(1/\varepsilon^{2})$ for the non-
strongly convex setting, is without loss of generality. Finally, we remark
that the lower bound for smooth functions is information theoretic (i.e.,
holds against _any_ algorithm) for $\varepsilon\lesssim\delta^{2}$.

###  4 Reproducibility with Non-Stochastic Inexact Gradient Oracles

In this section, we study reproducibility for optimizing a function
$f\in\mathcal{F}$ with non-stochastic inexact gradient oracle access (Section
2). In particular, we establish lower and upper bounds for first-order
algorithms (FOI) on the $(\varepsilon,\delta)$-deviation (Section 2) despite
initialization with the same point $x_{0}\in\text{dom}f$. Recall that this
setting allows us to capture reproducibility challenges due to non-
deterministic numerical errors introduced by a computing device (like GPUs)
during floating point computations. The main high-level message here is that
unlike in the stochastic gradient oracle case, in the non-stochastic gradient
oracle setting the iteration complexity $T$ has little to no effect on the
$(\varepsilon,\delta)$-deviation: intuitively this is because unlike the
stochastic setting, it is not possible to reduce the error in the gradients by
taking more samples and averaging.

For $D=O(1)$, $\varepsilon,\delta>0$ such that
$\delta\leq\nicefrac{{\varepsilon}}{{(2D)}}$, and number of iterations $T$,
the $(\varepsilon,\delta)$-deviation for optimizing convex functions whose
optimum has norm at most $D$ with a non-stochastic inexact gradient oracle is
as follows. Unless indicated otherwise, the lower bounds hold for any FOI
algorithm, and the upper bound is achieved by projected gradient descent for
$T=\Omega(\nicefrac{{1}}{{\varepsilon}})$ in the smooth settings and non-
smooth strongly convex setting, and
$T=\Omega(\nicefrac{{1}}{{\varepsilon^{2}}})$ in the non-smooth setting.

  * $\blacksquare$

Smooth functions: $(\varepsilon,\delta)$-deviation is
$\Theta(\nicefrac{{\delta^{2}}}{{\varepsilon^{2}}})$.111Note that Subsection
E.2 provides a similar upper bound result without the assumption that the
optimum lies in a ball of radius $O(1)$.

  * $\blacksquare$

Smooth and strongly convex functions: $(\varepsilon,\delta)$-deviation is
$\Theta(\delta^{2}\wedge\varepsilon)$.

  * $\blacksquare$

Lipschitz (non-smooth) functions: $(\varepsilon,\delta)$-deviation is
$\Theta(\nicefrac{{1}}{{(T\varepsilon^{2})}}+\nicefrac{{\delta^{2}}}{{\varepsilon^{2}}})$.

  * $\blacksquare$

Lipschitz (non-smooth) and strongly convex functions:
$(\varepsilon,\delta)$-deviation is
$\Theta((\nicefrac{{1}}{{T}}+\delta^{2})\wedge\varepsilon)$.

Like Section 3, this theorem is proved in several pieces; see Appendix A. Some
remarks follow:

Acceleration v.s. reproducibility? Note that for smooth functions, accelerated
methods can get $\varepsilon$-suboptimality in only
$T=O(\nicefrac{{1}}{{\sqrt{\varepsilon}}})$ iterations. A natural question is
if we can achieve similar $(\varepsilon,\delta)$-deviation for such
accelerated methods. However, it is known that accelerated methods are
unstable (see e.g., (Devolder et al., 2014; Attia and Koren, 2021)), and
cannot even achieve the desired $\varepsilon$-accuracy under the inexact
oracle model. Hence, we conjecture that iteration complexity of
$\Omega(1/\varepsilon)$ is necessary to achieve the desired reproducibility.

Nonsmoothness: Note that the lower bound shows that if gradient descent style
methods are run for the standard $T=1/\varepsilon^{2}$ iterations, then the
solutions can be forced to have $\Omega(1)$ deviation. Hence, to ensure
reproducibility, any method would need to have a slower convergence rate than
$T=1/\varepsilon^{2}$.

###  5 Reproducibility with Inexact Initialization Oracles

We now study reproducibility for optimizing a function $f\in\mathcal{F}$ with
an inexact initialization oracle (Section 2) and establish tight bounds on the
$(\varepsilon,\delta)$-deviation for FOI algorithms. An inexact initialization
oracle models situations where the initial values of the parameters are set
randomly, or incur some non-deterministic numerical error.

For any $\varepsilon,\delta>0$ and number of iterations $T$, the
$(\varepsilon,\delta)$-deviation for optimizing convex functions with an
inexact initialization oracle is as follows. In all cases, the lower bounds
hold for any FOI algorithm, and the upper bound is achieved by gradient
descent for $T=\Omega(\nicefrac{{1}}{{\varepsilon}})$ in the smooth settings
and non-smooth strongly convex setting, and
$T=\Omega(\nicefrac{{1}}{{\varepsilon^{2}}})$ in the non-smooth setting.

  * $\blacksquare$

Smooth functions: $(\varepsilon,\delta)$-deviation is $\Theta(\delta^{2})$.

  * $\blacksquare$

Smooth and strongly convex functions: $(\varepsilon,\delta)$-deviation is
$\Theta((\exp(-\Omega(T))\delta^{2}\wedge\varepsilon)$.

  * $\blacksquare$

Lipschitz (non-smooth) functions: $(\varepsilon,\delta)$-deviation is
$\Theta(\nicefrac{{1}}{{(T\varepsilon^{2})}}+\delta^{2})$.

  * $\blacksquare$

Lipschitz (non-smooth) and strongly convex functions:
$(\varepsilon,\delta)$-deviation is
$\Theta(\nicefrac{{1}}{{T}}\wedge\varepsilon)$.

Like Section 3, this theorem is proved in several pieces; see Appendix A. A
remark follows:

Nonsmoothness: As in the case of non-smooth optimization with inexact gradient
oracles, here as well we see that even the slightest amount of inexactness in
initialization can lead to non-negligible irreproducibility. Intuitively this
is a consequence of non-smoothness: If the function is not differentiable at
the reference point, even a slight bit of inexactness can lead to drastically
different trajectories right from the beginning.

###  6 Reproducibility in Optimization for Machine Learning

So far, we have studied reproducibility in general convex optimization with
different sources of perturbation. In machine learning, however, optimization
problems come with more structure, and hence a more nuanced analysis of
reproducibility is called for. In this section, we study reproducibility in
two specific optimization settings of interest in machine learning: finite sum
minimization which corresponds to minimizing training loss, and stochastic
convex optimization (SCO), which corresponds to minimizing population (or
test) loss. Instead of a detailed study of different types of convex functions
as done in the previous sections, here we focus on a few specific cases that
yield particularly interesting insights.

####  6.1 Optimizing Training Loss (Finite Sum Minimization)

Optimizing the training loss in ML can be cast as minimizing a function that
can be written as a finite sum of component functions:
$f({\bm{x}}):=\frac{1}{m}\sum_{i=1}^{m}f_{i}({\bm{x}})\,.$ Typical
optimization methods such as stochastic gradient descent can be implemented to
solve such problems by iteratively sampling one of the component functions and
taking its gradient. So randomness in sampling as well as non-deterministic
numerical errors in computing of gradient can lead to irreproducibility.
Randomness in sampling is a property of the algorithm and hence under the
control of the algorithm designer; however non-deterministic numerical errors
in gradient computations are beyond the control of the designer, and hence we
aim to quantify irreproducibility caused by this specific source (numerical
errors). To capture this intuition, we consider the following inexact gradient
oracle.

[$\delta$-bounded inexact component gradient oracle] A $\delta$-bounded
inexact component gradient oracle takes as input $i\in\\{1,2,\ldots,m\\}$ and
a point $x$, and outputs a vector $g_{i}({\bm{x}})$ such that

$\displaystyle\left\lVert g_{i}({\bm{x}})-\nabla
f_{i}({\bm{x}})\right\rVert^{2}\leq\delta^{2}\text{ for some }\nabla
f_{i}({\bm{x}})\in\partial f_{i}({\bm{x}}).$

The number of calls to the component gradient oracle is now our primary
measure of complexity. Then, a natural question that arises is whether we can
we tightly characterize reproduciblity of typical first-order methods when
given access to an inexact component gradient oracle. We first observe that in
the special case that all $f_{i}$ are identical, the inexact component
gradient oracle reduces to the non-stochastic inexact oracle (Section 2).
Hence, for the smooth and non-smooth function settings, we immediately obtain
lower bounds on the $(\varepsilon,\delta)$-deviation via our previous
analysis; see Subsection C.2 and Subsection D.4, respectively.

The non-smooth case in particular is interesting, as the deviation lower bound
is
$\Omega(\nicefrac{{1}}{{(T\varepsilon^{2})}}+\nicefrac{{\delta^{2}}}{{\varepsilon^{2}}})$,
where $T$ is the number of oracle calls. Now, if we use projected (full-batch)
gradient descent (as used for the matching upper bound in Subsection F.2),
then each iteration would require $m$ oracle calls, hence the
$(\varepsilon,\delta)$-deviation for such a method is
$O(\nicefrac{{m}}{{(T\varepsilon^{2})}}+\nicefrac{{\delta^{2}}}{{\varepsilon^{2}}})$,
which is worse than the lower bound mentioned above. Now, the key question is
whether we can reduce this deviation by other means. It turns out that
stochastic gradient descent (SGD), where in each step we sample a component
function to query randomly, when run with an appropriate learning rate and
using averaging, matches the optimal deviation up to constant factors. We
provide a detailed proof of the result in Appendix G. This provides another
justification for using SGD instead of gradient descent in practice.

[] For $G=O(1)$ and $D=O(1)$, let $f_{i}$ be an $G$-Lipschitz convex cost
function for each $i\in[m]$, and assume that the optimum of $f$ lies in a ball
of radius $D$. Let $\varepsilon,\delta>0$ be given parameters such that
$\delta\leq\nicefrac{{\varepsilon}}{{(2D)}}$, and
$T=\Omega(1/\varepsilon^{2})$ be a given number of iterations. Define the SGD
updates as follows: initialize ${\bm{x}}_{0}=0$, and for $t=0,1,\ldots,T-1$,
set ${\bm{x}}_{t+1}={\bm{x}}_{t}-\eta_{t}g_{i_{t}}({\bm{x}}_{t})$ where
$i_{t}\sim[n]$ uniformly at random. Under the inexact component gradient
oracle (Subsection 6.1), the average iterate $\bar{{\bm{x}}}_{T}$ of SGD with
stepsize $\eta=\Theta(\nicefrac{{1}}{{(\varepsilon T)}})$ satisfies
$\operatorname*{\mathbb{E}}f(\bar{{\bm{x}}}_{T})-\inf_{{\bm{x}}\in\text{dom}f}f({\bm{x}})\leq\varepsilon$
and
$\operatorname*{\mathbb{E}}\left\lVert\bar{{\bm{x}}}_{T}-\bar{{\bm{x}}}^{\prime}_{T}\right\rVert^{2}=O(\nicefrac{{1}}{{(T\varepsilon^{2})}}+\nicefrac{{\delta^{2}}}{{\varepsilon^{2}}})$,
where $\bar{{\bm{x}}}^{\prime}_{T}$ is the output of an independent run of
SGD.

####  6.2 Optimizing Population Loss (Stochastic Convex Optimization)

The fundamental problem of machine learning is to find the solution to the
following population (or test) loss minimization problem: minimize
$F({\bm{x}})=\operatorname*{\mathbb{E}}_{\xi\sim\Xi}f({\bm{x}},\xi)$, where
$\Xi$ is an unknown distribution on the examples $\xi$, given access to an
oracle that sample from $\Xi$. When the function $F$ is convex, this is also
known as Stochastic Convex Optimization (SCO).

In this setting, the sampling oracle is a natural source of irreproducibility
of optimization methods. To model this, we consider the _stochastic global
oracle_ , inspired by Foster et al. (2019). [$\delta$-bounded stochastic
global oracle] Given a function $F$, a $\delta$-bounded stochastic global
oracle for $F$ is an algorithm that, on each query, draws an independent
sample $\xi\sim\Xi$ from some distribution $\Xi$, and outputs a function
$f(\cdot,\xi):\text{dom}F\rightarrow\mathbb{R}$ such that for all
$x\in\text{dom}F$, we have
$F({\bm{x}})=\operatorname*{\mathbb{E}}_{\xi\sim\Xi}f({\bm{x}},\xi)$, and,
$\operatorname*{\mathbb{E}}_{\xi\sim\Xi}\left\lVert\nabla
f({\bm{x}},\xi)-\nabla F({\bm{x}})\right\rVert^{2}\leq\delta^{2}\,.$ Since in
each query, the stochastic global oracle returns the complete specification of
$f(\cdot,\xi)$, an algorithm using such an oracle has access to
$({\bm{x}},f({\bm{x}},\xi),\nabla
f({\bm{x}},\xi),\nabla^{2}f({\bm{x}},\xi),\cdots)$ for all ${\bm{x}}$.

Note that the inexact stochastic gradient oracle of Section 2 is a weaker form
of the stochastic global oracle. Thus, the key question is: under the powerful
stochastic global oracle, can we provide better reproducibility than the lower
bounds in Section 3? Surprisingly, the answer to the above question is
negative. That is, despite the stronger oracle setting, there is a function
class for which the lower bound on $(\varepsilon,\delta)$-deviation matches
that of Section 3.

[] Assume that $\varepsilon<1/200$ and $\varepsilon\lesssim\delta^{2}$. Then
there exists a family of population costs
${F}_{\theta}({\bm{x}})=\operatorname*{\mathbb{E}}_{\xi\sim
P_{\theta}}f_{\theta}({\bm{x}},\xi)$ parametrized by $\theta\in[1,2]$, and a
$\delta$-bounded stochastic global oracle for ${F}_{\theta}$, satisfying the
following property. Suppose that $A$ is any algorithm that for each
$\theta\in[1,2]$ uses at most $T$ queries to the stochastic global oracle and
outputs ${\bm{x}}_{T}^{\theta}$ that is $\varepsilon$-accurate, i.e.,
$\operatorname*{\mathbb{E}}{F}_{\theta}({\bm{x}}_{T}^{\theta})-\inf_{{\bm{x}}}{F}_{\theta}({\bm{x}})\leq\varepsilon$.
Then, there exists $\theta_{\sf bad}\in[1,2]$ such that
$\operatorname*{\mathrm{Var}}({\bm{x}}_{T}^{\theta_{\sf
bad}})\geq\Omega(\nicefrac{{\delta^{2}}}{{(T\varepsilon^{2})}})$. See
Subsection B.2 for a detailed proof. On the other hand, Section 3 also shows
that standard SGD with the weaker inexact stochastic gradient oracle achieves
this lower bound.

###  7 Conclusions

We presented a framework for studying reproducibility in convex optimization,
and explored limits of reproducibility for optimizing various function classes
under different sources of errors. For each of these settings, we provide
tight lower and upper bounds on reproducibility of first order iterative
algorithms. Overall, our results provide the following insights: a) Non-smooth
functions can be highly susceptible to even tiny errors which can creep in due
to say numerical errors in GPU. Thus introducing smoothness in deep learning
models might help with reproducibility. b) Generally, gradient descent type
methods with small learning rate are more reproducible. c) In finite-sum
settings, despite more randomness, SGD is more reproducible than full-batch
gradient descent.

The study in this paper is a first step towards addressing the challenging
problem of reproducibility in a rigorous framework. So, many important
directions remain unexplored. For example, study of reproducibility of
adaptive methods like AdaGrad that are not captured by (FOI) is interesting.
We study reproducibility in a strict form where we measure deviation in terms
of the learned parameters of the model. However, in practice, one might care
for reproducibility in model predictions only. So, extending our results to
such ML-driven scenarios should be relevant in practice. Finally, while our
lower bounds on reproducibility obviously also hold for non-convex
optimization, extension of the upper bounds to non-convex settings and
potentially, designing novel rigorous methods in this setting is a fascinating
direction for future work.

### Acknowledgments and Disclosure of Funding

Kwangjun Ahn was supported by graduate assistantship from the NSF Grant
(CAREER: 1846088), the ONR grant (N00014-20-1-2394) and MIT-IBM Watson as well
as a Vannevar Bush fellowship from Office of the Secretary of Defense. Part of
this work was done while Kwangjun Ahn was visiting the Simons Institute for
the Theory of Computing, Berkeley, California.

### References

  * Achille et al. [2017] A. Achille, M. Rovere, and S. Soatto.  Critical learning periods in deep neural networks.  _arXiv preprint arXiv:1711.08856_ , 2017.
  * Allen-Zhu and Li [2020] Z. Allen-Zhu and Y. Li.  Towards understanding ensemble, knowledge distillation and self-distillation in deep learning.  _arXiv preprint arXiv:2012.09816_ , 2020.
  * Anil et al. [2018] R. Anil, G. Pereyra, A. Passos, R. Ormandi, G. E. Dahl, and G. E. Hinton.  Large scale distributed neural network training through online distillation.  _arXiv preprint arXiv:1804.03235_ , 2018.
  * Attia and Koren [2021] A. Attia and T. Koren.  Algorithmic instabilities of accelerated gradient descent.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Aybat et al. [2020] N. S. Aybat, A. Fallah, M. Gurbuzbalaban, and A. Ozdaglar.  Robust accelerated gradient methods for smooth strongly convex functions.  _SIAM Journal on Optimization_ , 30(1):717–751, 2020.
  * Bansal and Gupta [2019] N. Bansal and A. Gupta.  Potential-function proofs for gradient methods.  _Theory of Computing_ , 15(1):1–32, 2019.
  * Bassily et al. [2020] R. Bassily, V. Feldman, C. Guzmán, and K. Talwar.  Stability of stochastic gradient descent on nonsmooth convex losses.  _Advances in Neural Information Processing Systems_ , 33:4381–4391, 2020.
  * Bauschke et al. [2011] H. H. Bauschke, P. L. Combettes, et al.  _Convex analysis and monotone operator theory in Hilbert spaces_ , volume 408.  Springer, 2011.
  * Bhojanapalli et al. [2021] S. Bhojanapalli, K. Wilber, A. Veit, A. S. Rawat, S. Kim, A. Menon, and S. Kumar.  On the reproducibility of neural network predictions.  _arXiv preprint arXiv:2102.03349_ , 2021.
  * Bousquet and Elisseeff [2001] O. Bousquet and A. Elisseeff.  Algorithmic stability and generalization performance.  _Advances in Neural Information Processing Systems_ , pages 196–202, 2001.
  * Bubeck [2014] S. Bubeck.  Convex optimization: algorithms and complexity.  _arXiv preprint arXiv:1405.4980_ , 2014.
  * Chen et al. [2020] Z. Chen, Y. Wang, D. Lin, D. Cheng, L. Hong, E. Chi, and C. Cui.  Beyond point estimate: Inferring ensemble prediction variation from neuron activation strength in recommender systems.  _arXiv preprint arXiv:2008.07032_ , 2020.
  * Cohen et al. [2018] M. Cohen, J. Diakonikolas, and L. Orecchia.  On acceleration with noise-corrupted gradients.  In _International Conference on Machine Learning_ , pages 1019–1028. PMLR, 2018.
  * D’Amour et al. [2020] A. D’Amour, K. Heller, D. Moldovan, B. Adlam, B. Alipanahi, A. Beutel, C. Chen, J. Deaton, J. Eisenstein, M. D. Hoffman, F. Hormozdiari, N. Houlsby, S. Hou, G. Jerfel, A. Karthikesalingam, M. Lucic, Y. Ma, C. McLean, D. Mincu, A. Mitani, A. Montanari, Z. Nado, V. Natarajan, C. Nielson, T. F. Osborne, R. Raman, K. Ramasamy, R. Sayres, J. Schrouff, M. Seneviratne, S. Sequeira, H. Suresh, V. Veitch, M. Vladymyrov, X. Wang, K. Webster, S. Yadlowsky, T. Yun, X. Zhai, and D. Sculley.  Underspecification presents challenges for credibility in modern machine learning.  _arXiv preprint arXiv:2011.03395_ , 2020.
  * d’Aspremont [2008] A. d’Aspremont.  Smooth optimization with approximate gradient.  _SIAM Journal on Optimization_ , 19(3):1171–1183, 2008.
  * Devolder et al. [2014] O. Devolder, F. Glineur, and Y. Nesterov.  First-order methods of smooth convex optimization with inexact oracle.  _Mathematical Programming_ , 146(1):37–75, 2014\.
  * Dietterich [2000] T. G. Dietterich.  Ensemble methods in machine learning.  _Lecture Notes in Computer Science_ , pages 1–15, 2000.
  * Du [2019] S. Du.  _Gradient Descent for Non-convex Problems in Modern Machine Learning_.  PhD thesis, Carnegie Mellon University, 2019.
  * Dusenberry et al. [2020] M. W. Dusenberry, D. Tran, E. Choi, J. Kemp, J. Nixon, G. Jerfel, K. Heller, and A. M. Dai.  Analyzing the role of model uncertainty for electronic health records.  In _Proceedings of the ACM Conference on Health, Inference, and Learning_ , pages 204–213, 2020.
  * Dwork et al. [2006] C. Dwork, F. McSherry, K. Nissim, and A. D. Smith.  Calibrating noise to sensitivity in private data analysis.  In _TCC_ , volume 3876 of _Lecture Notes in Computer Science_ , pages 265–284. Springer, 2006.
  * Eldar [2004] Y. C. Eldar.  Minimum variance in biased estimation: Bounds and asymptotically optimal estimators.  _IEEE Transactions on Signal Processing_ , 52(7):1915–1930, 2004.
  * Fort et al. [2020] S. Fort, H. Hu, and B. Lakshminarayanan.  Deep ensembles: A loss landscape perspective.  _arXiv preprint arXiv:1912.02757_ , 2020.
  * Foster et al. [2019] D. J. Foster, A. Sekhari, O. Shamir, N. Srebro, K. Sridharan, and B. Woodworth.  The complexity of making the gradient small in stochastic convex optimization.  In _Conference on Learning Theory_ , pages 1319–1345. PMLR, 2019\.
  * Frankle et al. [2020] J. Frankle, G. K. Dziugaite, D. Roy, and M. Carbin.  Linear mode connectivity and the lottery ticket hypothesis.  In _International Conference on Machine Learning_ , pages 3259–3269. PMLR, 2020.
  * Impagliazzo et al. [2022] R. Impagliazzo, R. Lei, T. Pitassi, and J. Sorrell.  Reproducibility in learning.  _To appear at STOC 2022 arXiv:2201.08430_ , 2022.
  * Kutin and Niyogi [2002] S. Kutin and P. Niyogi.  Almost-everywhere algorithmic stability and generalization error.  In _UAI_ , pages 275–282, 2002.
  * Lakshminarayanan et al. [2017] B. Lakshminarayanan, A. Pritzel, and C. Blundell.  Simple and scalable predictive uncertainty estimation using deep ensembles.  In _Advances in neural information processing systems_ , pages 6402–6413, 2017.
  * Lehmann and Romano [2005] E. Lehmann and J. Romano.  _Testing statistical hypotheses_.  Springer, 2005.
  * McSherry and Talwar [2007] F. McSherry and K. Talwar.  Mechanism design via differential privacy.  In _FOCS_ , pages 94–103. IEEE Computer Society, 2007.
  * Mhaskar [1997] H. Mhaskar.  On smooth activation functions.  In _Mathematics of Neural Networks_ , pages 275–279. Springer, 1997\.
  * Nair and Hinton [2010] V. Nair and G. E. Hinton.  Rectified linear units improve restricted Boltzmann machines.  In _ICML_ , 2010.
  * Nemirovski and Yudin [1983] A. S. Nemirovski and D. B. Yudin.  _Problem complexity and method efficiency in optimization_.  Wiley-Interscience, 1983.
  * Nesterov [2018] Y. Nesterov.  _Lectures on convex optimization_ , volume 137.  Springer, 2018.
  * Pineau et al. [2021] J. Pineau, P. Vincent-Lamarre, K. Sinha, V. Larivière, A. Beygelzimer, F. d’Alché Buc, E. Fox, and H. Larochelle.  Improving reproducibility in machine learning research: a report from the NeurIPS 2019 reproducibility program.  _Journal of Machine Learning Research_ , 22, 2021.
  * Shallue et al. [2018] C. J. Shallue, J. Lee, J. Antognini, J. Sohl-Dickstein, R. Frostig, and G. E. Dahl.  Measuring the effects of data parallelism on neural network training.  _arXiv preprint arXiv:1811.03600_ , 2018.
  * Shamir [2018] G. I. Shamir.  Systems and methods for improved generalization, reproducibility, and stabilization of neural networks via error control code constraints, 2018.
  * Shamir and Coviello [2020] G. I. Shamir and L. Coviello.  Anti-distillation: Improving reproducibility of deep networks.  _arXiv preprint arXiv:2010.09923_ , 2020.
  * Shamir et al. [2020] G. I. Shamir, D. Lin, and L. Coviello.  Smooth activations and reproducibility in deep networks.  _arXiv preprint arXiv:2010.09931_ , 2020.
  * Snapp and Shamir [2021] R. R. Snapp and G. I. Shamir.  Synthesizing irreproducibility in deep networks.  _arXiv preprint arXiv:2102.10696_ , 2021.
  * Sontag and Wang [1995] E. D. Sontag and Y. Wang.  On characterizations of the input-to-state stability property.  _Systems & Control Letters_, 24(5):351–359, 1995.
  * Summers and Dinneen [2021] C. Summers and M. J. Dinneen.  On nondeterminism and instability in neural network optimization, 2021\.
  * Tu et al. [2022] S. Tu, A. Robey, T. Zhang, and N. Matni.  On the sample complexity of stability constrained imitation learning.  _Proceedings of the 4rd Conference on Learning for Dynamics and Control, arXiv preprint arXiv:2102.09161_ , 2022.
  * Yu et al. [2021] H. Yu, Z. Chen, D. Lin, G. Shamir, and J. Han.  Dropout prediction variation estimation using neuron activation strength.  _arXiv preprint arXiv:2110.06435_ , 2021.
  * Zhuang et al. [2021] D. Zhuang, X. Zhang, S. L. Song, and S. Hooker.  Randomness in neural network training: Characterizing the impact of tooling.  _arXiv preprint arXiv:2106.11872_ , 2021.

## Appendix

###  Appendix A Summary of results in Appendix

In order to help readers navigates the results in the appendix, we summarize
the results in the appendix in the following table. In the main paper, we have
three main theorems: Theorems 3, 4 and 5, each corresponding to one column in
the following table. Further each of these theorems have four components
corresponding to four settings: smooth, smooth & strongly convex, nonsmooth
and finally nonsmooth & strongly convex. Each cell in the below table lists
the corresponding theorems which prove the lower and upper bounds
corresponding to a given setting. Theorems 3, 4 and 5 follow immediately from
the constituent theorems. In each cell of results, “Theorems X$\parallel$ Y”
indicates that the lower bound appears in Theorem X, and the upper bound in
Theorem Y. All lower bounds are for first-order iterative algorithms (_à la_
Nesterov [2018]) that we formally defined in (FOI).

| Stochastic Inexact | Non-stochastic Inexact | Inexact Initialization
---|---|---|---
Gradient Oracle | Gradient Oracle | Oracle
(Section 3) | (Section 4) | (Section 5)
Smooth | $\Theta(\nicefrac{{\delta^{2}}}{{T\varepsilon^{2}}})$ | $\Theta(\nicefrac{{\delta^{2}}}{{\varepsilon^{2}}})$ | $\Theta(\delta^{2})$
Theorems B.1†&C.1$\parallel$ E.1 | Theorems C.2$\parallel$ E.2 | Theorems C.3$\parallel$ E.3
Smooth | $\Theta(\nicefrac{{\delta^{2}}}{{T}}\wedge\varepsilon)$ | $\Theta(\delta^{2}\wedge\varepsilon)$ | $\Theta(e^{-\Omega(T)}\delta^{2}\wedge\varepsilon)$
Strongly-Convex | Theorems C.4$\parallel$ E.4 | Theorems C.5$\parallel$ E.5 | Theorems C.6$\parallel$ E.6
Nonsmooth | $\Theta(\nicefrac{{1}}{{T\varepsilon^{2}}})$ | $\Theta(\nicefrac{{1}}{{T\varepsilon^{2}}}+\nicefrac{{\delta^{2}}}{{\varepsilon^{2}}})$ | $\Theta(\nicefrac{{1}}{{T\varepsilon^{2}}}+\delta^{2})$
Theorems D.3$\parallel$ F.1 | Theorems D.4$\parallel$ F.2 | Theorems D.5$\parallel$ F.3
Nonsmooth | $\Theta(\nicefrac{{1}}{{T}}\wedge{\varepsilon})$ | $\Theta((\nicefrac{{1}}{{T}}+{\delta^{2}})\wedge{\varepsilon})$ | $\Theta(\nicefrac{{1}}{{T}}\wedge{\varepsilon})$
Strongly-Convex | Theorems D.6$\parallel$ F.4 | Theorems D.7$\parallel$ F.5 | Theorems D.8$\parallel$ F.6

${\dagger}$ For smooth costs and stochastic inexact gradient oracle, we also
have an additional information-theoretic lower bound of
$\Omega(\nicefrac{{\delta^{2}}}{{T\varepsilon^{2}}})$ when
$\varepsilon\lesssim\delta^{2}$ in Subsection B.1.

####  A.1 General guidance for navigating Appendix

$\blacksquare$ The information-theoretic lower bounds are presented in
Appendix B and can be read independently: $\bullet$ The lower bound proof for
the smooth costs (Subsection B.1) and is presented in Subsection B.1.
$\bullet$ The lower bound proof for the stochastic global oracle case
(Subsection 6.2) is similar to that of Subsection B.1 and is presented in
Subsection B.2. $\blacksquare$ The FOI lower bounds are quite technical and
often rely on delicate constructions. Hence, we present the proofs as follows:
$\bullet$ In Subsection C.1, to help readers understand the general proof
strategy, we first present the lower bound proof for smooth costs with the
stochastic inexact oracle (Subsection C.1). $\bullet$ In Appendix C, we
present the proofs of other FOI lower bounds for smooth costs. $\bullet$ In
Appendix D, we present the proofs of FOI lower bounds for nonsmooth costs. The
lower bound constructions are more complicated than the case of smooth costs.
Hence, to help reader understand the crux of arguments, we first present the
proof of a (weaker) lower bound against gradient descent (as opposed to the
entire class of FOI) in Subsection D.1. $\blacksquare$ The upper bounds are
presented in Appendix E (smooth costs), Appendix F (non-smooth costs), and
Appendix G (finite-sum setting).

###  Appendix B Information-theoretic lower bounds

####  B.1  Information-theoretic lower bound for stochastic inexact gradient
model

We state and prove the information-theoretic lower bound for smooth costs.

[] (Information-theoretic Lower Bound) Assume that $\varepsilon<1/200$ and
$\varepsilon\lesssim\delta^{2}$. Then there exists a family of smooth cost
functions $\\{{F}_{\theta}:\mathbb{R}\to\mathbb{R}\\}$ parameterized by
$\theta\in[-1,1]$ with the following property. Suppose $A$ is any algorithm
that for each $\theta\in[-1,1]$ uses at most $T$ queries to stochastic inexact
gradient oracle and outputs ${\bm{x}}_{T}^{\theta}$ that is
$\varepsilon$-accurate, i.e.,
$\operatorname*{\mathbb{E}}{F}_{\theta}({\bm{x}}_{T}^{\theta})-\inf_{\bm{x}}{F}_{\theta}({\bm{x}})\leq\varepsilon$.
Then, there exists $\theta_{\sf bad}\in[-1,1]$ such that
$(\varepsilon,\delta)$-deviation (Section 2 (a)) is lower bounded by
$\Omega(\frac{\delta^{2}}{T\varepsilon^{2}})$.

###### Proof.

Let $\varepsilon>0$ be a fixed small constant. We consider the following
family of cost functions $\\{{F}_{\theta}:\mathbb{R}\to\mathbb{R}\\}$
parametrized by $\theta\in[-1,1]$:

$\displaystyle{F}_{\theta}(x)=\begin{cases}100\varepsilon\cdot(x-\theta)^{2},&\text{for}~{}x\in[-1,1],\\\
100\varepsilon\cdot(1-\theta)^{2}+200\varepsilon(1-\theta)(x-1)&\text{for}~{}x>1\\\
100\varepsilon\cdot(-1-\theta)^{2}+200\varepsilon(-1-\theta)(x+1)&\text{for}~{}x<-1.\end{cases}$
(B.1)

Note that ${F}_{\theta}$ is ${200\epsilon}$-smooth for each $\theta\in[-1,1]$.

We consider the following stochastic first order oracle for each
${F}_{\theta}$. For a queried point $x$, the oracle outputs $g_{\theta}(x)$
defined as

$\displaystyle g_{\theta}(x)$
$\displaystyle=\begin{cases}\frac{1}{200\varepsilon}\nabla{F}_{\theta}(x)+z\quad\text{for
}z\sim
N(0,\frac{\delta^{2}}{2\cdot{200\epsilon}}),&\text{w.p.}~{}{200\epsilon}\,,\\\
0,&\text{w.p.}~{}1-{200\epsilon}\,.\end{cases}$ (B.2)

Since we define ${F}_{\theta}$ as a linear extension outside $[-1,1]$, below
we may assume that all gradient queries are made within $[-1,1]$.

Assuming that all queries are made in $[-1,1]$, we can rewrite (B.2) simply as

$\displaystyle g_{\theta}(x)=(x-\theta+z)\cdot{\sf Sample}\,,$ (B.3)

where ${\sf Sample}\sim\operatorname{Bernoulli}({200\epsilon})$. Let us verify
that this is a valid stochastic first order oracle. First, it is clear that

$\displaystyle\operatorname*{\mathbb{E}}[g_{\theta}(x)]={200\epsilon}\cdot(x-\theta+0)+(1-{200\epsilon})\cdot
0={200\epsilon}\cdot(x-\theta)=\nabla{F}_{\theta}(x).$

Next, for the variance, note that

$\displaystyle\operatorname*{\mathbb{E}}\left[g_{\theta}(x)-\nabla{F}_{\theta}(x)\right]^{2}$
$\displaystyle=\operatorname*{\mathbb{E}}\left[(x-\theta+z)\cdot{\sf
Sample}-\nabla{F}_{\theta}(x)\right]^{2}$
$\displaystyle={200\epsilon}\cdot\operatorname*{\mathbb{E}}\left[x-\theta+z-{200\epsilon}\cdot(x-\theta)\right]^{2}+(1-{200\epsilon})\cdot\left[-{200\epsilon}\cdot(x-\theta)\right]^{2}$
$\displaystyle={200\epsilon}(1-{200\epsilon})^{2}\cdot(x-\theta)^{2}+({200\epsilon})^{2}(1-{200\epsilon})\cdot(x-\theta)^{2}+{200\epsilon}\cdot\operatorname*{\mathbb{E}}[z^{2}]$
$\displaystyle={200\epsilon}(1-{200\epsilon})\cdot(x-\theta)^{2}+\delta^{2}/2.$

Hence, the variance is always upper bounded by
$\delta^{2}/2+{200\epsilon}(1-{200\epsilon})2^{2}\leq\delta^{2}$, as long as
$\varepsilon\lesssim\delta^{2}$.

We now prove the theorem. From the fact that the output $x^{\theta}_{T}$ is
$\varepsilon$-accurate, we have

$\displaystyle\forall\theta\in[-1,1],\quad
100\varepsilon\cdot\operatorname*{\mathbb{E}}\left[x^{\theta}_{T}-\theta\right]^{2}\leq\varepsilon.$

Using Jensen’s inequality, we know
$\left|\operatorname*{\mathbb{E}}[x^{\theta}_{T}]-\theta\right|^{2}=\left|\operatorname*{\mathbb{E}}[x^{\theta}_{T}-\theta]\right|^{2}\leq\operatorname*{\mathbb{E}}[x^{\theta}_{T}-\theta]^{2}$
which implies the following condition:

$\displaystyle\forall\theta\in[-1,1],\quad\left|\operatorname*{\mathbb{E}}[x^{\theta}_{T}]-\theta\right|\leq
0.1.$ (B.4)

This condition says if we regard $A$ as an estimator of $\theta$ for each
$\theta$, the bias is less than equal to $0.1$. For the following argument, we
hence change our perspective and regard $A$ as an estimator of $\theta$ based
on $T$ inexact gradient queries rather than an optimization algorithm.

Let us fix $\theta\in[-1,1]$. We first that we may assume that all gradient
queries are made at the point $0$. Indeed, from the expression for the
stochastic oracle (B.3), we know

$\displaystyle g_{\theta}(x)\overset{d}{=}g_{\theta}(0)+x\cdot{\sf Sample}.$

In particular, this implies that one can reconstruct a gradient query at any
point $x$ based on a gradient query at the point $0$. Hence, without loss of
generality, we may assume that all gradient queries are made at $x=0$.

Hence, we can regard $A$ as an estimator of $\theta$ based on $T$ independent
measurements $y_{1},\dots,y_{T}$ of form

$\displaystyle y_{i}=(-\theta+z_{i})\cdot{\sf Sample}_{i},\quad
i=1,\dots,T\,,$ (B.5)

where $z_{i}\sim N(0,\frac{\delta^{2}}{{200\epsilon}})$ and ${\sf
Sample}_{i}\sim\operatorname{Bernoulli}({200\epsilon})$. Now with this new
perspective in mind, we can lower bound the variance of the estimator
$\operatorname*{\mathrm{Var}}({\bm{x}}_{T}^{\theta})$ using the Cramer-Rao
lower bound. To that end, we first calculate the fisher information of the
measurement distribution.

Recall that each measurement is of form

$\displaystyle y=(-\theta+z)\cdot{\sf Sample}\,,$

where $z\sim N(0,\frac{\delta^{2}}{2\cdot{200\epsilon}})$ and ${\sf
Sample}\sim\operatorname{Bernoulli}({200\epsilon})$. Then, the log likelihood
is

$\displaystyle\ell(\theta;y)$
$\displaystyle=\ln\left[\frac{1}{\sqrt{2\pi\cdot\frac{\delta^{2}}{2\cdot{200\epsilon}}}}\exp\left(-\frac{2\cdot{200\epsilon}}{\delta^{2}}\cdot\frac{(y+\theta)^{2}}{2}\right)\right]\mathbbm{1}\left\\{{\sf
Sample}=1\right\\}+\delta_{[y=0]}\cdot\mathbbm{1}\left\\{{\sf
Sample}=0\right\\}\,.$

Taking derivatives of the log likelihood, we get

$\displaystyle\nabla_{\theta}\ell(\theta;y)$
$\displaystyle=-\frac{2\cdot{200\epsilon}}{\delta^{2}}\cdot(\theta+y)\cdot\mathbbm{1}\left\\{{\sf
Sample}=1\right\\}\,,$ $\displaystyle\nabla^{2}_{\theta}\ell(\theta;y)$
$\displaystyle=-\frac{2\cdot{200\epsilon}}{\delta^{2}}\cdot\mathbbm{1}\left\\{{\sf
Sample}=1\right\\}\,.$

Hence, the Fisher information is equal to

$\displaystyle
I(\theta)=\operatorname*{\mathrm{Cov}}_{\theta}\nabla_{\theta}\ell(\theta;y)$
$\displaystyle=-\operatorname*{\mathbb{E}}_{\theta}\nabla^{2}_{\theta}\ell(\theta;y)=\frac{2\cdot{200\epsilon}}{\delta^{2}}\cdot\Pr[{\sf
Sample}=1]=\frac{2\cdot({200\epsilon})^{2}}{\delta^{2}}\,.$

Hence, the fisher information $I_{T}(\theta)$ for the $T$ independent
measurements is equal to $T\frac{2\cdot({200\epsilon})^{2}}{\delta^{2}}\cdot$.

Let us recall the Cramér-Rao lower bound for biased estimators. (see, e.g.,
[Eldar, 2004, (3)]). Let
$b(\theta):=\operatorname*{\mathbb{E}}[\hat{\theta}]-\theta$ be the bias of an
estimator $\hat{\theta}$. Then, the following bounds hold:

$\displaystyle\operatorname*{\mathrm{Var}}(\hat{\theta})\geq\frac{[1+b^{\prime}(\theta)]^{2}}{I(\theta)}.$

From (B.4), it must be that $b^{\prime}(\theta_{\sf bad})\geq-\frac{1}{2}$ for
some $\theta_{\sf bad}\in[-1,1]$. To see this, suppose to the contrary that
$b^{\prime}(\theta)<-\frac{1}{2}$ for all $\theta\in[-1,1]$. Then, it must be
that

$\displaystyle b(1)\leq-\frac{1}{2}\cdot
2+b(-1)\overset{\eqref{condition:1}}{\leq}-1+0.1=-0.9,$

which is a contradiction since (B.4) ensures that $b(1)\geq-0.1$. Thus,
Proposition B.1 gives

$\displaystyle\operatorname*{\mathrm{Var}}({\bm{x}}_{T}^{\theta_{\sf
bad}})\geq\frac{(1+b^{\prime}(\theta_{\sf bad}))^{2}}{I_{T}(\theta_{\sf
bad})}\geq\Omega\left(\frac{\delta^{2}}{T\varepsilon^{2}}\right).$

This concludes the proof of the lower bound. ∎

####  B.2 Proof of lower bound (stochastic global oracle)

Recall Subsection 6.2 from the main text. See 6.2

###### Proof.

The construction and argument are analogous to the proof of Subsection B.1
(Subsection B.1). Fix $\varepsilon>0$ and consider the following family of
cost functions $\\{{F}_{\theta}:\mathbb{R}\to\mathbb{R}\\}$ parametrized by
unknown ${\theta\in[1,2]}$:

$\displaystyle{F}_{\theta}(x)=\begin{cases}200\varepsilon\cdot\left\\{\frac{1}{2}x^{2}-\theta
x\right\\},&\text{for}~{}x\in[1,2]\\\ \text{linear
extension}&\text{for}~{}x\not\in[1,2].\end{cases}$ (B.6)

Note that ${F}_{\theta}$ is ${200\epsilon}$-smooth for all $\theta\in[-1,1]$
and the minimum is achieved at $x=\theta$. Below, let us fix a ground truth
parameter $\theta$ and let $F(x)={F}_{\theta}(x)$.

Now define $f(x,\xi)$ as follows: with probability $200\varepsilon$,

$\displaystyle
f(x,\xi)=\begin{cases}\frac{1}{200\varepsilon}{F}(x)+zx\quad\text{for }z\sim
N(0,\frac{\delta^{2}}{2\cdot{200\epsilon}}),&\text{if}~{}x\in[1,2]\,,\\\
\frac{1}{200\varepsilon}{F}(x)+2z.&\text{if}~{}x>2\,,\\\
\frac{1}{200\varepsilon}{F}(x)+z&\text{if}~{}x<1\,,\end{cases}$ (B.7)

and $f(x,\xi)=0$ with probability $1-200\varepsilon$. Then clearly we have
$\operatorname*{\mathbb{E}}_{\xi}f(x,\xi)=F(x)$.

We first check that this construction satisfies Subsection 6.2. It is
sufficient to check the condition for $x\in[1,2]$ since outside the interval
the cost is defined as the linear extension. For $x\in[1,2]$, we have

$\displaystyle\nabla f(x,\xi)$
$\displaystyle=\begin{cases}x-\theta+z\quad\text{for }z\sim
N(0,\frac{\delta^{2}}{2\cdot{200\epsilon}}),&\text{w.p.}~{}{200\epsilon}\,,\\\
0,&\text{w.p.}~{}1-{200\epsilon}\,.\end{cases}$ (B.8)

This is precisely the expression (B.2) in the proof of Subsection B.1
(Subsection B.1), and hence, this clearly satisfies Subsection 6.2.

Now the key fact of the proof is that one can reconstruct the complete
specification of the function $f(\cdot,\xi)$ based on a gradient query,
provided that it is nonzero. This is because if nonzero, the gradient query at
$x$ is equal to $(x-\theta+z)$. This reveals $(\theta-z)$, from which one can
reconstruct the complete characterization
$f(x,\xi)=\frac{1}{2}x^{2}-(\theta-z)x$.

Hence, the information revealed by a single query to the stochastic global
oracle is as good as that revealed by a single query to the stochastic global
oracle. Thus, the setting is reduced to that of Subsection B.1 (Subsection
B.1), and using the same argument, the proof follows. ∎

###  Appendix C Proof of lower bounds (smooth costs)

We first introduce a helper function we will use throughout the proofs of FOI
lower bounds in the remaining sections.

##### Helper function for smooth costs lower bounds.

We will frequently use the following function for the FOI lower bounds for
smooth cost. Let $\mathcal{F}:\mathbb{R}\to\mathbb{R}$ be an one-dimensional
function defined as

$\displaystyle\mathcal{F}(x):=\begin{cases}x^{2}&\text{if }x\in[0,1],\\\
2x-1&\text{if }x\geq 1,\\\ 0&\text{if }x\leq 0.\end{cases}$ (C.1)

For reader’s convenience, we illustrate the helper function $\mathcal{F}$ in
Figure 1 below.

$-1$$0$$1$$2$$3$$0$$1$$2$$3$$4$$5$$0.5$$0.25$ Figure 1: Illustration of the
helper function $\mathcal{F}$ for the smooth costs lower bounds.

####  C.1 Stochastic inexact gradient model

The proofs of FOI lower bounds are quite technical and rely on delicate
constructions, and to illustrate our general proof strategy, we first present
a proof that is relatively simpler, yet captures the essence of the later
complicated constructions. More specifically, in this section, we will prove
the following FOI lower bounds for smooth costs against stochastic inexact
oracle.

[] (Lower Bound) Let $\varepsilon>0$, and $T$ be a given number of iterations.
There exists an $O(1)$-smooth convex function
$f:\mathbb{R}^{O(T)}\to\mathbb{R}$ and a stochastic inexact gradient oracle
such that any FOI algorithm $\mathcal{A}$ that starts at
${\bm{x}}_{0}={\bm{0}}$ has its $(\varepsilon,\delta)$-deviation lower bounded
by $\Omega(\frac{\delta^{2}}{T\varepsilon^{2}})$.

###### Proof.

Now consider the cost $f:{\bm{x}}=({\bm{x}}^{\sf
dum},y)\in\mathbb{R}^{T}\times\mathbb{R}\to\mathbb{R}$ defined as

$\displaystyle f({\bm{x}})=4\varepsilon\cdot\mathcal{F}(y+1)\,,$ (C.2)

where $\mathcal{F}$ is defined in (C.1). Here note that ${\bm{x}}^{\sf
dum}\in\mathbb{R}^{T}$ is dummy coordinates which do not appear in the cost
$f$. Next, we define the stochastic inexact oracle for $t=0,1,\dots,T-1$ as

$\displaystyle g({\bm{x}}_{t})=\nabla f({\bm{x}}_{t})+\delta
r_{t}{\bm{e}}_{1+t}\quad\text{where}~{}r_{t}\sim\operatorname{Unif}\\{\pm
1\\}\,,$ (C.3)

where ${\bm{e}}_{j}$ is the $j$-th coordinate vector. Then, clearly this
stochastic gradient fulfills the definition of stochastic inexact gradient
oracle. Here the stochastic gradient noises are designed such that during the
$t$-th iteration, the noise is added to the coordinate ${\bm{x}}^{\sf
dum}[1+t]$. In other words, the noises will be added to each coordinate of $u$
incrementally.

From here one, let us write iterates ${\bm{x}}_{t}=({\bm{x}}^{\sf
dum}_{t},y_{t})$. Note first that for $\varepsilon$-accuracy, it must be that
$|y_{T}-y_{0}|\geq 1/2$; otherwise
$f({\bm{x}}_{T})>4\varepsilon\cdot\mathcal{F}(0.5)=\varepsilon$; see Figure 1.
Based on the definition of FOI (see (FOI)), let us write

$\displaystyle{\bm{x}}_{T}={\bm{x}}_{0}-\sum_{t=0}^{T-1}\lambda_{t}^{(T)}g({\bm{x}}_{t})\,.$

Then from the construction (C.2), we know that for any ${\bm{x}}$, we know
$\frac{\partial f}{\partial y}({\bm{x}})\in[0,8\varepsilon]$. On the other
hand, as we discussed, we need $|y_{T}-y_{0}|\geq 1/2$. Hence, in order for
iterates to move far enough from the starting point, the coefficients have to
add up to a sufficiently large number:

$\displaystyle\sum_{t=0}^{T-1}|\lambda_{t}^{(T)}|\geq\frac{1}{16\varepsilon}\,,$
(C.4)

since otherwise, $|y_{T}-y_{0}|<8\varepsilon\cdot\frac{1}{16\varepsilon}=1/2$.

Now we will make use of the condition (C.4) to show that there is a large
deviation in the coordinates ${\bm{x}}^{\sf dum}$. More specifically, let us
lower bound $\operatorname*{\mathbb{E}}\left\lVert{\bm{x}}^{\sf
dum}_{T}-\operatorname*{\mathbb{E}}[{\bm{x}}^{\sf dum}_{T}]\right\rVert^{2}$.
From the construction of inexact oracle (C.3), it follows that

$\displaystyle\operatorname*{\mathbb{E}}\left\lVert{\bm{x}}^{\sf
dum}_{T}-\operatorname*{\mathbb{E}}[{\bm{x}}^{\sf dum}_{T}]\right\rVert^{2}$
$\displaystyle=\operatorname*{\mathbb{E}}\left\lVert\sum_{t=0}^{T-1}\lambda_{t}^{(T)}\delta
r_{t}\cdot{\bm{e}}_{1+t}\right\rVert^{2}=\sum_{t=0}^{T-1}(\lambda_{t}^{(T)})^{2}\delta^{2}\operatorname*{\mathbb{E}}[r_{t}^{2}]$
$\displaystyle=\sum_{t=0}^{T-1}(\lambda_{t}^{(T)})^{2}\delta^{2}\overset{(a)}{\geq}\delta^{2}\cdot\frac{1}{T}\cdot\left(\sum_{t=0}^{T-1}|\lambda_{t}^{(T)}|\right)^{2}\gtrsim\frac{\delta^{2}}{T\varepsilon^{2}},$

where ($a$) is due to Cauchy-Schwarz inequality. This concludes the proof. ∎

####  C.2 Non-stochastic inexact gradient model

[] (Lower Bound) Let $\varepsilon>0$ be a small constant, and $T$ be a given
number of iterations. There exists a $O(1)$-smooth convex function
$f:\mathbb{R}^{2}\to\mathbb{R}$ with a non-stochastic inexact gradient model
such that for any FOI algorithm $\mathcal{A}$ that starts at
${\bm{x}}_{0}={\bm{0}}$ has the $(\varepsilon,\delta)$-deviation lower bounded
by $\Omega(\frac{\delta^{2}}{\varepsilon^{2}})$.

###### Proof.

With $\mathcal{F}$ defined as (C.1), this time we consider a simpler
construction: the cost
$f:{\bm{x}}=(x,y)\in\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ is defined as

$\displaystyle f(x,y)=4\varepsilon\cdot\mathcal{F}(y+1).$ (C.5)

Let us write the iterate as ${\bm{x}}_{t}=(x_{t},y_{t})$. Note that for
$\varepsilon$-accuracy, it must be that $y_{T}\geq 1/2$; otherwise
$f(y_{T})>\varepsilon$. This means that in order to achieve
$\varepsilon$-suboptimality, the $y$-component of the iterate has to move at
least constant distance away from the starting point.

Now consider the following non-stochastic inexact oracle

$\displaystyle g({\bm{x}}_{t})=\nabla
f({\bm{x}}_{t})+\delta\cdot\frac{\partial}{\partial
y}\mathcal{F}(y_{t})\cdot{\bm{e}}_{1},$ (C.6)

where ${\bm{e}}_{1}$ is the first coordinate vector. Note that this is a valid
oracle because $0\leq\frac{\partial}{\partial
y}\mathcal{F}(y)\in[0,8\varepsilon]\in[0,1]$ for all $y$. Then from the
construction of the inexact gradient oracle (C.6), it follows that

$\displaystyle\frac{\delta}{4\varepsilon}y_{T}=x_{T}$ (C.7)

Letting ${\bm{x}}_{T}^{\sf exact}=(x_{t}^{\sf exact},y_{t}^{\sf exact})$ be
the iterate with exact gradients (without the gradient noises), since we know
$x_{T}^{\sf exact}=0$, (C.7) implies

$\displaystyle\left\lVert{\bm{x}}_{T}-{\bm{x}}_{T}^{\sf
exact}\right\rVert^{2}$ $\displaystyle\geq\left|x_{T}-x_{T}^{\sf
exact}\right|^{2}=x_{T}^{2}\gtrsim\frac{\delta^{2}}{\varepsilon^{2}}.$

This is the desired lower bound. ∎

####  C.3 Inexact initialization model

[] (Lower Bound) Let $\varepsilon>0$ be a small constant, and $T$ be a given
number of iterations. There exists a $O(1)$-smooth convex function
$f:\mathbb{R}^{O(T)}\to\mathbb{R}$ such that for any FOI algorithm
$\mathcal{A}$ the $(\varepsilon,\delta)$-deviation lower bounded by
$\Omega(\delta^{2})$ w.r.t. the reference point ${\bm{x}}_{0}^{\sf
ref}={\bm{0}}$.

###### Proof.

Consider $f:\mathbb{R}^{2}\to\mathbb{R}$ defined as $f(x,y)=(y-1)^{2}$. Choose
$x_{0}^{\sf ref}=0$ and the inexact initialization to be $x_{0}=(\delta,0)$.
Then, any first order algorithm only updates the second coordinate, which
implies that after $T$ iterations, we still have $\left\lVert x_{T}-x_{T}^{\sf
ref}\right\rVert\geq\delta^{2}$. ∎

####  C.4 Stochastic inexact gradient model (strongly convex costs)

[] (Lower Bound) Let $\varepsilon>0$ and $T$ be a given number of iterations.
There exists a $O(1)$-smooth and $\mu$-strongly convex function
$f:\mathbb{R}^{O(T)}\to\mathbb{R}$ with a stochastic inexact gradient model
such that any FOI algorithm $\mathcal{A}$ that starts at
${\bm{x}}_{0}={\bm{0}}$ has its $(\varepsilon,\delta)$-deviation lower bounded
by $\Omega(\frac{\delta^{2}}{T\mu^{2}}\wedge\frac{\varepsilon}{\mu})$.

###### Proof.

For ${\bm{x}}=({\bm{x}}^{\sf dum},y)$ where ${\bm{x}}^{\sf
dum}\in\mathbb{R}^{T}$ and $y\in\mathbb{R}$, consider the cost

$\displaystyle
f(x,y):=y+\frac{\mu}{2}y^{2}+\frac{\mu}{2}\left\lVert{\bm{x}}\right\rVert^{2}\,.$

We consider the initialization ${\bm{x}}_{0}=(0,0,0,\dots,0)$.

Next, we define the stochastic inexact oracle for $t=0,1,\dots,T-1$ as

$\displaystyle\bm{g}({\bm{x}}_{t})=\nabla f({\bm{x}}_{t})+\delta
r_{t}\cdot{\bm{e}}_{1+t}\quad\text{for}~{}r_{t}\sim\operatorname{Unif}\\{\pm
1\\}\,,$ (C.8)

where ${\bm{e}}_{j}$ is the $j$-th coordinate vector. Also, throughout the
proof we use the notation:

$\displaystyle{\bm{g}}^{\sf dum}_{t}:=\bm{g}({\bm{x}}^{\sf
dum}_{t},y_{t})[1,2,\dots,T]\quad\text{and}\quad
g^{y}_{t}:=\bm{g}({\bm{x}}^{\sf dum}_{t},y_{t})[T+1]\,.$

##### Warm-up: the case of simplified gradient noises.

For a moment, we assume that the inexact gradient oracle is non-stochastic
with

$\displaystyle\bm{g}({\bm{x}}_{t})=\nabla
f({\bm{x}}_{t})+\delta\cdot{\bm{e}}_{1+t}\,.$

We consider this case first to build the key intuition of the proof. The first
prove the following result that is crucial for the proof.

For each $t=0,1,\dots,T-1$, the output of a FOI algorithm satisfies

$\displaystyle\delta\cdot
y_{t}=\sum_{i=1}^{T}x_{t}[i]\quad\text{and}\quad\delta\cdot
g^{y}_{t}=\sum_{i=1}^{T}{\bm{g}}^{\sf dum}_{t}[i]\,,\quad\text{for each
}t=0,1,2,\dots,T.$ (C.9)

###### Proof.

We prove by induction on $t$. The statement trivially holds for $t=0$. Assume
that the conclusion holds for some $t$. We will first show that

$\displaystyle\delta\cdot y_{t+1}$ $\displaystyle=\sum_{i=1}^{T}x_{t+1}[i]\,.$
(C.10)

By the definition of FOI (see (FOI)), we have

$\displaystyle\delta\cdot y_{t+1}$
$\displaystyle=-\delta\cdot\sum_{j=0}^{t}\lambda^{(t+1)}_{j}g^{y}_{j}$
$\displaystyle=-\sum_{j=0}^{t}\lambda^{(t+1)}_{j}\left(\sum_{i=1}^{T}{\bm{g}}^{\sf
dum}_{j}[i]\right)=-\sum_{i=1}^{T}\sum_{j=0}^{t}\lambda^{(t+1)}_{j}{\bm{g}}^{\sf
dum}_{j}[i]$ $\displaystyle=\sum_{i=1}^{T}x_{t+1}[i].$

This completes the proof of (C.10). Next, we will show that

$\displaystyle\delta\cdot g^{y}_{t+1}=\sum_{i=1}^{T}{\bm{g}}^{\sf
dum}_{t+1}[i].$ (C.11)

This follows because

$\displaystyle\delta\cdot g^{y}_{t+1}$ $\displaystyle=\delta(1+\mu
y_{t+1})=\delta+\mu\left(\sum_{i=1}^{T}{\bm{x}}^{\sf dum}_{t+1}[i]\right)$
$\displaystyle\overset{(a)}{=}{\bm{g}}^{\sf
dum}_{t+1}[t+2]+\sum_{i=1}^{t+1}{\bm{g}}^{\sf
dum}_{t+1}[i]=\sum_{i=1}^{T}{\bm{g}}^{\sf dum}_{t+1}[i],$

where ($a$) uses the fact that ${\bm{x}}^{\sf dum}_{t+1}[i]=0$ for all
$i>t+1$. ∎

By Lemma C.4, it holds that

$\displaystyle\delta\cdot y_{T}$ $\displaystyle=\sum_{i=1}^{T}{\bm{x}}^{\sf
dum}_{T}[i]$ (C.12)

On the other hand, in order to achieve $\varepsilon$-suboptimality, we need
$y_{T}^{2}\gtrsim\frac{1}{\mu^{2}}$. Moreover, for
$\varepsilon$-suboptimality, we also need $\sum_{i=1}^{T}{\bm{x}}^{\sf
dum}_{T}[i]^{2}\lesssim\frac{\varepsilon}{\mu}$. Therefore, letting
${\bm{x}}_{T}^{\sf exact}=(({\bm{x}}^{\sf dum}_{T})^{\sf exact},y_{T}^{\sf
exact})$ be the iterate with exact gradients, since $({\bm{x}}^{\sf
dum}_{T})^{\sf exact}=0$, the condition (C.12) implies the following: whenever
$\sum_{i=1}^{T}{\bm{x}}^{\sf dum}_{T}[i]^{2}\lesssim\frac{\varepsilon}{\mu}$,

$\displaystyle\left\lVert{\bm{x}}_{T}-{\bm{x}}_{T}^{\sf
exact}\right\rVert^{2}$ $\displaystyle\geq\left\lVert{\bm{x}}^{\sf
dum}_{T}-({\bm{x}}^{\sf dum}_{T})^{\sf
exact}\right\rVert^{2}=\sum_{i=1}^{T}{\bm{x}}^{\sf dum}_{T}[i]^{2}$
$\displaystyle\geq\frac{1}{T}\left(\sum_{i=1}^{T}{\bm{x}}^{\sf
dum}_{T}[i]\right)^{2}=\frac{1}{T}\cdot\delta^{2}\cdot\left(y_{T}\right)^{2}\gtrsim\frac{\delta^{2}}{T\mu^{2}}.$

This is precisely equal to the desired lower bound.

##### Actual proof for the stochastic noise case.

Now coming back to the stochastic inexact gradient (C.8), one can prove the
following analog of Lemma C.4:

$\displaystyle\begin{split}\delta\cdot|y_{t}|\leq\sum_{i=1}^{T}|{\bm{x}}^{\sf
dum}_{t}[i]|\quad\text{and}\quad\delta\cdot|g^{y}_{t}|\leq\sum_{i=1}^{T}|{\bm{g}}^{\sf
dum}_{t}[i]|\,,\end{split}\quad\text{for each }t=0,1,2,\dots,T.$ (C.13)

Here we note that the construction ensures that $|{\bm{x}}^{\sf
dum}_{T}[i]|$’s are deterministic quantities (because regardless of whether
$r_{t}=\pm 1$ the absolute value is the same), and that is why we do not write
the expectation operators next to them.

The above result holds for the following reason, when $r_{t}=+1$ for all $t$,
the stochastic inexact gradient reduces to the non-stochastic inexact
gradient, in which case the equality holds in (C.13) without absolute values.
With $r_{t}=\pm 1$, one can no longer argue this. On the other hand, one can
apply triangle inequalities to obtain (C.13).

Again, in order to achieve $\varepsilon$-suboptimality, we need
$y_{T}^{2}\gtrsim\frac{1}{\mu^{2}}$ and $\sum_{i=1}^{T}{\bm{x}}^{\sf
dum}_{T}[i]^{2}\lesssim\frac{\varepsilon}{\mu}$. Therefore, whenever
$\sum_{i=1}^{T}{\bm{x}}^{\sf dum}_{T}[i]^{2}\lesssim\frac{\varepsilon}{\mu}$,
we have

$\displaystyle\operatorname*{\mathbb{E}}\left\lVert{\bm{x}}_{T}-\operatorname*{\mathbb{E}}[{\bm{x}}_{T}]\right\rVert^{2}$
$\displaystyle\geq\sum_{i=1}^{T}|{\bm{x}}^{\sf
dum}_{T}[i]|^{2}\geq\frac{1}{T}\left(\sum_{i=1}^{T}|{\bm{x}}^{\sf
dum}_{T}[i]|\right)^{2}\geq\frac{1}{T}\cdot\delta^{2}\cdot|y_{T}|^{2}\gtrsim\frac{\delta^{2}}{T\mu^{2}}.$

This completes the proof. ∎

####  C.5 Non-stochastic inexact gradient model

[] (Lower Bound) Let $\varepsilon>0$ be a small constant, and $T$ be a given
number of iterations. There exists a $O(1)$-smooth $\mu$-strongly convex
function $f:\mathbb{R}^{2}\to\mathbb{R}$ with a non-stochastic inexact
gradient model such that for any FOI algorithm $\mathcal{A}$ that starts at
${\bm{x}}_{0}={\bm{0}}$ has the $(\varepsilon,\delta)$-deviation lower bounded
by $\Omega(\frac{\delta^{2}}{\mu^{2}}\wedge\frac{\varepsilon}{\mu})$.

###### Proof.

For simplicity, we assume throughout the proof that $D=1$ and for
${\bm{x}}=(x,y)$ where $x,y\in\mathbb{R}$, consider the cost

$\displaystyle f(x,y):=y+\frac{\mu}{2}y^{2}+\frac{\mu}{2}x^{2}$

We consider the initialization ${\bm{x}}_{0}=(0,0)$. Next, consider the
following non-stochastic inexact oracle

$\displaystyle g({\bm{x}}_{t})=\nabla f({\bm{x}}_{t})+\delta{\bm{e}}_{1}$
(C.14)

where ${\bm{e}}_{1}$ is the first coordinate vector. Then from this
construction, one can verify similarly to Subsection C.4 that

$\displaystyle g({\bm{x}}_{t})[1]=\delta\cdot
g({\bm{x}}_{t})[2]\quad\text{and}\quad x_{t}=\delta\cdot y_{t}$

for $t=0,1,\dots,T$.

Now from the $\varepsilon$-suboptimality, it must be that
$y_{T}^{2}\gtrsim\frac{1}{\mu^{2}}$ and
$x_{T}^{2}\lesssim\frac{\varepsilon}{\mu}$. Hence, whenever
$x_{T}^{2}\lesssim\frac{\varepsilon}{\mu}$ holds, we have the following
deviation bound since $x_{T}^{\sf exact}=0$:

$\displaystyle\left\lVert{\bm{x}}_{T}-{\bm{x}}_{T}^{\sf
exact}\right\rVert^{2}$ $\displaystyle\geq|x_{T}-x_{T}^{\sf
exact}|^{2}=\delta^{2}\cdot y_{t}^{2}\gtrsim\frac{\delta^{2}}{\mu^{2}}.$

This completes the proof. ∎

####  C.6 Inexact initialization model (strongly convex costs)

[] (Lower Bound) Let $\varepsilon>0$ be a small constant, and $T$ be a given
number of iterations. There exists a $O(1)$-smooth $\mu$-strongly convex
function $f:\mathbb{R}^{\Omega(T)}\to\mathbb{R}$ such that for any FOI
algorithm $\mathcal{A}$ the $(\varepsilon,\delta)$-deviation lower bounded by
$\Omega(\exp(-\Omega(T))\delta^{2}\wedge\frac{\varepsilon}{\mu})$ w.r.t. the
reference point ${\bm{x}}_{0}^{\sf ref}={\bm{0}}$.

###### Proof.

We use the construction in [Nesterov, 2018, Theorem 2.1.13]. In particular,
for simplicity, we consider the construction for the infinite dimensional
Hilbert space $\ell_{2}$ as it simplifies the proof; in fact, a similar
argument works for $\mathbb{R}^{\Omega(T)}$. Let us recall the construction
(we follow the presentation in [Bubeck, 2014, Theorem 3.15]). Let
$A:\ell_{2}\rightarrow\ell_{2}$ be the linear operator that corresponds to the
infinite tri-diagonal matrix with $2$ on the diagonal and $-1$ on the upper
and lower diagonals. For some constant $\kappa\geq 1$, consider the following
$\mu$-strongly convex cost:

$\displaystyle f^{\sf Nes}(x)=\frac{\mu(\kappa-1)}{8}\left(\langle
Ax,x\rangle-2\langle{\bm{e}}_{1},x\rangle\right)+\frac{\mu}{2}\left\lVert
x\right\rVert^{2}\quad\text{and}\quad
q:=\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}.$

For the zero initialization $x_{0}=(0,0,\dots)\in\ell_{2}$, the cost satisfies
the following properties (see the proof of [Nesterov, 2018, Theorem 2.1.13]):

  * •

Output of any FOI satisfies $x_{t}[i]=0,\forall i\geq t$.

  * •

$x^{*}[i]=q^{i}$.

  * •

$\left\lVert
x_{0}-x^{*}\right\rVert^{2}=\sum_{i=1}^{\infty}(x^{*}[i])^{2}=\sum_{i=1}^{\infty}q^{2i}=\frac{q^{2}}{1-q^{2}}$.

  * •

$\left\lVert
x_{t}-x^{*}\right\rVert^{2}\geq\sum_{i=k+1}^{\infty}q^{2i}=\frac{q^{2(t+1)}}{1-q^{2}}=q^{2t}\left\lVert
x_{0}-x^{*}\right\rVert^{2}$.

Now we consider the following cost function: for
${\bm{x}}=(x,y)\in\ell_{2}\times\ell_{2}$

$\displaystyle f({\bm{x}})=f^{\sf Nes}(x)+f^{\sf Nes}(y)\,,$

and we consider the two initializations:

$\displaystyle{\bm{x}}_{0}^{\sf ref}$
$\displaystyle=\Big{(}(0,0,0,\dots),(q,q^{2},q^{3},\dots)\Big{)}$
$\displaystyle{\bm{x}}_{0}$
$\displaystyle=\Big{(}(0,0,0,\dots),(q,q^{2},q^{3},\dots,q^{L},0,0,\dots)\Big{)}$

for $L=\Omega(\log(1/\delta))$ is chosen such that $q^{L+1}/(1-q)\leq\delta$.
Then, it follow that $\left\lVert{\bm{x}}_{0}-{\bm{x}}_{0}^{\sf
ref}\right\rVert=\frac{q^{L+1}}{1-q}\leq\delta$. On the other hand, it follows
from the above property that

$\displaystyle\left\lVert{\bm{x}}_{t}-{\bm{x}}_{t}^{\sf
ref}\right\rVert^{2}\geq q^{t}\frac{q^{2L+1}}{1-q^{2}}\approx
q^{t}\cdot\delta^{2}.$

Hence, as long as
$\left\lVert{\bm{x}}_{t}-{\bm{x}}^{*}\right\rVert^{2},\left\lVert{\bm{x}}_{t}^{\sf
ref}-{\bm{x}}^{*}\right\rVert^{2}\lesssim\frac{\varepsilon}{\mu}$, the
deviation lower bound follows. ∎

###  Appendix D Proof of lower bounds (nonsmooth costs)

In this section, we present the proofs of lower bounds for nonsmooth costs.
The proof will be based on more complicated constructions than those for the
case of smooth costs, so before we dive into the proofs, we first present some
intuition behind the constructions.

####  D.1 Warm-up: lower bound against GD

In this section, as a warm-up, we will prove a (weaker) lower bound for a
simplified setting. In particular, we prove the lower bound against gradient
descent (GD). Formally, in the definition of FOI, we restrict that
$\lambda^{(t)}_{i}\equiv\lambda_{i}$ (i.e., the coefficient is a positive
number does not depend on the iterations $t$). In other words, for
$\lambda_{i}$, $i=0,1,\dots,T-1$,

$\displaystyle{\bm{x}}_{t}={\bm{x}}_{0}-\sum_{i=0}^{t-1}\lambda_{i}g({\bm{x}}_{i})\quad\text{for
each $t=1,2,\dots,T$}.$ (D.1)

Note that this is precisely GD with step sizes $\lambda_{t}$’s. For the lower
bound construction, let ${\bm{x}}=({\bm{x}}^{\sf
err},w)\in\mathbb{R}^{T}\times\mathbb{R}$ and consider the cost

$\displaystyle f({\bm{x}}^{\sf
err},w)=\underbrace{\max\left\\{0,~{}\max_{i=1,\dots,T}\left\\{{\bm{x}}^{\sf
err}[i]\right\\}\right\\}}_{=:G({\bm{x}}^{\sf
err})}+\underbrace{2\varepsilon\cdot\max\\{w+1,0\\}}_{=:\ell(w)}\,.$

Since the above cost function is nonsmooth, we specify the subgradient oracle
as follows: for both max terms above, we consider the subgradient oracle that
outputs the subgradient corresponding to the first argument that achieves the
maximum. Note that $\inf_{{\bm{x}}}f({\bm{x}})=0$. Consider the zero
initialization, i.e., $({\bm{x}}^{\sf err}_{0},w_{0})=({\bm{0}},0)$, and we
write the iterates as ${\bm{x}}_{t}:=({\bm{x}}^{\sf
err}_{t},w_{t})\in\mathbb{R}^{T}\times\mathbb{R}$.

For intuition, we describe the role of each coordinate:

  * $\blacksquare$

The first $T$ coordinates, ${\bm{x}}^{\sf err}\in\mathbb{R}^{T}$, correspond
to the part where the errors due to inexact oracle are added.

  * $\blacksquare$

The last coordinate, $w\in\mathbb{R}$, governs the overall cost; in order to
achieve $\varepsilon$-accuracy, the optimization algorithm has to decrease
coordinate $w$ by at least $1/2$.

The proof proceeds by considering two different scenarios:

##### Scenario 1 (exact gradients).

Consider the case where there is no noise in the gradients, i.e.,
$g({\bm{x}}_{t})=\nabla f({\bm{x}}_{t})$ for all $t$. Then, since
${\bm{x}}^{\sf err}_{0}=\bm{0}$, it follows that $\nabla G({\bm{x}}^{\sf
err}_{t})=\bm{0}$ for all $t$. Hence the algorithm will only update coordinate
$w_{t}$. Note that $\ell(w_{0})=2\varepsilon$, and hence in order to achieve
$\ell(w_{T})\leq\varepsilon$, it must be that $w_{T}\leq-1/2$. On the other
hand, we have

$\displaystyle\frac{\partial}{\partial
w}\ell(w)=0~{}~{}\text{or}~{}~{}2\varepsilon\quad\text{for
any}~{}~{}w\in\mathbb{R}.$

Hence, in order to achieve $w_{T}\leq-1/2$, it must be that

$\displaystyle\sum_{t=0}^{T-1}\lambda_{t}\geq\Omega(1/\varepsilon).$ (D.2)

This condition is analogous to (C.4) from the lower bound proof for smooth
costs (Subsection C.1).

##### Scenario 2 (inexact gradients).

Now let us consider the case where the gradient error during the $t$-th
iteration is non-stochastic and equal to $-\delta{\bm{e}}_{t}$, i.e.,

$\displaystyle g({\bm{x}}_{t})=\nabla
f({\bm{x}}_{t})-\delta{\bm{e}}_{t+1}\quad\text{for all $t=0,1,\dots,T-1$}.$

Here ${\bm{e}}_{t}$ denotes the $t$-th coordinate vector. Let us assume that
$\delta$ is much smaller than all the step sizes $\lambda_{t}$, in particular,
such that $\lambda_{i}\delta\ll\lambda_{i+1}$ for all $i=0,\dots,T-2$. Then
from GD iterations defined as (D.1), one can deduce that

$\displaystyle{\bm{x}}^{\sf err}_{t}=(-\lambda_{1}+\lambda_{0}\delta,\
-\lambda_{2}+\lambda_{1}\delta,\ \cdots,\ -\lambda_{t-1}+\lambda_{t-2}\delta,\
+\lambda_{t}\delta,0,\dots,0)\,.$

Thus, the following estimate on the deviation holds:

$\displaystyle\left\lVert{\bm{x}}^{\sf
err}_{T}\right\rVert^{2}=\left\lVert\sum_{t=1}^{T-1}(\lambda_{t}-\lambda_{t-1}\delta){\bm{e}}_{t}+\lambda_{T-1}\delta{\bm{e}}_{T}\right\rVert^{2}=\sum_{t=1}^{T-1}(\lambda_{t}-\lambda_{t-1}\delta)^{2}\approx\sum_{t=1}^{T-1}\lambda_{t}^{2}\,.$
(D.3)

##### Combining the two scenarios.

Thus far, we have obtained (D.2) and (D.3) from the two different scenarios.
The condition (D.2) shows that in order to achieve
$\varepsilon$-suboptimality, stepsizes have to add up to a large number, more
precisely, $\sum_{t=0}^{T-1}\lambda_{t}=\Omega(1/\varepsilon)$. On the other
hand, (D.3) characterizes that the deviation is on the order of the quantity
$\sum_{t=1}^{T-1}\lambda_{t}^{2}$. In order to formally connect these two
conditions, we make the following assumption:

$\displaystyle\lambda_{0}\leq O\left(\sum_{t=1}^{T-1}\lambda_{t}\right)\,.$
(D.4)

Then with this assumption, one obtain the following deviation bound:

$\displaystyle\left\lVert{\bm{x}}^{\sf
err}_{T}\right\rVert^{2}\approx\sum_{t=1}^{T-1}(\lambda^{(T)}_{t})^{2}\overset{(a)}{\geq}\frac{1}{T-1}\cdot\left(\sum_{t=1}^{T-1}\lambda^{(T)}_{t}\right)^{2}\gtrsim\frac{1}{T-1}\cdot\left(\sum_{t=0}^{T-1}\lambda^{(T)}_{t}\right)^{2}\gtrsim\frac{1}{T\varepsilon^{2}}\,,$

where $(a)$ is due to the Cauchy-Schwartz inequality. This is precisely the
desired lower bound.

For the lower bound against the entire class of FOI, there are some other
technical challenges arising from the fact that the coefficients
$\lambda^{(t)}_{i}$’s not only depend on $t$, but also could take negative
values. We need a more elaborate lower bound construction, as we explain in
the subsequent subsections.

####  D.2 Helper function

Before actual proofs, we introduce a helper function that we will use
throughout the proofs of FOI lower bounds. Let $\chi:\mathbb{R}\to\mathbb{R}$
is a non-smooth convex function defined as $\chi(x):=\max\\{x,0\\}$ and the
subgradients are defined as

$\displaystyle\begin{split}\nabla_{x}\chi(x)&:=\begin{cases}+1,&\text{if
}x\geq 0,\\\ 0,&\text{if }x<y.\end{cases}\end{split}$ (D.5)

The choice of subgradient $+1$ at the origin will play a crucial role in the
later proofs. For ${\bm{x}},{\bm{y}},{\bm{z}}\in\mathbb{R}^{T}$ and
$v\in\mathbb{R}$, let
$\mathcal{G}:({\bm{x}},{\bm{y}},{\bm{z}})\in\mathbb{R}^{3T}\to\mathbb{R}$ be
defined as

$\displaystyle\mathcal{G}({\bm{x}},{\bm{y}},{\bm{z}})$
$\displaystyle:=\max\\{0,~{}\mathcal{K}({\bm{x}},{\bm{y}},{\bm{z}})\\}\quad\text{and}$
(D.6) $\displaystyle\mathcal{K}({\bm{x}},{\bm{y}},{\bm{z}})$
$\displaystyle:=\max_{i=1,\dots,T}\Big{\\{}\chi({\bm{y}}[i])+\sum_{j=1}^{i-1}\frac{\left|{\bm{x}}[j]\right|}{2^{j-1}}+\frac{{\bm{x}}[i]}{2^{i-1}},~{}~{}\chi({\bm{z}}[i])+\sum_{j=1}^{i-1}\frac{\left|{\bm{x}}[j]\right|}{2^{j-1}}-\frac{{\bm{x}}[i]}{2^{i-1}}\Big{\\}}\,.$
(D.7)

Then $\mathcal{G}$ is clearly convex, as it is the maximum of convex
functions.

We specify the subgradients of $\mathcal{G}$ as follows: for all max terms in
(D.8), we get the subgradient of the first argument that achieves the maximum.
Then $\mathcal{G}$ is $O(1)$-Lipschitz: for any
${\bm{x}},{\bm{y}},{\bm{z}}\in\mathbb{R}^{T}$,

$\displaystyle\left\lVert\nabla\mathcal{G}({\bm{x}},{\bm{y}},{\bm{z}})\right\rVert^{2}\leq
1+\sum_{j=1}^{T}(\frac{1}{2^{j-1}})^{2}\leq
1+\sum_{j=1}^{\infty}\frac{1}{4^{j-1}}\leq 1+4/3\,.$

####  D.3 Stochastic inexact gradient model

[] (Lower Bound) Let $\varepsilon>0$ and $T$ be a given number of iterations.
There exists a $O(1)$-Lipschitz (nonsmooth) convex function
$f:\mathbb{R}^{O(T)}\to\mathbb{R}$ with a stochastic inexact gradient model
such that any FOI algorithm $\mathcal{A}$ that satisfies
$|\lambda_{0}^{(T)}|\leq
O\left(\left|\sum_{t=1}^{T-1}\lambda^{(T)}_{t}\right|\right)$ and starts at
${\bm{x}}_{0}={\bm{0}}$ has its $(\varepsilon,\delta)$-deviation lower bounded
by $\Omega(\frac{1}{T\varepsilon^{2}})$.

###### Proof.

For ${\bm{x}}=({\bm{x}}^{\sf err},{\bm{y}},{\bm{z}},w)$ where ${\bm{x}}^{\sf
err},{\bm{y}},{\bm{z}}\in\mathbb{R}^{T}$ and $w\in\mathbb{R}$, consider the
cost

$\displaystyle f({\bm{x}}^{\sf
err},{\bm{y}},{\bm{z}},w)=\mathcal{G}({\bm{x}}^{\sf
err},{\bm{y}},{\bm{z}})+\underbrace{2\varepsilon\cdot\max\\{w+1,0\\}}_{=:\ell(w)}\,,$
(D.8)

where $\mathcal{G}$ is defined in (D.6). Then, $f$ is convex since both
$\mathcal{G}$ and $\ell$ are convex, and $O(1)$-Lipschitz since both
$\mathcal{G}$ and $\ell$ are $O(1)$-Lipschitz.

For intuition, we describe the role of each coordinate as we did in the warm-
up section (Subsection D.1):

  * $\blacksquare$

The first $T$ coordinates, ${\bm{x}}^{\sf err}\in\mathbb{R}^{T}$, correspond
to the part where the errors due to inexact oracle are added.

  * $\blacksquare$

The next $2T$ coordinates, ${\bm{y}},{\bm{z}}\in\mathbb{R}^{T}$ will
contribute to large deviation when there are errors in the gradients.

  * $\blacksquare$

The last coordinate, $w\in\mathbb{R}$, governs the overall cost; in order to
achieve $\varepsilon$-accuracy, the optimization algorithm has to decrease
coordinate $w$ by at least $1/2$.

We use the following notation throughout the proof:
${\bm{x}}_{t}=({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t},w_{t})\in(\mathbb{R}^{T})^{3}\times\mathbb{R}$.
Consider the zero initialization ${\bm{x}}_{0}=({\bm{0}},{\bm{0}},{\bm{0}},0)$
and the following inexact gradient error for $t=0,1,2,\dots,T-1$:

$\displaystyle g({\bm{x}}_{t})=\nabla f({\bm{x}}_{t})+\delta
r_{t}\cdot{\bm{e}}_{1+t}\quad\text{for
$r_{t}\overset{iid}{\sim}\operatorname{Unif}\\{\pm 1\\}$},$ (D.9)

where ${\bm{e}}_{j}$ is the $j$-th coordinate vector. The following lemma
characterizes the key feature of the above construction. Note that Subsection
D.3 is the place where we use the following additional assumption that we made
for the case of nonsmooth costs: “ _for the case of nonsmooth costs, we
additionally assume that the coefficient of the latest gradient is nonzero,
i.e., $\lambda^{(t)}_{t-1}\neq 0$ for all $t$._”

Under the inexact gradient (D.9), the subgradient $\nabla\mathcal{G}$ has the
following properties:

  * $\blacksquare$

For each $t=1,2,\dots,T-1$, there exists $i_{t}\in\\{1,\dots,t\\}$ such that
the following holds:

$\displaystyle\begin{cases}\frac{\partial}{\partial{\bm{y}}[i_{t}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=1,~{}~{}\frac{\partial}{\partial{\bm{z}}[i_{t}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=0&\text{with probability }$1/2$,\\\
\frac{\partial}{\partial{\bm{y}}[i_{t}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=0,~{}~{}\frac{\partial}{\partial{\bm{z}}[i_{t}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=1&\text{with probability
}$1/2$.\end{cases}$ (D.10)

Moreover, for $i\neq i_{t}$,
$\frac{\partial}{\partial{\bm{y}}[i]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=\frac{\partial}{\partial{\bm{z}}[i]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=0$.

  * $\blacksquare$

If $i_{t}\neq t$ (i.e., $i_{t}<t$), then $i_{t}=i_{t^{\prime}}$ for some
$t^{\prime}<t$, and it holds that
$\frac{\partial}{\partial{\bm{y}}[i_{t}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=\frac{\partial}{\partial{\bm{y}}[i_{t^{\prime}}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t^{\prime}},{\bm{y}}_{t^{\prime}},{\bm{z}}_{t^{\prime}})$ and
$\frac{\partial}{\partial{\bm{z}}[i_{t}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=\frac{\partial}{\partial{\bm{z}}[i_{t^{\prime}}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t^{\prime}},{\bm{y}}_{t^{\prime}},{\bm{z}}_{t^{\prime}})$.

###### Proof of Subsection D.3.

Let us recall the definition of $\mathcal{G}({\bm{x}},{\bm{y}},{\bm{z}})$:

$\displaystyle\max\left\\{0~{}~{},\max_{i=1,\dots,T}\Big{\\{}\chi({\bm{y}}[i])+\sum_{j=1}^{i-1}\frac{\left|{\bm{x}}[j]\right|}{2^{j-1}}+\frac{{\bm{x}}[i]}{2^{i-1}},~{}~{}\chi({\bm{z}}[i])+\sum_{j=1}^{i-1}\frac{\left|{\bm{x}}[j]\right|}{2^{j-1}}-\frac{{\bm{x}}[i]}{2^{i-1}}\Big{\\}}\right\\}\,.$
(D.11)

From this, it is clear that there must be at most one $i\in\\{1,\dots,T\\}$
for which either
$\frac{\partial}{\partial{\bm{y}}[i]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})\neq 0$ or
$\frac{\partial}{\partial{\bm{z}}[i]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})\neq 0$. This proves the “ _Moreover, for
$i\neq i_{t}$, $\frac{\partial}{\partial{\bm{y}}[i]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=\frac{\partial}{\partial{\bm{z}}[i]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=0$_” part of the first bullet point.

Next, we prove the expression (D.10). We begin with $t=1$. Since
$g({\bm{x}}_{0})=\nabla f({\bm{x}}_{0})+\delta r_{0}\cdot{\bm{e}}_{1}$ and
$\lambda^{(1)}_{0}\neq 0$, it follows that ${\bm{x}}^{\sf err}_{1}[1]\neq 0$.
Then, we claim that the first bullet point holds for $t=1$. Since we know
${\bm{x}}^{\sf err}_{1}[1]\neq 0$ and ${\bm{x}}^{\sf
err}_{1}[2],\dots,{\bm{x}}^{\sf err}_{1}[T]=0$, for ${\bm{x}}_{1}$, the
maximum in (D.11) is achieved by $i=1$. This implies that $i_{1}=1$. Moreover,
depending on the sign of $r_{0}$, we either have
$\frac{\partial}{\partial{\bm{y}}[1]}\mathcal{G}({\bm{x}}^{\sf
err}_{1},{\bm{y}}_{1},{\bm{z}}_{1})=1$ or
$\frac{\partial}{\partial{\bm{z}}[1]}\mathcal{G}({\bm{x}}^{\sf
err}_{1},{\bm{y}}_{1},{\bm{z}}_{1})=1$ with equal probability. Thus, (D.10)
holds for $t=1$.

Next, consider $t>1$. Since $g({\bm{x}}_{t-1})=\nabla f({\bm{x}}_{t-1})+\delta
r_{t-1}\cdot{\bm{e}}_{t}$ and $\lambda^{(t)}_{t-1}\neq 0$, it follows that
${\bm{x}}^{\sf err}_{t}[t]\neq 0$. Moreover, we know ${\bm{x}}^{\sf
err}_{t}[t+1],\dots,{\bm{x}}^{\sf err}_{t}[T]=0$. Hence, we have the following
two scenarios:

  * $\bullet$

Case 1: ${\bm{y}}[i],{\bm{z}}[i]\leq 0$ for all $i\in\\{1,2,\dots,t-1\\}$.
Note that this hold—for instance—if the coefficients FOI are all non-negative,
i.e., $\lambda^{(t)}_{i}\geq 0$ of for all $i\in\\{1,2,\dots,t-1\\}$ (most
first order optimization algorithms usually follow this). In that case, the
maximum in (D.11) is achieved by $i=t$ This implies that $i_{t}=t$. Moreover,
depending on the sign of $r_{t-1}$, we either have
$\frac{\partial}{\partial{\bm{y}}[t]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=1$ or
$\frac{\partial}{\partial{\bm{z}}[t]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=1$ with equal probability. Thus, again
(D.10) holds for $t$.

  * $\bullet$

Case 2: Somehow FOI chooses to follow positive gradient directions (which is
unlikely in practice) and it happens that ${\bm{y}}[i]>0$ or ${\bm{z}}[i]>0$
for some $i\in\\{1,2,\dots,t-1\\}$. In such a case, the maximum in (D.11)
could be achieved by $i\in\\{1,2,\dots,t-1\\}$, i.e.,
$i_{t}\in\\{1,2,\dots,t-1\\}$. Then it must be that ${\bm{y}}[i_{t}]>0$ or
${\bm{z}}[i_{t}]>0$. This can happen only if $i_{t}=i_{t^{\prime}}$ for some
$t^{\prime}<t$. Hence, it follows that
$\frac{\partial}{\partial{\bm{y}}[i_{t}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=\frac{\partial}{\partial{\bm{y}}[i_{t^{\prime}}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t^{\prime}},{\bm{y}}_{t^{\prime}},{\bm{z}}_{t^{\prime}})$ and
$\frac{\partial}{\partial{\bm{z}}[i_{t}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=\frac{\partial}{\partial{\bm{z}}[i_{t^{\prime}}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t^{\prime}},{\bm{y}}_{t^{\prime}},{\bm{z}}_{t^{\prime}})$. This proves
the second bullet point in the statement. In particular, (D.10) holds for $t$.

This completes the proof of Subsection D.3. ∎

Now we use Subsection D.3 to prove Subsection D.3. From the construction
(D.8), we know that $f({\bm{x}}_{0})=2\varepsilon$. In order to achieve
$\varepsilon$-accuracy, we need $w_{T}\leq-1/2$. Note that
$\frac{\partial}{\partial w}\ell(w_{t})=2\varepsilon$ for all $t\geq 0$. From
the fact that $w_{T}\leq-\frac{1}{2}$, it follows that

$\displaystyle\sum_{t=0}^{T-1}\lambda^{(T)}_{t}\frac{\partial}{\partial
w}\ell(v_{t},w_{t})\geq\frac{1}{2}\quad\Longleftrightarrow\quad\sum_{t=0}^{T-1}\lambda^{(T)}_{t}\geq\frac{1}{4\varepsilon}\,.$

Now using the assumption that $|\lambda_{0}^{(T)}|\leq
O\left(\left|\sum_{t=1}^{T-1}\lambda^{(T)}_{t}\right|\right)$, we obtain

$\displaystyle\left|\sum_{t=1}^{T-1}\lambda^{(T)}_{t}\right|\gtrsim|\lambda^{(T)}_{0}|+\left|\sum_{t=1}^{T-1}\lambda^{(T)}_{t}\right|\geq\left|\sum_{t=0}^{T-1}\lambda^{(T)}_{t}\right|\geq\frac{1}{4\varepsilon}\,,$

which leads to the following conditon:

$\displaystyle\left|\sum_{t=1}^{T-1}\lambda^{(T)}_{t}\right|\geq\Omega\left(\frac{1}{\varepsilon}\right)$
(D.12)

This condition is analogous to (D.2) from Subsection D.1. Now to better
illustrate our proof strategy for the remaining part, we first consider a
special case.

##### Warm-up: proof for the special case.

As a warm-up, we first consider the special case where in the definition of
FOI, all the coefficients $\lambda^{(t)}_{i}$ are non-negative, i.e.,

$\displaystyle{\bm{x}}_{t}={\bm{x}}_{0}-\sum_{i=0}^{t-1}\lambda^{(t)}_{i}g({\bm{x}}_{i})\quad\text{for
some }\lambda_{i}^{(t)}\geq 0,~{}i=0,\dots,t-1,$ (D.13)

Then this case belongs to Case 1 in the proof of Subsection D.3. As a
consequence, $i_{t}=t$ in Subsection D.3 and the following conclusion holds:

  * $\blacksquare$

For each $t=1,2,\dots,T-1$,

$\displaystyle\begin{cases}\frac{\partial}{\partial{\bm{y}}[t]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=1,~{}~{}\frac{\partial}{\partial{\bm{z}}[t]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=0&\text{with probability }$1/2$,\\\
\frac{\partial}{\partial{\bm{y}}[t]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=0,~{}~{}\frac{\partial}{\partial{\bm{z}}[t]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=1&\text{with probability
}$1/2$.\end{cases}$ (D.14)

Moreover, for $i\neq t$,
$\frac{\partial}{\partial{\bm{y}}[i]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=\frac{\partial}{\partial{\bm{z}}[i]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=0$.

We use (D.14) to lower bound
$\operatorname*{\mathbb{E}}\left\lVert{\bm{x}}_{T}-\operatorname*{\mathbb{E}}[{\bm{x}}_{T}]\right\rVert^{2}$.
From (D.14), it holds that for $t=1,2,\dots,T$,

$\displaystyle{\bm{y}}_{T}[t]=\begin{cases}\lambda^{(T)}_{t},&\text{with
probability}~{}$1/2$,\\\ 0,&\text{with probability}~{}$1/2$\,.\end{cases}$

Hence, we have

$\displaystyle\operatorname*{\mathbb{E}}\left\lVert{\bm{y}}_{T}-\operatorname*{\mathbb{E}}[{\bm{y}}_{T}]\right\rVert^{2}$
$\displaystyle\geq\sum_{t=1}^{T-1}\operatorname*{\mathbb{E}}\left({\bm{y}}_{T}[t]-\operatorname*{\mathbb{E}}{\bm{y}}_{T}[t]\right)^{2}=\sum_{t=1}^{T-1}\operatorname*{\mathbb{E}}\left({\bm{y}}_{T}[t]-\frac{1}{2}\lambda^{(T)}_{t}\right)^{2}=\frac{1}{4}\sum_{t=1}^{T-1}(\lambda^{(T)}_{t})^{2}$
$\displaystyle\overset{(a)}{\geq}\frac{1}{4}\cdot\frac{1}{T-1}\cdot\left(\sum_{t=1}^{T-1}\lambda^{(T)}_{t}\right)^{2}\overset{\eqref{large_coeff_nonsmooth}}{\gtrsim}\frac{1}{4}\cdot\frac{1}{T-1}\cdot\left(\frac{1}{\varepsilon}\right)^{2}\gtrsim\frac{1}{T\varepsilon^{2}}.$

Here $(a)$ follows from Cauchy-Schwarz inequality. Therefore, we get the
desired deviation bound as follows:

$\displaystyle\operatorname*{\mathbb{E}}\left\lVert{\bm{x}}_{T}-\operatorname*{\mathbb{E}}[{\bm{x}}_{T}]\right\rVert^{2}$
$\displaystyle\geq\operatorname*{\mathbb{E}}\left\lVert{\bm{y}}_{T}-\operatorname*{\mathbb{E}}[{\bm{y}}_{T}]\right\rVert^{2}\gtrsim\frac{1}{T\varepsilon^{2}}.$

##### The proof for the general case.

Now we consider the case of general FOI where the coefficients
$\lambda^{(t)}_{i}$ are not necessarily non-negative. With $i_{t}$ defined in
the statement of Subsection D.3, let

$\displaystyle\mathcal{I}:=\left\\{i_{t}\in\\{1,2,\dots,T\\}~{}:~{}t=1,2,\dots,T-1\right\\}\,.$

Then Subsection D.3 ensures that that for each $i\in\mathcal{I}$,

$\displaystyle{\bm{y}}_{T}[i]=\begin{cases}\sum_{\begin{subarray}{c}t=1,2,\dots,T-1\\\
~{}\text{s.t.}~{}i_{t}=i\end{subarray}}\lambda^{(T)}_{t},&\text{with
probability}~{}$1/2$,\\\ 0,&\text{with probability}~{}$1/2$\,.\end{cases}$

Hence, we have

$\displaystyle\operatorname*{\mathbb{E}}\left\lVert{\bm{y}}_{T}-\operatorname*{\mathbb{E}}[{\bm{y}}_{T}]\right\rVert^{2}$
$\displaystyle\geq\sum_{t=1}^{T-1}\operatorname*{\mathbb{E}}\left({\bm{y}}_{T}[t]-\operatorname*{\mathbb{E}}{\bm{y}}_{T}[t]\right)^{2}=\frac{1}{4}\sum_{i\in\mathcal{I}}\left(\sum_{\begin{subarray}{c}t=1,2,\dots,T-1\\\
~{}\text{s.t.}~{}i_{t}=i\end{subarray}}\lambda^{(T)}_{t}\right)^{2}$
$\displaystyle\overset{(a)}{\geq}\frac{1}{4}\cdot\frac{1}{|\mathcal{I}|}\cdot\left(\sum_{i\in\mathcal{I}}\sum_{\begin{subarray}{c}t=1,2,\dots,T-1\\\
~{}\text{s.t.}~{}i_{t}=i\end{subarray}}\lambda^{(T)}_{t}\right)^{2}=\frac{1}{4}\cdot\frac{1}{|\mathcal{I}|}\cdot\left(\sum_{t=1}^{T-1}\lambda^{(T)}_{t}\right)^{2}\overset{\eqref{large_coeff_nonsmooth}}{\gtrsim}\frac{1}{T\varepsilon^{2}}.$

Here $(a)$ follows from Cauchy-Schwarz inequality. Therefore, we get the
desired deviation bound as follows:

$\displaystyle\operatorname*{\mathbb{E}}\left\lVert{\bm{x}}_{T}-\operatorname*{\mathbb{E}}[{\bm{x}}_{T}]\right\rVert^{2}\geq\operatorname*{\mathbb{E}}\left\lVert{\bm{y}}_{T}-\operatorname*{\mathbb{E}}[{\bm{y}}_{T}]\right\rVert^{2}\gtrsim\frac{1}{T\varepsilon^{2}}.$

∎

####  D.4 Non-stochastic inexact gradient model

[] (Lower Bound) Let $\varepsilon>0$, and $T$ be a given number of iterations.
There exists a $O(1)$-Lipschitz and nonsmooth convex function
$f:\mathbb{R}^{O(T)}\to\mathbb{R}$ with a non-stochastic inexact gradient
oracle such that for any FOI algorithm $\mathcal{A}$ that satisfies
$|\lambda_{0}^{(T)}|\leq
O\left(\left|\sum_{t=1}^{T-1}\lambda^{(T)}_{t}\right|\right)$ and starts at
${\bm{x}}_{0}={\bm{0}}$ has a minimum of $(\varepsilon,\delta)$-deviation of
$\Omega(\frac{1}{T\varepsilon^{2}}+\frac{\delta^{2}}{\varepsilon^{2}})$.

###### Proof.

We consider an almost identical construction to the one considered in the
proof of Subsection D.3, namely (D.8). The only difference is that now we add
an extra dummy coordinate, namely the $(3T+2)$-th coordinate, which does not
appear in the cost. Let us denote this dummy coordinate by $u$. Concretely, we
consider the following cost: For ${\bm{x}}=({\bm{x}}^{\sf
err},{\bm{y}},{\bm{z}},w,u)$ where ${\bm{x}}^{\sf
err},{\bm{y}},{\bm{z}}\in\mathbb{R}^{T}$ and $w,u\in\mathbb{R}$, consider the
cost

$\displaystyle f({\bm{x}}^{\sf
err},{\bm{y}},{\bm{z}},w,u)=\mathcal{G}({\bm{x}}^{\sf
err},{\bm{y}},{\bm{z}})+\underbrace{2\varepsilon\cdot\max\\{w+1,0\\}}_{=:\ell(w)}\,.$
(D.15)

We denote the iterates of FOI due to _exact_ gradients by

${\bm{x}}_{t}^{\sf exact}=(({\bm{x}}^{\sf err}_{t})^{\sf
exact},{\bm{y}}_{t}^{\sf exact},{\bm{z}}_{t}^{\sf exact},w_{t}^{\sf
exact},u_{t}^{\sf exact})\,.$

Now we define the inexact gradient oracle. The noise in the inexact gradient
oracle consists of two parts. For $t=0,1,2,\dots,T-1$,

$\displaystyle g({\bm{x}}_{t})=\nabla
f({\bm{x}}_{t})+\frac{\delta}{\sqrt{2}}{\bm{e}}_{1+t}+\frac{\delta}{\sqrt{2}}{\bm{e}}_{3T+2}\,,$
(D.16)

where ${\bm{e}}_{j}$ is the $j$-th coordinate vector. Note that the first part
of the error is similar to the error for the stochastic error case, and the
second part of the error is added to the dummy coordinate $u$.

Then analogous to Subsection D.3, one can establish the following result. We
skip the proof since it is very analogous to that of Subsection D.3. Under the
inexact gradient (D.9), the subgradient $\nabla\mathcal{G}$ has the following
properties:

  * $\blacksquare$

For each $t=1,2,\dots,T-1$, there exists $i_{t}\in\\{1,\dots,t\\}$ such that
either one of the following holds:

$\displaystyle\begin{cases}\frac{\partial}{\partial{\bm{y}}[i_{t}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=1,~{}~{}\frac{\partial}{\partial{\bm{z}}[i_{t}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=0,\quad\text{or}\\\
\frac{\partial}{\partial{\bm{y}}[i_{t}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=0,~{}~{}\frac{\partial}{\partial{\bm{z}}[i_{t}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=1.\end{cases}$

Moreover, for $i\neq i_{t}$,
$\frac{\partial}{\partial{\bm{y}}[i]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=\frac{\partial}{\partial{\bm{z}}[i]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=0$.

  * $\blacksquare$

If $i_{t}\neq t$ (i.e., $i_{t}<t$), then $i_{t}=i_{t^{\prime}}$ for some
$t^{\prime}<t$, and it holds that
$\frac{\partial}{\partial{\bm{y}}[i_{t}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=\frac{\partial}{\partial{\bm{y}}[i_{t^{\prime}}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t^{\prime}},{\bm{y}}_{t^{\prime}},{\bm{z}}_{t^{\prime}})$ and
$\frac{\partial}{\partial{\bm{z}}[i_{t}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})=\frac{\partial}{\partial{\bm{z}}[i_{t^{\prime}}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t^{\prime}},{\bm{y}}_{t^{\prime}},{\bm{z}}_{t^{\prime}})$.

The rest of the proof is similar to that of Subsection D.3. From the
construction (D.15), we know that $f({\bm{x}}_{0})=2\varepsilon$. In order to
achieve $\varepsilon$-accuracy, one can similarly deduce that (D.12) holds.

With $i_{t}$ defined in the statement of Subsection D.4, let

$\displaystyle\mathcal{I}:=\left\\{i_{t}\in\\{1,2,\dots,T\\}~{}:~{}t\in\mathcal{T}\right\\}\,.$

Then one can similarly argue using Subsection D.4 that for each
$i\in\mathcal{I}$, either

$\displaystyle{\bm{y}}_{T}[i]=\sum_{\begin{subarray}{c}t=1,2,\dots,T-1\\\
~{}\text{s.t.}~{}i_{t}=i\end{subarray}}\lambda^{(T)}_{t}\quad\text{or}\quad{\bm{z}}_{T}[i]=\sum_{\begin{subarray}{c}t=1,2,\dots,T-1\\\
~{}\text{s.t.}~{}i_{t}=i\end{subarray}}\lambda^{(T)}_{t}\,.$

This, together with the fact that ${\bm{y}}_{T}^{\sf exact}={\bm{0}}$ and
${\bm{z}}_{T}^{\sf exact}={\bm{0}}$, implies that

$\displaystyle\left\lVert{\bm{y}}_{T}-{\bm{y}}_{T}^{\sf
exact}\right\rVert^{2}+\left\lVert{\bm{z}}_{T}-{\bm{z}}_{T}^{\sf
exact}\right\rVert^{2}\geq\sum_{i\in\mathcal{I}}\left(\sum_{\begin{subarray}{c}t=1,2,\dots,T-1\\\
~{}\text{s.t.}~{}i_{t}=i\end{subarray}}\lambda^{(T)}_{t}\right)^{2}$
$\displaystyle\quad\overset{(a)}{\geq}\frac{1}{|\mathcal{I}|}\cdot\left(\sum_{i\in\mathcal{I}}\sum_{\begin{subarray}{c}t=1,2,\dots,T-1\\\
~{}\text{s.t.}~{}i_{t}=i\end{subarray}}\lambda^{(T)}_{t}\right)^{2}=\frac{1}{4}\cdot\frac{1}{|\mathcal{I}|}\cdot\left(\sum_{t=1}^{T-1}\lambda^{(T)}_{t}\right)^{2}\overset{\eqref{large_coeff_nonsmooth}}{\gtrsim}\frac{1}{T\varepsilon^{2}}.$

On the other hand, since $u_{T}^{\sf exact}=0$, the deviation in the dummy
coordinate can be lower bounded as follows:

$\displaystyle|u_{T}-u_{T}^{\sf
exact}|^{2}=\left(\frac{\delta}{\sqrt{2}}\cdot\sum_{t=0}^{T-1}\lambda^{(T)}_{t}\right)^{2}\gtrsim\frac{\delta^{2}}{\varepsilon^{2}}\,.$

Therefore, combining all together, we get

$\displaystyle\left\lVert{\bm{x}}_{T}-{\bm{x}}_{T}^{\sf
exact}\right\rVert^{2}$
$\displaystyle\geq\left\lVert{\bm{y}}_{T}-{\bm{y}}_{T}^{\sf
exact}\right\rVert^{2}+\left\lVert{\bm{z}}_{T}-{\bm{z}}_{T}^{\sf
exact}\right\rVert^{2}+\left\lVert u_{T}-u_{T}^{\sf exact}\right\rVert^{2}$
$\displaystyle\gtrsim\frac{1}{T\varepsilon^{2}}+\frac{\delta^{2}}{\varepsilon^{2}}\,,$

as desired. ∎

####  D.5 Inexact initialization model

[] (Lower Bound) Let $\varepsilon>0$ be a small constant, and $T$ be a given
number of iterations. There exists a $O(1)$-Lipschitz (nonsmooth) convex
function $f:\mathbb{R}^{O(T)}\to\mathbb{R}$ such that for any FOI algorithm
$\mathcal{A}$ the $(\varepsilon,\delta)$-deviation lower bounded by
$\Omega(\frac{1}{T\varepsilon^{2}}+\delta^{2})$ w.r.t. the reference point
${\bm{x}}_{0}^{\sf ref}={\bm{0}}$.

###### Proof.

For the initialization error model, we use a simpler construction. For
${\bm{x}}=({\bm{x}}^{\sf
err},{\bm{y}},w,u)\in(\mathbb{R}^{T})^{2}\times\mathbb{R}^{2}\to\mathbb{R}$,
consider the cost defined as

$\displaystyle f({\bm{x}}^{\sf err},{\bm{y}},w)=\max\\{0,~{}{\bm{x}}^{\sf
err}[1]+{\bm{y}}[1],~{}\dots,~{}{\bm{x}}^{\sf
err}[T]+{\bm{y}}[T]\\}+2\varepsilon\cdot\max\\{w+1,0\\}\,.$

Here $u$ is a dummy coordinate that does not appear in the cost. For both max
terms above, we consider the subgradient that outputs the gradient of the
first argument that achieves the maximum. Then, $f$ is clearly
$O(1)$-Lipschitz.

We set the reference and inexact intializations as follows:

$\displaystyle{\bm{x}}_{0}^{\sf
ref}=(0,0,\dots,0)\quad\text{and}\quad{\bm{x}}_{0}=(\underbrace{\delta/\sqrt{2T},\cdots,\delta/\sqrt{2T}}_{\text{first
$T$}},\underbrace{0,\cdots,0}_{\text{second $T$}},0,\delta/\sqrt{2})\,.$
(D.17)

Following the previous notations, we write the iterates as
${\bm{x}}_{t}=({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},w_{t})\in(\mathbb{R}^{T})^{2}\times\mathbb{R}$.
Moreover, we will write the iterates corresponding to the reference
initialization as ${\bm{x}}_{t}^{\sf ref}=(({\bm{x}}^{\sf err}_{t})^{\sf
ref},{\bm{y}}_{t}^{\sf ref},w_{t}^{\sf ref})$.

From the fact that the algorithm has to achieve $\varepsilon$-suboptimality,
it follows that $w_{T}\leq-1/2$. On the other hand, from the inexact
initialization (D.17), it holds that for all $t=0,\dots,T-1$,

$\displaystyle\exists
i_{t}\in\\{1,\dots,T\\}~{}s.t.~{}\begin{cases}\frac{\partial}{\partial{\bm{y}}[i_{t}]}f({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},w_{t})=1&\text{and}\\\
\frac{\partial}{\partial{\bm{y}}[i]}f({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},w_{t})=0&\text{for }i\neq i_{t}.\end{cases}$

Hence, it holds that
$\frac{1}{2\varepsilon}w_{T}=\sum_{i=1}^{T}{\bm{y}}_{T}[i]$. Since we know
${\bm{y}}_{T}^{\sf ref}=(0,0,\dots,0)$ and $u^{\sf ref}=0$, the following
deviation lower bound holds:

$\displaystyle\left\lVert{\bm{x}}_{T}-{\bm{x}}_{T}^{\sf ref}\right\rVert^{2}$
$\displaystyle\geq\left\lVert{\bm{y}}_{T}-{\bm{y}}_{T}^{\sf
ref}\right\rVert^{2}+\left\lVert u_{T}-u_{T}^{\sf
ref}\right\rVert^{2}=\sum_{i=1}^{T}({\bm{y}}_{T}[i])^{2}+u_{T}^{2}$
$\displaystyle\geq\frac{1}{T}\cdot\left(\sum_{i=1}^{T}{\bm{y}}_{T}[i]\right)^{2}+\frac{\delta^{2}}{2}=\frac{1}{T}\cdot\frac{1}{4\varepsilon^{2}}\cdot\left(w_{T}\right)^{2}+\frac{\delta^{2}}{2}\gtrsim\frac{1}{T\varepsilon^{2}}+\delta^{2}$

This completes the proof. ∎

####  D.6 Stochastic inexact gradient model (strongly-convex costs)

[] (Lower Bound) Let $\varepsilon>0$ and $T$ be a given number of iterations.
There exists a $O(1)$-Lipschitz (nonsmooth) and $\mu$-strongly convex function
$f:\mathbb{R}^{O(T)}\to\mathbb{R}$ with a stochastic inexact gradient model
such that any FOI algorithm $\mathcal{A}$ that starts at
${\bm{x}}_{0}={\bm{0}}$ has its $(\varepsilon,\delta)$-deviation lower bounded
by $\Omega(\frac{1}{T\mu^{2}}\wedge\frac{\varepsilon}{\mu})$.

###### Proof.

For ${\bm{x}}=({\bm{x}}^{\sf
err},{\bm{y}},{\bm{z}},w)\in(\mathbb{R}^{T})^{3}\times\mathbb{R}^{2}$,
consider the cost defined as

$\displaystyle f({\bm{x}}^{\sf
err},{\bm{y}},{\bm{z}},w)=\underbrace{\mathcal{G}({\bm{x}}^{\sf
err}+\delta{\bm{e}}_{1},{\bm{y}},{\bm{z}})+\frac{\mu}{2}\left\lVert({\bm{x}}^{\sf
err},{\bm{y}},{\bm{z}})\right\rVert^{2}}_{=:\mathcal{G}^{\mu}({\bm{x}}^{\sf
err},{\bm{y}},{\bm{z}})}+\underbrace{w+\frac{\mu}{2}w^{2}}_{=:\ell^{\mu}(w)}\,.$
(D.18)

We consider the same inexact gradient oracle defined in (D.9). Then similarly
to Subsection D.3, it holds that for each $t=1,2,\dots,T-1$, there exists
$i_{t}\in\\{1,\dots,t\\}$ such that the following holds:

$\displaystyle\begin{cases}\frac{\partial}{\partial{\bm{y}}[i_{t}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t},v_{t})=1,~{}~{}\frac{\partial}{\partial{\bm{z}}[i_{t}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t},v_{t})=0&\text{with probability }$1/2$,\\\
\frac{\partial}{\partial{\bm{y}}[i_{t}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t},v_{t})=0,~{}~{}\frac{\partial}{\partial{\bm{z}}[i_{t}]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t},v_{t})=1&\text{with probability
}$1/2$.\end{cases}$ (D.19)

Moreover, for $i\neq i_{t}$,
$\frac{\partial}{\partial{\bm{y}}[i]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t},v_{t})=\frac{\partial}{\partial{\bm{z}}[i]}\mathcal{G}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t},v_{t})=0$. Throughout the rest of the
proof, we use the following notations:

$\displaystyle\begin{cases}\bm{g}^{\bm{y}}_{t}=(\bm{g}^{\bm{y}}_{t}[1],\bm{g}^{\bm{y}}_{t}[2],\dots,\bm{g}^{\bm{y}}_{t}[T]):=\nabla_{\bm{y}}\mathcal{G}^{\mu}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t},v_{t}),\\\
\bm{g}^{\bm{z}}_{t}=(\bm{g}^{\bm{z}}_{t}[1],\bm{g}^{\bm{z}}_{t}[2],\dots,\bm{g}^{\bm{z}}_{t}[T]):=\nabla_{\bm{z}}\mathcal{G}^{\mu}({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t},v_{t})\\\
g^{w}_{t}:=\frac{\partial\ell^{\mu}}{\partial w}(w_{t})\end{cases}$

We prove the following crucial result for the proof.

For each $t=0,1,2,\dots,T$, the output of a FOI algorithm satisfies the
following:

$\displaystyle\begin{split}w_{t}=\sum_{i=1}^{T}{\bm{y}}_{t}[i]+{\bm{z}}_{t}[i]\quad\text{and}\quad
g^{w}_{t}=g^{v}_{t}+\sum_{i=1}^{T}\bm{g}^{\bm{y}}_{t}[i]+\bm{g}^{\bm{z}}_{t}[i]\,.\end{split}$
(D.20)

###### Proof.

We prove by induction on $t$. We first prove the statement for $t=0$. Recall
the definition $\mathcal{G}^{\mu}({\bm{x}}^{\sf
err},{\bm{y}},{\bm{z}}):=\mathcal{G}({\bm{x}}^{\sf
err}+\delta{\bm{e}}_{1},{\bm{y}},{\bm{z}})+\frac{\mu}{2}\left\lVert({\bm{x}}^{\sf
err},{\bm{y}},{\bm{z}})\right\rVert^{2}$. Since the first coordinate of
${\bm{x}}^{\sf err}$ is $\delta$, it follows that $g^{w}_{0}=1$,
$\bm{g}^{\bm{y}}_{0}={\bm{e}}_{1}$ and $\bm{g}^{\bm{z}}_{0}=\bm{0}$. Hence,
the statement holds for $t=0$.

Assume that the conclusion holds for some $t$. We will first show that
$w_{t+1}=\sum_{i=1}^{T}{\bm{y}}_{t+1}[i]+{\bm{z}}_{t+1}[i]$. Using the
definition of FOI together with the inductive hypothesis, we have

$\displaystyle w_{t+1}$
$\displaystyle=-\sum_{j=0}^{t}\lambda^{(t+1)}_{j}g^{w}_{j}=-\sum_{j=0}^{t}\lambda^{(t+1)}_{j}\left(\sum_{i=1}^{T}\bm{g}^{\bm{y}}_{j}[i]+\bm{g}^{\bm{z}}_{j}[i]\right)$
$\displaystyle=-\sum_{i=1}^{T}\sum_{j=0}^{t}\lambda^{(t+1)}_{j}\left(\bm{g}^{\bm{y}}_{j}[i]+\bm{g}^{\bm{z}}_{j}[i]\right)=\sum_{i=1}^{T}{\bm{y}}_{t+1}[i]+{\bm{z}}_{t+1}[i].$

Next, we show that
$g^{w}_{t+1}=\sum_{i=1}^{T}\bm{g}^{\bm{y}}_{t+1}[i]+\bm{g}^{\bm{z}}_{t+1}[i]$.
Using the conclusion we just proved, we obtain

$\displaystyle g^{w}_{t+1}$ $\displaystyle=1+\mu
w_{t+1}=1+\mu\left(\sum_{i=1}^{T}{\bm{y}}_{t+1}[i]+{\bm{z}}_{t+1}[i]\right)=\sum_{i=1}^{T}\left(\bm{g}^{\bm{y}}_{t+1}[i]+\bm{g}^{\bm{z}}_{t+1}[i]\right)\,,$

where in the last equality, we used the fact that
$\nabla\mathcal{G}^{\mu}({\bm{x}}^{\sf err}_{t},{\bm{y}}_{t},{\bm{z}}_{t})$ is
zero except for a single coordinate that is equal to $1$. ∎

Note that in order for $f({\bm{x}}_{T})$ to achieve $\varepsilon$-accuracy, it
must be that

$\displaystyle
w_{T}\in\left[-\frac{1}{2\mu}-O(\sqrt{\varepsilon}),-\frac{1}{2\mu}+O(\sqrt{\varepsilon})\right]\,.$

Note that by symmetry, ${\bm{y}}_{T}[i]+{\bm{z}}_{T}[i]$ is a deterministic
quantity. Subsection D.3 ensures that

$\displaystyle{\bm{y}}_{T}[i]=\begin{cases}{\bm{y}}_{T}[i]+{\bm{z}}_{T}[i]&\text{with
probability}~{}$1/2$,\\\ 0&\text{with
probability}~{}$1/2$,\end{cases}\quad\text{for all}~{}i=2,\dots,T.$

Hence, the deviation
$\operatorname*{\mathbb{E}}\left\lVert{\bm{x}}_{T}-\operatorname*{\mathbb{E}}[{\bm{x}}_{T}]\right\rVert^{2}$
is again lower bounded by
$\frac{1}{4}\sum_{i=2}^{T}({\bm{y}}_{T}[i]+{\bm{z}}_{T}[i])^{2}$. Since either
${\bm{y}}_{T}[i]$ or ${\bm{z}}_{T}[i]$ has to be zero, it follows that
$({\bm{y}}_{T}[i]+{\bm{z}}_{T}[i])^{2}=({\bm{y}}_{T}[i])^{2}+({\bm{z}}_{T}[i])^{2}$.
From the $\varepsilon$-suboptimality of ${\bm{x}}_{T}$, it also holds that
$\sum_{i=2}^{T}({\bm{y}}_{T}[i])^{2}+({\bm{z}}_{T}[i])^{2}\leq\frac{2}{\mu}\varepsilon$.
Hence either
$\sum_{i=2}^{T}({\bm{y}}_{T}[i])^{2}+({\bm{z}}_{T}[i])^{2}=\Omega(\frac{\varepsilon}{\mu})$
(in which case the deviation is lower bounded by
$\Omega(\frac{\varepsilon}{\mu})$), or we use Subsection D.6 to conclude that

$\displaystyle\sum_{i=2}^{T}({\bm{y}}_{T}[i]+{\bm{z}}_{T}[i])^{2}$
$\displaystyle\overset{(a)}{\geq}\frac{1}{T-1}\cdot\left(\sum_{i=2}^{T}{\bm{y}}_{T}[i]+{\bm{z}}_{T}[i]\right)^{2}$
$\displaystyle=\frac{1}{T-1}\cdot\left(w_{T}-({\bm{y}}_{T}[1]+{\bm{z}}_{T}[1])\right)^{2}\gtrsim\frac{1}{T\mu^{2}}\,,$

where $(a)$ follows form the Cauchy-Schwartz inequality, and the last
inequality is due to the fact that
$({\bm{y}}_{T}[1]+{\bm{z}}_{T}[1])^{2}\lesssim\varepsilon$. This completes the
proof. ∎

####  D.7 Non-stochastic inexact gradient model (strongly-convex costs)

[] (Lower Bound) Let $\varepsilon>0$ and $T$ be a given number of iterations.
There exists a $O(1)$-Lipschitz (nonsmooth) and $\mu$-strongly convex function
$f:\mathbb{R}^{O(T)}\to\mathbb{R}$ with a non-stochastic inexact gradient
model such that any FOI algorithm $\mathcal{A}$ that starts at
${\bm{x}}_{0}={\bm{0}}$ has its $(\varepsilon,\delta)$-deviation lower bounded
by
$\Omega((\frac{1}{T\mu^{2}}+\frac{\delta^{2}}{\mu^{2}})\wedge\frac{\varepsilon}{\mu})$.

###### Proof.

We consider the same construction as the one considered in the proof of
Subsection D.6. The only difference is that now we add an extra dummy
coordinate $u$ (the $(3T+2)$-th coordinate), and add $\frac{\mu}{2}u^{2}$ to
the overall cost (so that the overall cost is still $\mu$-strongly convex).
Following the previous convention, we will write the iterate as
${\bm{x}}_{t}=({\bm{x}}^{\sf
err}_{t},{\bm{y}}_{t},{\bm{z}}_{t},w_{t},u_{t})\in(\mathbb{R}^{T})^{3}\times\mathbb{R}^{3}$.
We define the inexact gradient as:

$\displaystyle g({\bm{x}}_{t})=\nabla
f({\bm{x}}_{t})+\frac{\delta}{\sqrt{2}}{\bm{e}}_{1+t}+\frac{\delta}{\sqrt{2}}{\bm{e}}_{3T+3}\,,$

where ${\bm{e}}_{j}$ is the $j$-th coordinate vector. Then following the same
argument as Subsection D.6, it holds that

$\displaystyle
w_{T}=\sum_{i=1}^{T}{\bm{y}}_{T}[i]+{\bm{z}}_{T}[i]\quad\text{and}\quad\delta/\sqrt{2}\cdot
w_{T}=u_{T}\,.$ (D.21)

From the fact that ${\bm{x}}_{T}$ achieves $\varepsilon$-suboptimality, it
must be that

$\displaystyle
w_{T}\in\left[-\frac{1}{2\mu}-O(\sqrt{\varepsilon}),-\frac{1}{2\mu}+O(\sqrt{\varepsilon})\right]\,.$

Now the rest of the proof follows similarly to that of Subsection D.6. Using
the facts ${\bm{y}}_{T}^{\sf exact}=0$, ${\bm{z}}_{T}^{\sf exact}=0$ and
$u_{T}^{\sf exact}=0$, we have the following lower bound on the deviation:

$\displaystyle\left\lVert{\bm{x}}_{T}-{\bm{x}}_{T}^{\sf
exact}\right\rVert^{2}$
$\displaystyle\geq\left\lVert{\bm{y}}_{T}-{\bm{y}}_{T}^{\sf
exact}\right\rVert^{2}+\left\lVert{\bm{z}}_{T}-{\bm{z}}_{T}^{\sf
exact}\right\rVert^{2}+\left\lVert u_{T}-u_{T}^{\sf exact}\right\rVert^{2}$
$\displaystyle\geq
u_{T}^{2}+\sum_{i=2}^{T}[({\bm{y}}_{T}[i])^{2}+({\bm{z}}_{T}[i])^{2}]$

From the $\varepsilon$-suboptimality of ${\bm{x}}_{T}$, it must be that
$u_{T}^{2}+\sum_{i=2}^{T}({\bm{y}}_{T}[i])^{2}+({\bm{z}}_{T}[i])^{2}\leq\frac{2}{\mu}\varepsilon$.
Thus, either
$u_{T}^{2}+\sum_{i=2}^{T}({\bm{y}}_{T}[i])^{2}+({\bm{z}}_{T}[i])^{2}=\Omega(\frac{\varepsilon}{\mu})$
(in which case the deviation is lower bounded by
$\Omega(\frac{\varepsilon}{\mu})$), or we use the fact that
$(\sum_{i=2}^{T}{\bm{y}}_{T}[i]+{\bm{z}}_{T}[i])^{2}\gtrsim\frac{1}{\mu^{2}}$
to conclude that

$\displaystyle
u_{T}^{2}+\sum_{i=2}^{T}[({\bm{y}}_{T}[i])^{2}+({\bm{z}}_{T}[i])^{2}]\geq\frac{\delta^{2}}{2}w_{T}^{2}+\frac{1}{2(T-1)}\left(\sum_{i=1}^{T}{\bm{y}}_{T}[i]+{\bm{z}}_{T}[i]\right)^{2}$