# BiSLS/SPS: Auto-tune Step Sizes for
Stable Bi-level Optimization

Chen Fan
Department of Computer Science
University of British Columbia
&Gaspard Choné-Ducasse
Department of Mathematics and Applications
Ecole Normale Supérieure
Mark Schmidt
Department of Computer Science
University of British Columbia
Canada CIFAR AI Chair (Amii)
&Christos Thrampoulidis
Department of Electrical and Computer Engineering
University of British Columbia

###### Abstract

The popularity of bi-level optimization (BO) in deep learning has spurred a
growing interest in studying gradient-based BO algorithms. However, existing
algorithms involve two coupled learning rates that can be affected by
approximation errors when computing hypergradients, making careful fine-tuning
necessary to ensure fast convergence. To alleviate this issue, we investigate
the use of recently proposed adaptive step-size methods, namely stochastic
line search (SLS) and stochastic Polyak step size (SPS), for computing both
the upper and lower-level learning rates. First, we revisit the use of SLS and
SPS in single-level optimization without the additional interpolation
condition that is typically assumed in prior works. For such settings, we
investigate new variants of SLS and SPS that improve upon existing suggestions
in the literature and are simpler to implement. Importantly, these two
variants can be seen as special instances of general family of methods with an
envelope-type step-size. This unified envelope strategy allows for the
extension of the algorithms and their convergence guarantees to BO settings.
Finally, our extensive experiments demonstrate that the new algorithms, which
are available in both SGD and Adam versions, can find large learning rates
with minimal tuning and converge faster than corresponding vanilla SGD or Adam
BO algorithms that require fine-tuning.

##  1 Introduction

Bi-level optimization has found its applications in various fields of machine
learning, such as hyperparameter optimization [14, 17, 30, 40], adversarial
training [51], data distillation [2, 53], neural architecture search [28, 39],
neural-network pruning [52], and meta-learning [13, 37, 11]. Specifically, it
is used widely for problems that exhibit a hierarchical structure of the
following form:

$\displaystyle\min_{x\in
X}F(x)=\mathbb{E}_{\phi}[f(x,y^{*}(x);\phi)]\qquad\text{s.t.}\qquad
y^{*}(x)=\operatorname*{argmin}_{y\in Y}\mathbb{E}_{\psi}[g(x,y;\psi)].$ (1)

Here, the solution to the lower-level objective $g$ becomes the input to the
upper-level objective $f$, and in (1) the upper-level variable $x$ is fixed
when optimizing the lower-level variable $y$. To solve such bi-level problems
using gradient-based methods requires computing the hypergradient of $F$,
which based on the chain rule is given as [15]:

$\displaystyle\nabla
F(x)=\nabla_{x}f(x,y^{*}(x))+\nabla_{xy}^{2}g(x,y^{*}(x))[\nabla_{yy}^{2}g(x,y^{*}(x))]^{-1}\nabla_{y}f(x,y^{*}(x)).$
(2)

In practice, the closed-form solution $y^{*}(x)$ is difficult to obtain, and
one strategy is to run a few steps of (stochastic) gradient descent on $g$
w.r.t. $y$ to get an approximation $\bar{y}$, and use $\bar{y}$ in places of
$y^{*}(x)$. We denote the stochastic hypergradient based on $\bar{y}$ as
$h_{f}(x,\bar{y})$ and the stochastic gradient of $g$ w.r.t. $y$ as $h_{g}$.
This leads to a general gradient-based framework for solving bi-level
optimization [15, 19, 4]. At each iteration $k$, run T (can be one or more)
steps of SGD on $y$, i.e. $y^{k,t+1}=y^{k,t}-\beta h_{g}^{k,t}$, then run one
step on $x$ using the approximated hypergradient:

$\displaystyle x^{k+1}=x^{k}-\alpha
h_{f}(x^{k},y^{k+1}),\quad\text{where}\quad y^{k+1}=y^{k,T}.$ (3)

Based on this framework, a series of stochastic algorithms have been developed
to achieve the optimal or near-optimal rate of their deterministic
counterparts [7, 8]. These algorithms can be broadly divided into single-loop
($T=1$) or double-loop ($T>1$) categories [23].

Figure 1: Results based on hyper-representation learning task (see Sec 4 for
details). Validation loss against upper-level iterations for different values
of $\beta$ (left, $\alpha=0.005$) and $\alpha$ (right, $\beta=0.01$). Unless
carefully tuned, vanilla SGD-based methods for BO are very unstable.

Unlike minimizing the single-level finite-sum (convex) problem

$\displaystyle F(x):=\min_{x\in\mathcal{C}}\frac{1}{N}\sum_{i=1}^{N}f_{i}(x),$
(4)

where only one learning rate is involved when using SGD, bi-level optimization
involves tuning both the lower and upper-level learning rates ($\beta$ and
$\alpha$ respectively). This poses a significant challenge due to the
potential correlation between these learning rates [19]. Thus, as observed in
Figure 1, algorithm divergence can occur when either $\alpha$ or $\beta$ is
large. While there is considerable literature on achieving faster rates in bi-
level optimization [24, 5, 7, 8], only a few studies have focused on
stabilizing its training and automating the tuning of $\alpha$ and $\beta$.
This work addresses the question: Is it possible to utilize large $\alpha$ and
$\beta$ without manual tuning?

Algorithm 1 BiSLS-Adam/SGD

1:$x^{0}$, $y^{0}$, $K$, $T$, $\delta$, $\alpha_{b,0}$, $\beta_{b,0}$, $w$,
$\eta$

2:$x$

3:for $k=0,1,\dots,K-1$ do

4: $y^{k,0}=y^{k}$

5: for $t=0,1,\dots,T-1$ do

6: $\beta_{b,k}^{t}\leftarrow$ reset($\beta$, $\beta_{b,0}$, $\eta$, opt)
$\triangleright$ see Algorithm 2

7: $\beta\leftarrow$ linesearch based on (8) starting from $\beta_{b,k}^{t}$

8: $y^{k,t+1}=y^{k,t}\,-\,\beta\,h_{g}^{k,t}$,

9: end for

10: $y^{k+1}=y^{k,T-1}$;$\hat{x}^{k}=x^{k}$; $\hat{y}^{k+1}=y^{k+1}$

11: $\alpha\leftarrow$ reset($\alpha$, $\alpha_{b,0}$, $\eta$, opt)

12: while (14) based on ($\hat{x}^{k}$, $\hat{y}^{k+1}$, $\alpha$, $\delta$)
does not hold. do

13: $\alpha=\alpha*w$

14: $\hat{x}^{k}=x^{k}\,-\,\alpha\,h_{f}(x^{k},y^{k+1})$ or

15: $\hat{x}^{k}=x^{k}-\alpha A_{k}^{-1}h_{f}(x^{k},y^{k+1})$

16: $\hat{y}^{k+1}=y^{k+1}\,-\,\beta\,h_{g}(\hat{x}^{k},y^{k+1})$

17: end while

18: $x^{k+1}=x^{k}\,-\,\alpha\,h_{f}^{k}$

19:end for

In doing so, we explore the use of stochastic adaptive-step size methods,
namely stochastic Polyak step size (SPS) and stochastic line search (SLS),
which utilize gradient information to adjust the learning rate at each
iteration [44, 29]. These methods have been demonstrated to perform well in
interpolation settings with strong convergence guarantees [44, 29]. However,
applying them to bi-level optimization (BO) introduces significant challenges,
as follows. ① BO requires tuning two correlated learning rates (for lower and
upper-level). ② The bias in the stochastic approximation of the hypergradient
complicates the practical performance and convergence analysis of SLS and SPS.
③ Other algorithmic challenges arise for both algorithms: For SLS, verifying
the stochastic Armijo condition at the upper-level involves evaluating the
objective at a new $(x,y^{*}(x))$ pair, while $y^{*}(x)$ is only approximately
known; For SPS, most existing variants guarantee good performance only in
interpolating settings, which are typically not satisfied for the upper-level
objective in BO [22]. Before presenting our solutions to the challenges above
in Sec 2, we first review the most closely related literature.

###  1.1 Related Work

##### Gradient-Based Bi-level Optimization

Penalty or gradient-based approaches have been used for solving bi-level
optimization problems [10, 45, 21]. Here we focus our discussions on
stochastic gradient-based methods as they are closely related to this work.
For double-loop algorithms, an early work (BSA) by Ghadimi and Wang [15] has
derived the sample complexity of $\phi$ in achieving an $\epsilon$-stationary
point to be $\mathcal{O}(\epsilon^{-2})$, but require the number of lower-
level steps to satisfy $T\sim\mathcal{O}(\epsilon^{-1})$. Using a warm start
strategy (stocBiO), Ji et al. [22] removed this requirement on $T$. However,
to achieve the same sample complexity, the batch size of stocBiO grows as
$\mathcal{O}(\epsilon^{-1})$. Chen et al. [4] removed both requirements on $T$
and batch size by using the smoothness properties of $y^{*}(x)$ and setting
the step sizes $\alpha$ and $\beta$ at the same scale. For single-loop
algorithms, a pioneering work by Hong et al. [19] gave a sample complexity of
$\mathcal{O}(\epsilon^{-2.5})$, provided $\alpha$ and $\beta$ are on two
different scales (TTSA). By making corrections to the $y$ variable update
(STABLE), Chen et al. [5] improved the rate to $\mathcal{O}(\epsilon^{-2})$.
However, extra matrix projections required by STABLE can incur high
computation cost [5, 4]. By incorporating momentum into the updates of $x$ and
$y$ (SUSTAIN), Khanduri et al. [24] further improved the rate to
$\mathcal{O}(\epsilon^{-1.5})$ [6]. Besides these single or double-loop
algorithms, a series of works have drawn ideas from variance reduction to
achieve faster convergence rates for BO. For example, Yang et al. [49]
designed the VRBO algorithm based on SPIDER [12]. Dagréou et al. [7, 8]
designed the SABA and SRBA algorithms based on SAGA and SARAH respectively,
and demonstrate that they can achieve the optimal rate of
$\mathcal{O}(\epsilon^{-1})$ [9, 35].

Huang et al. [20] proposes to use Adam-type step sizes in BO. However, it
introduces three sequences of learning rates $(\alpha_{k},\beta_{k},\eta_{k})$
that require tuning, which limits its practical usage. To our knowledge, none
of these works have explicitly addressed the fundamental problem of how to
select $\alpha$ and $\beta$ in bi-level optimization. In this work, we focus
on the alternating SGD framework (T can be $1$ or larger), and design
efficient algorithms that find large $\alpha$ and $\beta$ without tuning,
while ensuring the stability of training.

Figure 2: Experiments on quadratic functions adapted from [29]. The objective
is the sum of two-dimensional functions
$f_{i}=\frac{1}{2}(x-x_{i}^{*})^{T}H_{i}(x-x_{i}^{*})$, where $H_{i}$ is
positive definite and $i=1,2$ (see Appendix B for more details). From left to
right, we show: the objective value, distance to optimum, step size, and
iterate trajectories.

##### Adaptive Step Size

Adaptive step-size such as Adam has found great success in modern machine
learning, and different variants have been proposed [25, 38, 47, 31, 32].
Here, we limit our discussions on two adaptive step sizes that are most
relevant to this work. The Armijo line search is a classic way for finding
step sizes for gradient descent [48]. Vaswani et al. [44] extends it to the
stochastic setting (SLS) and demonstrates that the algorithm works well with
minimal tuning required under interpolation, where model fits the data
perfectly. Hence, the method is adaptive to local smoothness of the objective,
which is typically difficult to predict a priori. However, the theoretical
guarantee of SLS in the non-interpolating regime is lacking. In fact, the
results in Figure 3 suggest that SLS can perform poorly for convex losses when
interpolation is not satisfied. Besides SLS, another adaptive method derived
from Polyak step size is proposed by Loizou et al. [29] with the name
stochastic Polyak step size (SPS). Loizou et al. [29] further places an upper
bound on the step size resulting in the $\text{SPS}_{\max}$ variant. Similar
to SLS, the algorithm performs well when the model is over-parametrized.
Without interpolation, the algorithm converges to a neighborhood of the
solution whose size depends on this upper bound. In a later work, Orvieto et
al. [36] make the SPS converge to the exact solution by ensuring the step size
and its upper bound are both non-increasing (DecSPS ). However, enforcing
monotonicity may result in the step size being smaller than decaying-step SGD
and losing the adaptive features of SPS (see Figure 2, 3). In this work, we
propose new versions of SLS and SPS that do not require monotonicity and
extend them into the alternating SGD bi-level optimization framework (3).

##  2 Summary of Contributions

We discuss our main contributions in this section, which is organized as
follows. First, we discuss our variants of SPS and SLS, and unify them under
the umbrella of “envelope-type step-size”. Then, we extend the envelope-type
step size to the bi-level setting. Finally, we discuss our bi-level line-
search algorithms based on Adam and SGD.

##### Converging SPSB and SLSB by Envelope Approach

We first propose simple variants of SLS and SPS that converge in the non-
interpolating setting while not requiring the step size to be monotonic. To
this end, we introduce a new stochastic Polyak step size (SPSB). For
comparison, we also recall the step-sizees of $\text{SPS}_{\max}$ and DecSPS .
For all methods, the iterate updates are given as
$x_{k+1}=x_{k}-\gamma_{k}\nabla f_{i_{k}}(x^{k})$ where $i_{k}$ is sampled
uniformly from $[n]=\\{1,\ldots,n\\}$ at each iteration $k$. The step-sizes
$\gamma_{k}$ are then defined as follows:

$\displaystyle\text{$\text{SPS}_{\max}$~{}\cite[cite]{[\@@bibref{Number}{loizou2021stochastic}{}{}]}:
}\quad\gamma_{k}=\min\\{\frac{f_{i_{k}}(x^{k})-f_{i_{k}}^{*}}{c\lVert\nabla
f_{i_{k}}(x^{k})\rVert^{2}},\gamma_{b,0}\\}$ (5)
$\displaystyle\text{DecSPS~{}\cite[cite]{[\@@bibref{Number}{orvieto2022dynamics}{}{}]}:}\quad\gamma_{0}=\bar{\gamma}\quad\gamma_{k}=\frac{1}{c_{k}}\min\\{\frac{f_{i_{k}}(x^{k})-l_{i_{k}}^{*}}{\lVert\nabla
f_{i_{k}}(x^{k})\rVert^{2}},c_{k-1}\gamma_{k-1}\\}\quad\forall k\geq 1$ (6)
$\displaystyle\textbf{SPSB~{}(ours):
}\quad\gamma_{k}=\min\\{\frac{f_{i_{k}}(x^{k})-l_{i_{k}}^{*}}{c_{k}\lVert\nabla
f_{i_{k}}(x^{k})\rVert^{2}},\gamma_{b,k}\\},$ (7)

where $f_{i}^{*}=\inf_{x}f_{i}(x)$,
$\bar{\gamma}=\frac{1}{c_{0}}\min\\{\frac{f_{i_{0}}(x^{0})-l_{i_{0}}^{*}}{\lVert\nabla
f_{i_{0}}(x^{0})\rVert^{2}},c_{0}\gamma_{b,0}\\}$ , $c_{k}$ is non-decreasing,
$\gamma_{b,k}$ is non-increasing, and $l_{i}^{*}\leq f_{i}^{*}$ is any lower
bound.

Figure 3: Binary linear classification on w8a dataset using logistic loss [3].
We choose $\gamma_{b,0}=1000$ for all algorithms; $c=1$ and $\bar{c}=1$ for
$\text{SPS}_{\max}$ and SLS respectively; $c_{k}=\sqrt{k+1}$ for DecSPS ;
$c_{k}=1$ and $\gamma_{b,k}=\frac{\gamma_{b,0}}{\sqrt{k+1}}$ for SPSB ;
$\bar{c}=0.1$ and $\gamma_{b,k}=\frac{\gamma_{b,0}}{\sqrt{k+1}}$ for SLSB ;
$\gamma_{b,k}=\frac{\gamma_{b,0}}{\sqrt{k+1}}$ for decaying-step SGD.

Unlike $\text{SPS}_{\max}$ in which $\gamma_{b,0}$ is a constant, our upper
bound $\gamma_{b,k}$ is non-increasing. Also, unlike DecSPS in which both the
step size and the upper bound are non-increasing (this is because
$\gamma_{k}\leq\frac{c_{k-1}}{c_{k}}\gamma_{k-1}$ and
$\min\\{\frac{1}{2cL_{\max}},\frac{c_{0}\gamma_{b,0}}{c_{k}}\\}\leq\gamma_{k}\leq\frac{c_{0}\gamma_{b,0}}{c_{k}}$
[36, Lemma 1]), we simplify the recursive structure and do not require the
step-size to be monotonic. As we empirically observe in Figure 3, the step
size of DecSPS is similar to that of decaying SGD and in fact can be much
smaller. Interestingly, the resulting performance of DecSPS is worse than
$\text{SPS}_{\max}$ despite $\text{SPS}_{\max}$ eventually becoming unstable
once iterates get closer to the neighborhood of a solution and the step-size
naturally behaves erratically. This is not unexpected due to small gradient
norms (note division by gradient-norm in (5)) and dissimilarity between
samples in the non-interpolating scenario. Moreover, note that the adaptivity
of SPS in the early stage seems to be lost in DecSPS due to monotonicity of
the latter. On the other hand, SPSB not only takes advantage of the large SPS
steps that leads to fast convergence, but also stays regularized due to the
non-increasing upper bound $\gamma_{b,k}$ in (19). These observations are
further supported by the experiments on quadratic functions given in Figure 2,
where we observe the fast convergence of SPSB and the instability of
$\text{SPS}_{\max}$ . Motivated by the good practical performance of SPSB , we
take a similar approach for SLS. The SLS proposed and analyzed by Vaswani et
al. [44] starts with $\gamma_{b,0}$ and in each iteration $k$ finds the
largest $\gamma_{k}\leq\gamma_{b,0}$ that satisfies:

$\displaystyle f_{i_{k}}(x_{k}-\gamma_{k}\nabla f_{i_{k}}(x_{k}))\leq
f_{i_{k}}(x_{k})-\bar{c}\cdot\gamma_{k}\lVert\nabla
f_{i_{k}}(x_{k})\rVert^{2},\quad 0<\bar{c}<1.$ (8)

To ensure its convergence without interpolation, we replace $\gamma_{b,0}$
with appropriate non-increasing sequence $\gamma_{b,k}$. We name this variant
of SLS as SLSB . Interestingly, the empirical performance and step size of
SLSB are similar to those of SPSB (see Figure 3). This can be explained by
observing that the step sizes of SPSB and SLSB share similar envelope
structures, as follows (see Lemma 1 in Appendix A):

$\displaystyle\text{SPSB~{}:
}\quad\min\\{\frac{1}{2cL_{\max}},\gamma_{b,k}\\}\leq\gamma_{k}=\min\\{\frac{f_{i_{k}}(x^{k})-l_{i_{k}}^{*}}{c\lVert\nabla
f_{i_{k}}(x^{k})\rVert^{2}},\gamma_{b,k}\\},\quad 0<c,$
$\displaystyle\text{SLSB~{}:
}\quad\min\\{\frac{2(1-\bar{c})}{L_{\max}},\gamma_{b,k}\\}\leq\gamma_{k}\leq\min\\{\frac{f_{i_{k}}(x^{k})-l_{i_{k}}^{*}}{\bar{c}\lVert\nabla
f_{i_{k}}(x^{k})\rVert^{2}},\gamma_{b,k}\\},\quad 0<\bar{c}<1.$

Therefore, we unify their analysis based on the following generic _envelope-
type step size_ :

$\displaystyle\gamma_{k}=\min\\{\max\\{\gamma_{l,k},\tilde{\gamma}_{k}\\},\gamma_{b,k}\\},\quad\gamma_{l,k}=\min\\{w,\gamma_{b,k}\\},$
(9)

where $\omega>0$, $\gamma_{b,k}$ is non-increasing, and $\tilde{\gamma}_{k}$
satisfies
$\gamma_{l,k}:=\min\\{\omega,\gamma_{b,k}\\}\leq\tilde{\gamma}_{k}\leq\gamma_{b,k}$.
We show that this envelope-type step size converges at a rate
$\mathcal{O}(\frac{1}{\sqrt{K}})$ and $\mathcal{O}(\frac{1}{K})$ for convex
and strongly-convex losses respectively.

##### Envelope Step Size for Bi-level Optimization (BiSPS)

We extend the analysis of envelope-type step size to the bi-level setting. The
step sizes for upper and lower-level objectives of our general envelope-type
method are:

$\displaystyle\text{Upper:
}\alpha_{k}=\min\\{\max\\{\alpha_{l,k},\tilde{\alpha}_{k}\\},\alpha_{b,k}\\}\quad\text{hence}\quad\alpha_{l,k}\leq\tilde{\alpha}_{k}\leq\alpha_{b,k}$
(10) $\displaystyle\text{Lower:
}\beta_{k,t}=\min\\{\frac{g(x^{k},y^{k,t};\psi)-g(x^{k},y^{*}_{x^{k},\psi};\psi)}{p\lVert\nabla_{y}g(x^{k},y^{k,t};\psi)\rVert^{2}},\beta_{b,k}\\}\quad\forall
t,$ (11)

where $y^{*}_{x^{k},\psi}$ is the minimizer of the function
$g(x^{k},\cdot;\psi)$, and $\alpha_{l,k}$, $\alpha_{b,k}$, and $\beta_{b,k}$
are three non-increasing sequences. Note that $\beta_{b,k}$ is fixed over the
lower-level iterations for a given $k$, therefore, this is equivalent to
running $T$ steps of $\text{SPS}_{\max}$ to minimize the function $g$ at each
upper iteration $k$. However, the decrease in the upper bound $\beta_{b,k}$
with $k$ is crucial to guarantee the overall convergence of the algorithm (see
Theorem 3). Starting from the general step-size rules in (10), (11), our bi-
level extension of SPS, which we call BiSPS, follow by setting $\alpha_{k}$ in
the form of SPS computed using stochastic hypergradient $h_{f}^{k}$. That is,

$\displaystyle\bar{\alpha}_{k}=\frac{f(x^{k},y^{k+1};\phi)-l_{f(\cdot,y^{k+1};\phi)}^{*}}{p\lVert
h_{f}^{k}\rVert^{2}},\quad\alpha_{l,k}=\frac{\alpha_{l,0}}{\sqrt{k+1}},\quad\alpha_{b,k}=\frac{\alpha_{b,0}}{\sqrt{k+1}},$
(12)

where $\alpha_{l,0}\leq\alpha_{b,0}$ and $l_{f(\cdot,y^{k+1};\phi)}^{*}$ is a
lower bound for $\inf_{x}f(x,y^{k+1};\phi)$. For computing $h_{f}^{k}$, we can
take a similar approach as previous works [15, 19, 4] that use Neumann series
setting

$\displaystyle
h_{f}^{k}=\nabla_{x}f(x^{k},y^{k+1};\phi)-\nabla_{xy}g(x^{k},y^{k+1};\psi_{0})\bigl{[}\frac{N}{L_{g}}\prod_{j=1}^{\bar{N}}(I-\nabla^{2}_{yy}g(x^{k},y^{k+1};\psi_{j}))\bigr{]}\nabla_{y}f(x^{k},y^{k+1};\phi),$
(13)

where $\bar{N}$ is sampled uniformly from $[N]$ and $N$ is the total number of
samples.

Figure 4: Results on data distillation experiments adapted from Lorraine et
al. [30] (see Sec 4 for details). We compare BiSPS and decaying-step SGD for
different values of $\alpha_{b,0}$ where Hessian inverse in (2) is computed
based on the Identity matrix (left) or Neumann series (right). The lower-level
learning rate is fixed at $10^{-4}$.

For BiSPS, we use the same sample for $f(x^{k},y^{k+1};\phi)$ and $\nabla
f_{x}(x^{k},y^{k+1};\phi)$ when evaluating $\bar{\alpha}_{k}$ in (12).
Interestingly, we also empirically observe that using independent samples for
computing $\bar{\alpha}_{k}$ and $h_{f}^{k}$ resulting in similar performance
as using the same sample. The optimal rate of SGD for non-convex bi-level
optimization is $\mathcal{O}(\frac{1}{\sqrt{K}})$ without a growing batch size
[4]. We show that BiSPS can obtain the same rate (see Theorem 3) by taking the
envelope-type step-size of the form (10) and (11). We implement BiSPS
according to (12) and observe that it has better performances over decaying-
step SGD with less variations across different values of $\alpha_{b,0}$ (see
Figure 4 and note that decaying-step SGD is of the form
$\frac{\alpha_{b,0}}{\sqrt{k+1}}$).

##### Stochastic Line-Search Algorithms for Bi-level Optimization

The challenge of extending SLS to bi-level optimization is rooted in the term
$y^{*}(x)$. In fact, we realize that some of the bi-level objectives are of
the form $F(x)=f(y^{*}(x))$. That is, $f$ does not have an explicit dependence
on $x$, e.g. the data hyper-cleaning task [22]. This implies that when SLS
takes a potential step on $x$, the approximation of $y^{*}(x)$ (i.e,
$\bar{y}(x)$) also needs to be updated, otherwise there is no change in
function values. Moreover, the use of approximation $\bar{y}(x)$ and the
stochastic estimation error in hypergradient would not gaurantee a step size
can be always found. To this end, we modify the Armijo line-search rule to be:

$\displaystyle\text{BiSLS-SGD: }\quad
f\bigl{(}x^{k}-\alpha_{k}h_{f}^{k},\hat{y}^{k+1}(x^{k}-\alpha_{k}h_{f}^{k})\bigr{)}\leq
f(x^{k},y^{k+1})-p\alpha_{k}\lVert h_{f}^{k}\rVert^{2}+\delta,$ (14)
$\displaystyle\text{BiSLS-Adam: }\quad
f\bigl{(}x^{k}-\alpha_{k}A_{k}^{-1}h_{f}^{k},\hat{y}^{k+1}(x^{k}-\alpha_{k}A_{k}^{-1}h_{f}^{k})\bigr{)}\leq
f(x^{k},y^{k+1})-p\alpha_{k}\lVert h_{f}^{k}\rVert_{A_{k}^{-1}}^{2}+\delta,$

where $p,\delta>0$ and $A_{k}$ is a positive definite matrix such that
$A_{k}^{2}=G_{k}$.

(a) Validation Loss

(b) Upper-level step sizes

(c) Lower-level step sizes

Figure 5: Results on hyper-representation learning task (see Sec 4 for
details). (a) Validation loss against upper-level iterations for comparing
BiSLS-Adam/SGD to fine-tuned Adam/SGD. (b)(c) Upper (left) and lower-level
(right) learning rates found by BiSLS-Adam. For the tuned Adam, the optimal
lower and upper-level learning rates are $\mathcal{O}(1)$ and
$\mathcal{O}(10^{-4})$, respectively. BiSLS-Adam outperforms tuned Adam/SGD
with a starting point that is $5$ orders of magnitude larger than the optimal
step size.

Similar to the single-level Adam case, the matrix $G_{k}$ in the bi-level
setting is defined as
$G_{k}=(\beta_{2}G_{k-1}+(1-\beta_{2})\operatorname*{diag}(h_{f}^{k}{h_{f}^{k}}^{T}))/(1-\beta_{2}^{k})$
[25, 43]. Moreover, BiSLS-Adam takes the following steps for updating the
variable $x$: $x^{k+1}=x^{k}-\alpha_{k}A_{k}^{-1}m_{k}$ where
$m^{k+1}=\beta_{1}m^{k}-(1-\beta_{1})h_{f}^{k}$. The details are given in
Algorithms 1 and 2. We denote the search starting point for the upper-level as
$\alpha_{b,k}$ at iteration $k$, and denote it as $\smash{\beta_{b,k}^{t}}$ at
step $t$ within iteration $k$ for the lower-level. We remark the following key
benefits of resetting $\alpha_{b,k}$ and $\beta_{b,k}^{t}$ (by using Algorithm
2) to larger values with reference to $\alpha_{k}$ and $\beta_{k}^{t}$
(respectively) at each step: (1) Avoid always searching from $\alpha_{b,0}$ or
$\beta_{b,0}^{0}$, thus, reducing computation cost, and, (2) preserving an
overall non-increasing (not necessarily monotonic) trend for $\alpha_{b,k}$
and $\smash{\beta_{b,k}^{t}}$, thus, improving training stability. We found
different values of $\eta$ all work well (see Appendix B). The key algorithmic
challenge we are facing is that during the backtracking process, for any
candidate $\alpha_{k}$, we need to compute
$\hat{x}^{k}:=x^{k}-\alpha_{k}h_{f}^{k}$ and approximate
$\smash{y^{*}(\hat{x}^{k})}$ with $\smash{\hat{y}^{k+1}}$ (see Algorithm 1).
To limit the cost induced by this nested loop, we limit the number of steps to
obtain $\smash{\hat{y}^{k+1}}$ to be 1. Moreover, $\delta$ in (14) plays the
role of a safeguard that ensures a step size can be found.

Algorithm 2 reset

1:$p$, $q$, $\eta\geq 1$, opt

2:$p$

3:if opt = 1 then

4: $p\leftarrow q$

5:else if opt = 2 then

6: $p\leftarrow p$

7:else if opt = 3 then

8: $p\leftarrow\eta\cdot p$

9:end if

We set it to be small to avoid finding unrealistically large learning rates
while tolerating some error in the hypergradient estimation (see Appendix B
for experiments on the sensitivity of $\delta$). In practice, we empirically
find that simply setting $\delta=0$ works well. In Figure 5(a), we observe
that BiSLS-Adam outperforms fine-tuned Adam or SGD. Surprisingly, its training
is stable even when the search starting point $\alpha_{b,0}$ is $5$ orders of
magnitude larger than a fine-tuned learning rate ($\mathcal{O}(10^{-4})$).
Importantly, BiSLS-Adam finds large upper and lower-level learning rates in
early phase (see Figure 5(b), 5(c)) for different values of $\alpha_{b,0}$ and
$\beta_{b,0}$ that span $3$ orders of magnitudes. Interestingly, the learning
rates naturally decay with training (also see Figure 6(c) and 6(d)). In
essence, BiSLS is a truly adaptive (no knowledge of initialization required)
and robust (different initialization works) method that finds large $\alpha$
and $\beta$ without tuning. In the next section, we give the convergence
results of the envelope-type step size.

##  3 Convergence Results

###  3.1 Envelope-type step size for single-level optimization

We first state the assumptions, which are standard in the literature, that
will be used for analyzing single-level problems. Assumption 1 is on the
Lipschitz continuity of $f$ and $f_{i}$ in Problem 4.

######  Assumption 1.

The individual function $f_{i}$ is convex and $L_{i}$-smooth such that
$\lVert\nabla f_{i}(x)-\nabla f_{i}(x^{\prime})\rVert\leq L_{i}\lVert
x-x^{\prime}\rVert,\forall i,\forall x\in\operatorname{dom}{f}$ and the
overall function $f$ is $L$-smooth. We denote
$L_{\max}\triangleq\textstyle\max_{i}L_{i}$. Furthermore, we assume there
exists $l_{i}^{*}$ such that $l_{i}^{*}\leq
f_{i}^{*}:=\inf_{x}f_{i}(x),\forall i$, and $f$ is lower bounded by $f^{*}$
obtained by some $x^{*}$ such that $f^{*}=f(x^{*})$.

The following bounded gradient assumption is also used in the analysis of
convex problems [41, 33].

######  Assumption 2.

There exists $G>0$ such that $\lVert\nabla f_{i}(x)\rVert^{2}\leq G,\forall
i$.

We first state the theorem for the envelop-type step size defined in (9) for
convex functions.

######  Theorem 1.

Suppose Assumption 1, 2 hold, each $f_{i}$ is convex, $C=\operatorname{dom}f$,
$\gamma_{k}$ is independent of the sample $\nabla f_{k}(x^{k})$, and choose
$\gamma_{b,k}=\frac{\gamma_{b,0}}{\sqrt{k+1}}$. Then, the envelope-type step
size in (9) achieves the following rate,

$\displaystyle\mathbb{E}[f(\bar{x}^{K})-f(x^{*})]\leq\frac{\lVert
x^{0}-x^{*}\rVert^{2}}{2\gamma_{l,K-1}K}+\frac{\gamma_{b,0}^{2}G^{2}\log(K)}{2\gamma_{l,K-1}K},$

where $\gamma_{l,K-1}=\min\\{\omega,\frac{\gamma_{b,0}}{\sqrt{K}}\\}$ and
$\bar{x}^{K}=\frac{1}{K}\sum_{k=0}^{K}x^{k}$.

We were not able to give a convergence result that uses the same sample for
computing the step size and the gradient. However, we empirically observe that
the performance is very similar when using either one or two independent
samples per iteration (see Figure 2 and Appendix B). When two independent
samples $i_{k}$ and $j_{k}$ are used per iteration, the first computes the
gradient sample $\nabla f_{i_{k}}(x^{k})$, and the other computes the step-
size $\gamma_{k}$. For example, for SPSB this gives
$\gamma_{k}=\min\\{\frac{f_{j_{k}}(x^{k})-l^{*}_{j_{k}}}{c_{k}\lVert\nabla
f_{j_{k}}(x^{k})\rVert^{2}},\gamma_{b,k}\\}$. This type of assumption has been
used in several other works for analyzing adaptive step sizes [27, 44, 29].
Under this assumption, we specialize the results of Theorem 1 to SPSB and SLSB
, where
$\gamma_{l,K-1}=\min\\{\frac{1}{2cL_{\max}},\frac{\gamma_{b,0}}{\sqrt{K}}\\}$
and
$\gamma_{l,K-1}=\min\\{\frac{2(1-\bar{c})}{L_{\max}},\frac{\gamma_{b,0}}{\sqrt{K}}\\}$
respectively. Concretely, for $K\geq\gamma_{b,0}^{2}L_{\max}^{2}$, SLSB and
SPSB with $\gamma_{b,k}=\frac{\gamma_{b,0}}{\sqrt{k+1}}$ and
$c=\bar{c}=\frac{1}{2}$ achieve the following rate:
$\mathbb{E}[f(\bar{x}^{K})-f(x^{*})]\leq\frac{\lVert
x^{0}-x^{*}\rVert^{2}}{2\gamma_{b,0}\sqrt{K}}+\frac{\gamma_{b,0}G^{2}\log(K)}{2\sqrt{K}}.$
Next, we state the result for the envelop-type step size when $f$ is
$\mu$-strongly convex.

######  Theorem 2.

Suppose a $\mu$-strongly convex function $f$ satisfying Assumptions 1 and 2,
assume $\mathcal{C}$ is a closed and convex set, and $\gamma_{k}$ is
independent of the sample $\nabla f_{k}(x^{k})$. Then an envelope-type step
size as in (9) with $\gamma_{b,k}=\frac{\gamma_{b,0}}{k+1}$,
$\gamma_{b,0}\geq\frac{1}{\mu}$, and $\omega\mu<1$ achieves the following rate

$\displaystyle\mathbb{E}[f(\bar{x}_{K})-f(x^{*})]\leq\frac{\mu
k_{0}}{2(K-k_{0})}\bigl{(}e^{-k_{0}\mu\omega}\lVert
x_{0}-x^{*}\rVert^{2}+\gamma_{b,0}^{2}G^{2}\bigr{)}+\frac{\gamma_{b,0}G^{2}\log
K}{2(K-k_{0})},$

where $\bar{x}_{K}=\frac{1}{K-k_{0}}\sum_{k=k_{0}}^{K-1}x^{k}$ and
$k_{0}=\max\\{1,\lceil\gamma_{b,0}/\omega\rceil-1\\}$.

We can again apply the result of Theorem 2 to SPSB and SLSB with
$\gamma_{b,k}=\frac{\gamma_{b,0}}{k+1}$, $\gamma_{b,0}\geq\frac{1}{\mu}$,
$\omega=1/L_{\max}$, and $c=\bar{c}=\frac{1}{2}$ to get an explicit rate:
$\mathbb{E}[f(\bar{x}_{K})-f(x^{*})]\leq\frac{\mu
k_{0}}{2(K-k_{0})}\bigl{(}e^{\frac{-k_{0}\mu}{L_{\max}}}\lVert
x_{0}-x^{*}\rVert^{2}+\gamma_{b,0}^{2}G^{2}\bigr{)}+\frac{\gamma_{b,0}G^{2}\log
K}{2(K-k_{0})}$, where $k_{0}=\max\\{1,\lceil\gamma_{b,0}L_{\max}\rceil-1\\}$.

######  Remark 1.

Under the envelop-type step size framework and the assumption of two
independent samples, SLSB and SPSB share the same convergence rates of
$\mathcal{O}(\frac{1}{\sqrt{K}})$ and $\mathcal{O}(\frac{1}{K})$ as SGD with
decaying step-size for convex and strongly-convex losses respectively. This is
not surprising because of the structure of the envelope step-size in (9).
Indeed, the proof is similar to the standard proof of analogous rate for SGD
with decaying step-size. Nonetheless, we include it here for completeness.

###  3.2 Envelope-type step size for bi-level optimization

We start with recalling standard assumptions in BO [22, 19, 15, 4]. We denote
$z=[x;y]$ and recall the bi-level problem in (1). The first assumption is on
the lower-level objective $g$.

######  Assumption 3.

The function $g(x,y)$ is $\mu_{g}$ strongly convex in $y$ for any given $x$.
Moreover, $\nabla g$ is Lipschitz continuous: $\lVert\nabla
g(x_{1},y_{1})-\nabla g(x_{2},y_{2})\rVert\leq L_{g}\lVert z_{1}-z_{2}\rVert$
(also assume that this holds true for each sampled function $g(x,y;\psi)$),
and $\nabla^{2}g$ is Lipschitz continuous:
$\lVert\nabla^{2}g(x_{1},y_{1})-\nabla^{2}g(x_{2},y_{2})\rVert\leq L_{G}\lVert
z_{1}-z_{2}\rVert$. We further assume that
$\lVert\nabla^{2}_{xy}g(x,y)\rVert\leq C_{g}$, and the condition number is
defined as $\kappa=\frac{L_{g}}{\mu_{g}}$.

Next, we state the assumptions on the upper objective $f$.

######  Assumption 4.

The function $f$ and its gradients are Lipschitz continuous. That is: $\lVert
f(x_{1},y_{1})-f(x_{2},y_{2})\rVert\leq L_{1}\lVert z_{1}-z_{2}\rVert$ and
$\lVert\nabla f(x_{1},y_{1})-\nabla f(x_{2},y_{2})\rVert\leq L_{f,1}\lVert
z_{1}-z_{2}\rVert$. We also assume that $\lVert\nabla_{y}f(x,y)\rVert\leq
C_{f}$.

Furthermore, we make the following standard assumptions on the estimates of
$\nabla f$, $\nabla g$, and $\nabla^{2}g$.

######  Assumption 5.

The stochastic gradients are unbiased: $\mathbb{E}_{\phi}[\nabla
f(x,y;\phi)]=\nabla f(x,y)$, $\mathbb{E}_{\psi}[\nabla g(x,y;\psi)]=\nabla
g(x,y)$, and $\mathbb{E}_{\psi}[\nabla^{2}g(x,y;\psi)]=\nabla^{2}g(x,y)$. The
variances of $\nabla f(x,y;\phi)$ and $\nabla^{2}g(x,y;\psi])$ are bounded:
$\mathbb{E}_{\phi}[\lVert\nabla f(x,y;\phi)-\nabla
f(x,y)\rVert^{2}]\leq\sigma_{f}^{2}$ and
$\mathbb{E}_{\psi}[\lVert\nabla^{2}g(x,y;\psi)-\nabla^{2}g(x,y)\rVert^{2}]\leq\sigma_{G}^{2}$.

Finally, we introduce the bounded optimal function value assumption in (15),
which is used specifically for analyzing step size of the form (11) in the bi-
level setting:

$\displaystyle\mathbb{E}_{\psi}[g(x,y^{*}(x);\psi)-g(x,y^{*}_{x,\psi};\psi)]\leq\sigma_{g}^{2},\forall
x,$ (15)
$\displaystyle\mathbb{E}[\lVert\nabla_{y}g(x,y)-\nabla_{y}g(x,y;\phi)\rVert^{2}]\leq\sigma_{g}^{2},\forall
x,y,$ (16)

where $y^{*}(x)=\min_{y}g(x,y)$ and $y^{*}_{x,\psi}=\min_{y}g(x,y;\psi)$ for a
given $x$ (recall that at any iteration $k$, the lower-level steps in BiSPS
are $\text{SPS}_{\max}$ with an upper bound $\beta_{b,k}$; furthermore,
$\beta_{b,k}$ is non-increasing w.r.t. upper iteration $k$). The one-variable
analogous assumption of (15) has been used in the analysis of
$\text{SPS}_{\max}$ [29]. Here, we extend it to a two-variable function.
Unlike the bounded variance assumption (16), which needs to hold true for all
$x$ and $y$, we require (15) to hold at $y^{*}(x)$ for any given $x$. As
mentioned previously, the closed form solution $y^{*}(x)$ is difficult to
obtain. Thus, we define the following expression by replacing $y^{*}(x)$ with
$y$ in (2):

$\displaystyle\bar{\nabla}f(x,y)$
$\displaystyle=\nabla_{x}f(x,y)+\nabla_{xy}^{2}g(x,y)[\nabla_{yy}^{2}g(x,y)]^{-1}\nabla_{y}f(x,y).$
(17)

A stochastic Neumann series in (13) approximates (17) with $x$ and $y$ being
$x^{k}$ and $y^{k+1}$ (respectively), also recall that $y^{k+1}$ is an
approximation of $y^{*}(x^{k})$ by running $T$ lower-level SGD steps to
minimize $g$ w.r.t. $y$ for a fixed $x^{k}$. Based on Assumptions 3, 4, and 5,
we have the following results to be used in the analysis [15]: $\lVert\nabla
F(x_{1})-\nabla F(x_{2})\rVert\leq L_{F}\lVert x_{1}-x_{2}\rVert$, $\lVert
y^{*}(x_{1})-y^{*}(x_{2})\rVert\leq L_{y}\lVert x_{1}-x_{2}\rVert$, and
$\lVert\bar{\nabla}f(x,y^{*}(x))-\bar{\nabla}f(x,y)\rVert\leq L_{f}\lVert
y^{*}(x)-y\rVert$. Furthermore, the bias in the stochastic hypergradient in
(13) (denoted as $B$) decays exponentially with $N$ and its variance is
bounded, i.e. $\mathbb{E}[\lVert
h_{f}^{k}-\mathbb{E}[h_{f}^{k}]\rVert^{2}]\leq\tilde{\sigma}_{f}^{2}$ (see
Appendix A for details) [19].

Now, we state our main theorem based on step size of the form (10) and (11).

######  Theorem 3.

Suppose $f$ and $g$ satisfy Assumptions 3, 4, and 5, learning-rate upper
bounds $\alpha_{b,k}=\frac{\alpha_{b,0}}{\sqrt{k+1}}$ and
$\alpha_{l,k}=\frac{\alpha_{l,0}}{\sqrt{k+1}}$ with $\alpha_{b,0}$ and
$\alpha_{l,0}$ satisfying
$\frac{1}{L_{F}+4L_{y}^{2}}\geq\frac{\alpha_{b,0}^{2}}{\alpha_{l,0}}$ and
$\alpha_{l,0}\leq\alpha_{b,0}$. Further assume that $\alpha_{k}$ is
independent of the stochastic hypergradient $h_{f}^{k}$, and each sampled
function $g(x,y;\psi)$ is convex. Then under the Assumption (15) with
$p\geq\frac{1}{2}$, $C_{k}=\min\\{\frac{1}{2pL_{g}},\beta_{b,k}\\}$,
$T\geq\frac{\log(\alpha_{b,0}L_{f}^{2}+2)}{-\log(1-\mu_{g}C_{K-1})}$, and
$\beta_{b,k}=\frac{\beta_{b,0}}{k+1}$, BiSPS achieves the rate:

$\displaystyle\frac{1}{K}\sum_{k=0}^{K-1}\mathbb{E}[\lVert\nabla
F(x^{k})\rVert^{2}]$
$\displaystyle\leq\mathcal{\tilde{O}}(\frac{\kappa^{3}}{\sqrt{K}}+\frac{\kappa^{2}\log
K}{\sqrt{K}}).$ (18)

######  Remark 2.

We further give the convergence result under the bounded variance assumption
(16) in Appendix A. Theorem 3 shows that BiSPS matches the optimal rate of SGD
up to a logarithmic factor without a growing batch size. We notice that the
assumption (15) largely simplifies the expression on $T$ and does not require
an explicit upper bound on $\beta_{b,0}$. As in the single-level case, whether
using one sample or two samples (which makes upper-level step-size independent
of gradient) gives similar empirical performances (see Appendix B). Note that
the independence assumption is only needed for the upper-level. Thus, the two-
sample requirement of theorem does not apply to the lower-level problem. This
is useful from computational standpoint as typical bi-level algorithms run
multiple lower-level updates for each upper-level iteration.

##  4 Additional Hyper-Representation and Data Distillation Experiments

Hyper-representation learning: The experiments are performed on MNIST dataset
using LeNet [26, 42]. We use conjugate gradient method for solving system of
equations when computing the hypergradient [17]. The upper and lower-level
objectives are to optimize the embedding layers and the classifier (i.e. the
last layer of the neural net), respectively (see Appendix B for details). For
constant-step SGD and Adam, we tune the lower-level learning rate
$\beta\in\\{10.0,5.0,1.0,0.5,0.1,0.05,0.01\\}$. For the upper-level learning
rate, we tune $\alpha\in\\{0.001,0.0025,0.005,0.01,0.05,0.1\\}$ for SGD, and
$\alpha\in\\{10^{-5},5\cdot 10^{-5},10^{-4},5\cdot 10^{-4},0.001,0.01\\}$ for
Adam (recall that $\delta$ in (14) is set to 0). Based on the results of
Figure 6, we make the following key observations: ① line-search at the upper-
level is essential for achieving the optimal performance (Figure 6(a)); ②
BiSLS-Adam/SGD not only converges fast but also generalizes well (Figure
6(b)); ③ BiSLS-Adam/SGD is highly robust to search starting points
$\alpha_{b,0}$ and $\beta_{b,0}$ (Figure 6(c), 6(d)). It addresses the
fundamental question of how to tune $\alpha$ and $\beta$ in bi-level
optimization (see Appendix B for additional results on search cost).

(a) Validation Loss

(b) Validation Accuracy

(c) Different $\alpha_{b,0}$

(d) Different $\beta_{b,0}$

Figure 6: Validation loss (a) and accuracy (b) against iterations. (a)
Comparisons between whether to use or not use line-search at the upper or
lower level; (b) Generalization performance of BiSLS-Adam/SGD and fine-tuned
Adam/SGD; (c) Validation loss against iterations for different values of
$\alpha_{b,0}$ ($\beta_{b,0}$ fixed at $100$). (d) Same plot as (c) but for
different values of $\beta_{b,0}$ ($\alpha_{b,0}$ fixed at $10$).

Data distillation: The goal of data distillation is to generate a small set of
synthetic data from an original dataset that preserves the performance of a
model when trained using the generated data [46, 50]. We adapted the
experiment set up from Lorraine et al. [30] to distill MNIST digits. We
present the results in Figure 7, where we observe that BiSLS-SGD converges
significantly faster than fine-tuned Adam or SGD, and generate realistic MNIST
images (see Appendix B for more results).

(a) Neumann

(b) Identity

(c) Distilled MNIST Images

Figure 7: (a)(b): Comparison between BiSLS-SGD and Adam/SGD for Data
Distillation on MNIST dataset. Validation loss plotted against iterations. (a)
Hypergradient computed using Neumann series; (b) Inverse Hessian in (2)
treated as the Identity [30] when computing the hypergradient; (c) Distilled
MNIST images after 3000 iterations of BiSLS-SGD.

##  5 Conclusion

In this work, we have given simple alternatives to SLS and SPS that show good
empirical performance in non-interpolating scenario without requiring the step
size to be monotonic. We unify their analysis based on a simplified envelope-
type step size, and extend the analysis to the bi-level setting while
designing a SPS-based bi-level algorithm. In the end, we propose bi-level
line-search algorithm BiSLS-Adam/SGD that is empirically truly robust and
adaptive to learning rate initialization. Our work opens several possible
future directions. Given the superior performance of BiSLS, we prioritize an
analysis of its convergence rates. The difficulty stems from: (a) the bias in
hypergradient estimation; (b) the dual updates in $x$ and $y^{*}(x)$
(incurring nested loop structures); (c) the error in estimating $y^{*}(x)$. On
single-level optimization, we remark as an important direction to relax the
two-sample assumption on SPSB /SLSB . On the other hand, we also deem
interesting that the two-sample variants of our algorithms do not appear to
deteriorate performance in practice. we hope our ideas motivate further
research on practical bi-level optimization algorithms that avoid the often
computationally exhausting need for tuning step-sizes. At the same time, we
hope to see more results alike extending the SPS/SLS beyond the originally
investigated setting of interpolating single-level optimization.

##### Acknowledgement

This work was partially supported by the Natural Sciences and Engineering
Research Council of Canada (NSERC) Discovery Grants RGPIN-2021-03677 and
GPIN-2022-03669.

## References

  * Beck [2017] Amir Beck.  _First-order methods in optimization_.  SIAM, 2017.
  * Borsos et al. [2020] Zalán Borsos, Mojmir Mutny, and Andreas Krause.  Coresets via bilevel optimization for continual learning and streaming.  _Advances in Neural Information Processing Systems_ , 33:14879–14890, 2020.
  * Chang and Lin [2011] Chih-Chung Chang and Chih-Jen Lin.  Libsvm: a library for support vector machines.  _ACM transactions on intelligent systems and technology (TIST)_ , 2(3):1–27, 2011.
  * Chen et al. [2021] Tianyi Chen, Yuejiao Sun, and Wotao Yin.  Tighter analysis of alternating stochastic gradient method for stochastic nested problems, 2021.
  * Chen et al. [2022] Tianyi Chen, Yuejiao Sun, Quan Xiao, and Wotao Yin.  A single-timescale method for stochastic bilevel optimization, 2022.
  * Cutkosky and Orabona [2020] Ashok Cutkosky and Francesco Orabona.  Momentum-based variance reduction in non-convex sgd, 2020.
  * Dagréou et al. [2022] Mathieu Dagréou, Pierre Ablin, Samuel Vaiter, and Thomas Moreau.  A framework for bilevel optimization that enables stochastic and global variance reduction algorithms, 2022.
  * Dagréou et al. [2023] Mathieu Dagréou, Thomas Moreau, Samuel Vaiter, and Pierre Ablin.  A lower bound and a near-optimal algorithm for bilevel empirical risk minimization, 2023.
  * Defazio et al. [2014] Aaron Defazio, Francis Bach, and Simon Lacoste-Julien.  Saga: A fast incremental gradient method with support for non-strongly convex composite objectives, 2014.
  * Falk and Liu [1995] James E Falk and Jiming Liu.  On bilevel programming, part i: general nonlinear cases.  _Mathematical Programming_ , 70:47–72, 1995.
  * Fan et al. [2021] Chen Fan, Parikshit Ram, and Sijia Liu.  Sign-maml: Efficient model-agnostic meta-learning by signsgd, 2021.
  * Fang et al. [2018] Cong Fang, Chris Junchi Li, Zhouchen Lin, and Tong Zhang.  Spider: Near-optimal non-convex optimization via stochastic path integrated differential estimator, 2018.
  * Finn et al. [2017] Chelsea Finn, Pieter Abbeel, and Sergey Levine.  Model-agnostic meta-learning for fast adaptation of deep networks, 2017\.
  * Franceschi et al. [2017] Luca Franceschi, Michele Donini, Paolo Frasconi, and Massimiliano Pontil.  Forward and reverse gradient-based hyperparameter optimization, 2017.
  * Ghadimi and Wang [2018] Saeed Ghadimi and Mengdi Wang.  Approximation methods for bilevel programming, 2018.
  * Gower et al. [2019] Robert Mansel Gower, Nicolas Loizou, Xun Qian, Alibek Sailanbayev, Egor Shulgin, and Peter Richtárik.  Sgd: General analysis and improved rates.  In _International conference on machine learning_ , pages 5200–5209. PMLR, 2019.
  * Grazzi et al. [2020] Riccardo Grazzi, Luca Franceschi, Massimiliano Pontil, and Saverio Salzo.  On the iteration complexity of hypergradient computation, 2020.
  * Grazzi et al. [2021] Riccardo Grazzi, Massimiliano Pontil, and Saverio Salzo.  Convergence properties of stochastic hypergradients, 2021.
  * Hong et al. [2022] Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang.  A two-timescale framework for bilevel optimization: Complexity analysis and application to actor-critic, 2022.
  * Huang et al. [2023] Feihu Huang, Junyi Li, and Shangqian Gao.  Biadam: Fast adaptive bilevel optimization methods, 2023.
  * Ishizuka and Aiyoshi [1992] Yo Ishizuka and Eitaro Aiyoshi.  Double penalty method for bilevel optimization problems.  _Annals of Operations Research_ , 34(1):73–88, 1992.
  * Ji et al. [2021] Kaiyi Ji, Junjie Yang, and Yingbin Liang.  Bilevel optimization: Convergence analysis and enhanced design, 2021.
  * Ji et al. [2022] Kaiyi Ji, Mingrui Liu, Yingbin Liang, and Lei Ying.  Will bilevel optimizers benefit from loops, 2022.
  * Khanduri et al. [2021] Prashant Khanduri, Siliang Zeng, Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang.  A near-optimal algorithm for stochastic bilevel optimization via double-momentum, 2021.
  * Kingma and Ba [2014] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * LeCun et al. [1998] Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner.  Gradient-based learning applied to document recognition.  _Proceedings of the IEEE_ , 86(11):2278–2324, 1998.
  * Li and Orabona [2019] Xiaoyu Li and Francesco Orabona.  On the convergence of stochastic gradient descent with adaptive stepsizes, 2019.
  * Liu et al. [2019] Hanxiao Liu, Karen Simonyan, and Yiming Yang.  Darts: Differentiable architecture search, 2019.
  * Loizou et al. [2021] Nicolas Loizou, Sharan Vaswani, Issam Laradji, and Simon Lacoste-Julien.  Stochastic polyak step-size for sgd: An adaptive learning rate for fast convergence, 2021.
  * Lorraine et al. [2020] Jonathan Lorraine, Paul Vicol, and David Duvenaud.  Optimizing millions of hyperparameters by implicit differentiation.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1540–1552. PMLR, 2020.
  * Loshchilov and Hutter [2019] Ilya Loshchilov and Frank Hutter.  Decoupled weight decay regularization, 2019.
  * Luo et al. [2019] Liangchen Luo, Yuanhao Xiong, Yan Liu, and Xu Sun.  Adaptive gradient methods with dynamic bound of learning rate, 2019.
  * Nemirovski et al. [2009] Arkadi Nemirovski, Anatoli Juditsky, Guanghui Lan, and Alexander Shapiro.  Robust stochastic approximation approach to stochastic programming.  _SIAM Journal on optimization_ , 19(4):1574–1609, 2009.
  * Nesterov [2003] Yurii Nesterov.  _Introductory lectures on convex optimization: A basic course_ , volume 87.  Springer Science & Business Media, 2003.
  * Nguyen et al. [2017] Lam M Nguyen, Jie Liu, Katya Scheinberg, and Martin Takáč.  Sarah: A novel method for machine learning problems using stochastic recursive gradient.  In _International Conference on Machine Learning_ , pages 2613–2621. PMLR, 2017.
  * Orvieto et al. [2022] Antonio Orvieto, Simon Lacoste-Julien, and Nicolas Loizou.  Dynamics of sgd with stochastic polyak stepsizes: Truly adaptive variants and convergence to exact solution, 2022.
  * Rajeswaran et al. [2019] Aravind Rajeswaran, Chelsea Finn, Sham Kakade, and Sergey Levine.  Meta-learning with implicit gradients, 2019.
  * Reddi et al. [2019] Sashank J Reddi, Satyen Kale, and Sanjiv Kumar.  On the convergence of adam and beyond.  _arXiv preprint arXiv:1904.09237_ , 2019.
  * Ren et al. [2021] Pengzhen Ren, Yun Xiao, Xiaojun Chang, Po-Yao Huang, Zhihui Li, Xiaojiang Chen, and Xin Wang.  A comprehensive survey of neural architecture search: Challenges and solutions.  _ACM Computing Surveys (CSUR)_ , 54(4):1–34, 2021\.
  * Shaban et al. [2019] Amirreza Shaban, Ching-An Cheng, Nathan Hatch, and Byron Boots.  Truncated back-propagation for bilevel optimization.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 1723–1732. PMLR, 2019.
  * Shalev-Shwartz et al. [2007] Shai Shalev-Shwartz, Yoram Singer, and Nathan Srebro.  Pegasos: Primal estimated sub-gradient solver for svm.  In _Proceedings of the 24th international conference on Machine learning_ , pages 807–814, 2007.
  * Sow et al. [2022] Daouda Sow, Kaiyi Ji, and Yingbin Liang.  On the convergence theory for hessian-free bilevel algorithms, 2022.
  * Vaswani et al. [2020] Sharan Vaswani, Issam Laradji, Frederik Kunstner, Si Yi Meng, Mark Schmidt, and Simon Lacoste-Julien.  Adaptive gradient methods converge faster with over-parameterization (but you should do a line-search).  _arXiv preprint arXiv:2006.06835_ , 2020.
  * Vaswani et al. [2021] Sharan Vaswani, Aaron Mishkin, Issam Laradji, Mark Schmidt, Gauthier Gidel, and Simon Lacoste-Julien.  Painless stochastic gradient: Interpolation, line-search, and convergence rates, 2021.
  * Vicente et al. [1994] Luis Vicente, Gilles Savard, and Joaquim Júdice.  Descent approaches for quadratic bilevel programming.  _Journal of optimization theory and applications_ , 81(2):379–399, 1994.
  * Wang et al. [2020] Tongzhou Wang, Jun-Yan Zhu, Antonio Torralba, and Alexei A. Efros.  Dataset distillation, 2020.
  * Ward et al. [2021] Rachel Ward, Xiaoxia Wu, and Leon Bottou.  Adagrad stepsizes: Sharp convergence over nonconvex landscapes, 2021.
  * Wright [2006] Jorge Nocedal Stephen J Wright.  Numerical optimization, 2006.
  * Yang et al. [2021] Junjie Yang, Kaiyi Ji, and Yingbin Liang.  Provably faster algorithms for bilevel optimization, 2021.
  * Yu et al. [2023] Ruonan Yu, Songhua Liu, and Xinchao Wang.  Dataset distillation: A comprehensive review, 2023.
  * Zhang et al. [2022] Yihua Zhang, Guanhua Zhang, Prashant Khanduri, Mingyi Hong, Shiyu Chang, and Sijia Liu.  Revisiting and advancing fast adversarial training through the lens of bi-level optimization.  In _International Conference on Machine Learning_ , pages 26693–26712. PMLR, 2022.
  * Zhang et al. [2023] Yihua Zhang, Yuguang Yao, Parikshit Ram, Pu Zhao, Tianlong Chen, Mingyi Hong, Yanzhi Wang, and Sijia Liu.  Advancing model pruning via bi-level optimization, 2023.
  * Zhou et al. [2022] Yongchao Zhou, Ehsan Nezhadarya, and Jimmy Ba.  Dataset distillation using neural feature regression.  _arXiv preprint arXiv:2206.00719_ , 2022.

Notes: A few small typos have been revised in the main body. The $>$ sign in
(14) is corrected to $\leq$, and the missing $A_{k}^{-1}$ for BiSLS-Adam in
(14) is added. In Theorem 3 and the remark that follows, the letter $\gamma$
is corrected to letter $\beta$. In Assumption 3, $L_{f}$ is changed to
$L_{f,1}$.

##  Appendix A Proofs of Theorems and Additional Convergence Results

###  A.1 Useful Lemmas

Lemma 1 provides more details on the envelope structure of SPSB and SLSB given
in (10) and (11). The lower bound in (19) will also be used in Lemma 5 (for
bounding the term $\lVert y^{k+1}-y^{*}(x^{k})\rVert^{2}$).

######  Lemma 1.

Under the Assumption 1, we have the following:

$\displaystyle\text{SPSB~{}:
}\quad\min\\{\frac{1}{2cL_{\max}},\gamma_{b,k}\\}\leq\gamma_{k}=\min\\{\frac{f_{i_{k}}(x^{k})-l_{i_{k}}^{*}}{c\lVert\nabla
f_{i_{k}}(x^{k})\rVert^{2}},\gamma_{b,k}\\},\quad 0<c,$ (19)
$\displaystyle\text{SLSB~{}:
}\quad\min\\{\frac{2(1-\bar{c})}{L_{\max}},\gamma_{b,k}\\}\leq\gamma_{k}\leq\min\\{\frac{f_{i_{k}}(x^{k})-l_{i_{k}}^{*}}{\bar{c}\lVert\nabla
f_{i_{k}}(x^{k})\rVert^{2}},\gamma_{b,k}\\},\quad 0<\bar{c}<1.$ (20)

###### Proof.

The bounds in (19) have been shown in [29, 36]. For (20), the first part of
the inequality has been shown in [44]. For the second part, recall the Armijo
condition (14):

$\displaystyle f_{i_{k}}(x_{k}-\gamma_{k}\nabla f_{i_{k}}(x_{k}))\leq
f_{i_{k}}(x_{k})-\bar{c}\cdot\gamma_{k}\lVert\nabla
f_{i_{k}}(x_{k})\rVert^{2},\quad 0<\bar{c}<1.$

We can then rearrange this to obtain

$\displaystyle\gamma_{k}\leq\frac{f_{i_{k}}(x_{k})-f_{i_{k}}(x_{k}-\gamma_{k}\nabla
f_{i_{k}}(x_{k}))}{\bar{c}\lVert\nabla
f_{i_{k}}(x^{k})\rVert^{2}}\leq\frac{f_{i_{k}}(x_{k})-f_{i_{k}}^{*}}{\bar{c}\lVert\nabla
f_{i_{k}}(x^{k})\rVert^{2}}\leq\frac{f_{i_{k}}(x_{k})-l_{i_{k}}^{*}}{\bar{c}\lVert\nabla
f_{i_{k}}(x^{k})\rVert^{2}},$ (21)

where $l_{i_{k}}^{*}$ is any lower bound for $f_{i_{k}}^{*}$. Also recall that
$\gamma_{b,k}$ is the search starting point at iteration $k$, hence (20) holds
for SLSB . ∎

Lemma 2 gives the expressions for the constants $L_{F}$, $L_{y}$, and $L_{f}$.
The proof can be found in [15, Lem 2.2].

######  Lemma 2.

Under Assumptions 3 and 4, we have the following:

$\displaystyle\lVert\nabla F(x_{1})-\nabla F(x_{2})\rVert\leq L_{F}\lVert
x_{1}-x_{2}\rVert,$ $\displaystyle\lVert y^{*}(x_{1})-y^{*}(x_{2})\rVert\leq
L_{y}\lVert x_{1}-x_{2}\rVert,$
$\displaystyle\lVert\bar{\nabla}f(x,y^{*}(x))-\bar{\nabla}f(x,y)\rVert\leq
L_{f}\lVert y^{*}(x)-y\rVert,$

where

$\displaystyle
L_{f}=L_{f,1}+\frac{L_{f,1}L_{g}}{\mu_{g}}+\frac{L_{1}}{\mu_{g}}(L_{G}+\frac{L_{g}L_{G}}{\mu_{g}})\sim\mathcal{O}(\kappa^{2})$
$\displaystyle L_{y}=\frac{L_{g}}{\mu_{g}}\sim\mathcal{O}(\kappa)$
$\displaystyle
L_{F}=L_{f,1}+\frac{L_{g}(L_{f,1}+L_{f})}{\mu_{g}}+\frac{L_{1}}{\mu_{g}}(L_{G}+\frac{L_{g}L_{G}}{\mu_{g}})\sim\mathcal{O}(\kappa^{3})$

Lemma 3 is on the bias and variance of the stochastic hypergradient in (13),
which has the following form [19]

$\displaystyle
h_{f}^{k}=\nabla_{x}f(x^{k},y^{k+1};\phi)-\nabla_{xy}g(x^{k},y^{k+1};\psi_{0})\bigl{[}\frac{N}{L_{g}}\prod_{j=1}^{\bar{N}}(I-\nabla^{2}_{yy}g(x^{k},y^{k+1};\psi_{j}))\bigr{]}\nabla_{y}f(x^{k},y^{k+1};\phi).$
(22)

Recall that the hypergradient surrogate defined in (17) based on
$(x^{k},y^{k+1})$ is

$\displaystyle\bar{\nabla}f(x^{k},y^{k+1})$
$\displaystyle=\nabla_{x}f(x^{k},y^{k+1})+\nabla_{xy}^{2}g(x^{k},y^{k+1})[\nabla_{yy}^{2}g(x^{k},y^{k+1})]^{-1}\nabla_{y}f(x^{k},y^{k+1}).$
(23)

Given a filtration $\mathcal{F}^{{}^{\prime}}_{k}$ up to and including $x^{k}$
and $y^{k+1}$, the bias of the stochastic hypergradient is defined as
$B=\bar{\nabla}f(x^{k},y^{k+1})-\mathbb{E}[h_{f}^{k}|\mathcal{F}^{{}^{\prime}}_{k}]$,
and the variance is defined as
$\mathbb{E}[\lVert\bar{\nabla}f(x^{k},y^{k+1})-\mathbb{E}[h_{f}^{k}|\mathcal{F}^{{}^{\prime}}_{k}]\rVert^{2}]$.
Lemma 3 has been proven in [19, Lem 1.] (also see [4, Lem 5.]). Lemmas 1, 2,
and 3 will be used in the proofs of Theorems 3 and 4.

######  Lemma 3.

Under Assumptions 3, 4, and 5, the bias and variance of the stochastic
hypergradient $h_{f}^{k}$ satisfy the following

$\displaystyle\text{Bias:
}\lVert\bar{\nabla}f(x^{k},y^{k+1})-\mathbb{E}[h_{f}^{k}|\mathcal{F}^{{}^{\prime}}_{k}]\rVert\leq\frac{C_{g}C_{f}}{\mu_{g}}(1-\frac{\mu_{g}}{L_{g}})^{N},\forall
k$ $\displaystyle\text{Variance:
}\mathbb{E}[\lVert\bar{\nabla}f(x^{k},y^{k+1})-\mathbb{E}[h_{f}^{k}|\mathcal{F}^{{}^{\prime}}_{k}]\rVert^{2}]\leq\tilde{\sigma}_{f}^{2},\forall
k,$

where N is the total number of samples, and
$\tilde{\sigma}_{f}^{2}=\sigma_{f}^{2}+[(\sigma_{f}^{2}+L_{1}^{2})(\sigma_{G}^{2}+2L_{g}^{2})+\sigma_{f}^{2}L_{g}^{2}]\frac{3}{\mu_{g}^{2}}\sim\mathcal{O}(\kappa^{2})$.

Lemma 4 is on the descent of the quantity $V^{k}:=F(x^{k})-F(x^{*})+\lVert
y^{k}-y^{*}(x^{k})\rVert^{2}$. It will be used in the proofs of Theorem 3 and
4.

######  Lemma 4.

Suppose $F$ satisfies Assumptions 3, 4, and 5, sequences
$\alpha_{b,k}=\frac{\alpha_{b,0}}{\sqrt{k+1}}$ and
$\alpha_{l,k}=\frac{\alpha_{l,0}}{\sqrt{k+1}}$ with $\alpha_{b,0}$ and
$\alpha_{l,0}$ satisfying
$\frac{1}{L_{F}+4L_{y}^{2}}\geq\frac{\alpha_{b,0}^{2}}{\alpha_{l,0}}$ and
$\alpha_{l,0}\leq\alpha_{b,0}$. Further assume that $\alpha_{k}$ is
independent of $h_{f}^{k}$. Then step sizes of the form (10) and (11) achieve
the following:

$\displaystyle E[V^{k+1}]$
$\displaystyle\leq\mathbb{E}[V^{k}]-\frac{\alpha_{l,k}}{2}\mathbb{E}[\lVert\nabla
F(x^{k})\rVert^{2}]+(\alpha_{b,k}L_{f}^{2}+2)\mathbb{E}[\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}]+$
$\displaystyle\quad\quad\alpha_{b,k}B^{2}+(2L_{y}^{2}\alpha_{b,k}^{2}+\frac{L_{F}\alpha_{b,k}^{2}}{2})\tilde{\sigma}_{f}^{2}-\mathbb{E}[\lVert
y^{k}-y^{*}(x^{k})\rVert^{2}],\forall k,$ (24)

where $V^{k}=F(x^{k})-F(x^{*})+\lVert y^{k+1}-y^{*}(x^{k+1})\rVert^{2}$ and
recall that $F$ is the upper-level loss defined in (1).

###### Proof.

We denote
$\mathbb{E}[h_{f}^{k}|\mathcal{F}^{{}^{\prime}}_{k}]=\bar{h}_{f}^{k}$. By the
$L_{F}$-smoothness of the objective $F$:

$\displaystyle F(x^{k+1})\leq F(x^{k})+\langle\nabla
F(x^{k}),x^{k+1}-x^{k}\rangle+\frac{L_{F}}{2}\lVert x^{k+1}-x^{k}\rVert^{2}.$

Take expectation conditioned on a filtration of past iterates
$\mathcal{F}_{k}^{{}^{\prime}}$ (up to and include $x^{k}$, $y^{k+1}$):

$\displaystyle\mathbb{E}[F(x^{k+1})|\mathcal{F}_{k}^{{}^{\prime}}]$
$\displaystyle\leq F(x^{k})+\mathbb{E}[\langle\nabla
F(x^{k}),x^{k+1}-x^{k}\rangle|\mathcal{F}_{k}^{{}^{\prime}}]+\frac{L_{F}}{2}\mathbb{E}[\lVert
x^{k+1}-x^{k}\rVert^{2}|\mathcal{F}_{k}^{{}^{\prime}}]$
$\displaystyle=F(x^{k})-\mathbb{E}[\alpha_{k}\langle\nabla
F(x^{k}),h_{f}^{k}\rangle|\mathcal{F}_{k}^{{}^{\prime}}]+\frac{L_{F}\mathbb{E}[\alpha_{k}^{2}|\mathcal{F}_{k}^{{}^{\prime}}]}{2}\mathbb{E}[\lVert
h_{f}^{k}-\bar{h}_{f}^{k}+\bar{h}_{f}^{k}\rVert^{2}|\mathcal{F}_{k}^{{}^{\prime}}]$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}F(x^{k})-\mathbb{E}[\alpha_{k}|\mathcal{F}_{k}^{{}^{\prime}}]\langle\nabla
F(x^{k}),\bar{h}_{f}^{k}\rangle+\frac{L_{F}\mathbb{E}[\alpha_{k}^{2}|\mathcal{F}_{k}^{{}^{\prime}}]}{2}\lVert\bar{h}_{f}^{k}\rVert^{2}+\frac{L_{F}\mathbb{E}[\alpha_{k}^{2}|\mathcal{F}_{k}^{{}^{\prime}}]}{2}\tilde{\sigma}_{f}^{2}$
$\displaystyle=F(x^{k})-\frac{\mathbb{E}[\alpha_{k}|\mathcal{F}_{k}^{{}^{\prime}}]}{2}\lVert\nabla
F(x^{k})\rVert^{2}-\frac{\mathbb{E}[\alpha_{k}|\mathcal{F}_{k}^{{}^{\prime}}]}{2}\lVert\bar{h}_{f}^{k}\rVert^{2}+\frac{\mathbb{E}[\alpha_{k}|\mathcal{F}_{k}^{{}^{\prime}}]}{2}\lVert\nabla
F(x^{k})-\bar{h}_{f}^{k}\rVert^{2}+$
$\displaystyle\quad\quad\frac{L_{F}\mathbb{E}[\alpha_{k}^{2}|\mathcal{F}_{k}^{{}^{\prime}}]}{2}\lVert\bar{h}_{f}^{k}\rVert^{2}+\frac{L_{F}\mathbb{E}[\alpha_{k}^{2}|\mathcal{F}_{k}^{{}^{\prime}}]}{2}\tilde{\sigma}_{f}^{2},$

where (a) is by the assumption that $\alpha_{k}$ is independent of
$h_{f}^{k}$. Then expand the term $\lVert\nabla
F(x^{k})-\bar{h}_{f}^{k}\rVert^{2}$ as follows:

$\displaystyle\lVert\nabla F(x^{k})-\bar{h}_{f}^{k}\rVert^{2}$
$\displaystyle=\lVert\nabla
F(x^{k})-\bar{\nabla}f(x^{k},y^{k+1})+\bar{\nabla}f(x^{k},y^{k+1})-\bar{h}_{f}^{k}\rVert^{2}$
$\displaystyle\leq 2\lVert\nabla
F(x^{k})-\bar{\nabla}f(x^{k},y^{k+1})\rVert^{2}+2\lVert\bar{\nabla}f(x^{k},y^{k+1})-\bar{h}_{f}^{k}\rVert^{2}$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}2L_{f}^{2}\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}+2B^{2},$

where (b) is by Lemma 2 and 3. Substituting this into the above:

$\displaystyle\mathbb{E}[F(x^{k+1})|F_{k}^{{}^{\prime}}]$ $\displaystyle\leq
F(x^{k})-\frac{\mathbb{E}[\alpha_{k}|F_{k}^{{}^{\prime}}]}{2}\lVert\nabla
F(x^{k})\rVert^{2}-\frac{\mathbb{E}[\alpha_{k}|F_{k}^{{}^{\prime}}]}{2}\lVert\bar{h}_{f}^{k}\rVert^{2}+\frac{L_{F}\mathbb{E}[\alpha_{k}^{2}|F_{k}^{{}^{\prime}}]}{2}\lVert\bar{h}_{f}^{k}\rVert^{2}+$
$\displaystyle\quad\quad
E[\alpha_{k}|\mathcal{F}_{k}^{{}^{\prime}}]L_{f}^{2}\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}+E[\alpha_{k}|\mathcal{F}_{k}^{{}^{\prime}}]B^{2}+\frac{L_{F}\mathbb{E}[\alpha_{k}^{2}|\mathcal{F}_{k}^{{}^{\prime}}]}{2}\tilde{\sigma}_{f}^{2}$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}F(x^{k})-\frac{\alpha_{l,k}}{2}\lVert\nabla
F(x^{k})\rVert^{2}-\frac{\alpha_{l,k}}{2}\lVert\bar{h}_{f}^{k}\rVert^{2}+\frac{L_{F}\alpha_{b,k}^{2}}{2}\lVert\bar{h}_{f}^{k}\rVert^{2}+$
$\displaystyle\quad\quad\alpha_{b,k}L_{f}^{2}\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}+\alpha_{b,k}B^{2}+\frac{L_{F}\alpha_{b,k}^{2}}{2}\tilde{\sigma}_{f}^{2},$

where (c) is by $\alpha_{l,k}\leq\mathbb{\alpha}_{k}$ and
$\alpha_{b,k}\geq\alpha_{k}$. Then take total expectations and subtract
$F(x^{*})$:

$\displaystyle\mathbb{E}[F(x^{k+1})-F(x^{*})]$
$\displaystyle\leq\mathbb{E}[F(x^{k})-F(x^{*})]-\frac{\alpha_{l,k}}{2}\mathbb{E}[\lVert\nabla
F(x^{k})\rVert^{2}]-\frac{\alpha_{l,k}}{2}\mathbb{E}[\lVert\bar{h}_{f}^{k}\rVert^{2}]+\frac{L_{F}\alpha_{b,k}^{2}}{2}\mathbb{E}[\lVert\bar{h}_{f}^{k}\rVert^{2}]+$
$\displaystyle\quad\quad\alpha_{b,k}L_{f}^{2}\mathbb{E}[\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}]+\alpha_{b,k}B^{2}+\frac{L_{F}\alpha_{b,k}^{2}}{2}\tilde{\sigma}_{f}^{2}$
(25)

Now define the potential function
$V^{k}:=F(x^{k})-F(x^{*})+\mathbb{\lVert}y^{k+1}-y^{*}(x^{k+1})\rVert^{2}$ and
expand the term $\lVert y^{k+1}-y^{*}(x^{k+1})\rVert^{2}$ as follows:

$\displaystyle\lVert y^{k+1}-y^{*}(x^{k+1})\rVert^{2}$ $\displaystyle=\lVert
y^{k+1}-y^{*}(x^{k})+y^{*}(x^{k})-y^{*}(x^{k+1})\rVert^{2}$ $\displaystyle\leq
2\lVert y^{k+1}-y^{*}(x^{k})\rVert^{2}+2\lVert
y^{*}(x^{k})-y^{*}(x^{k+1})\rVert^{2}$
$\displaystyle\stackrel{{\scriptstyle(d)}}{{\leq}}2\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}+2L_{y}^{2}\lVert x^{k+1}-x^{k}\rVert^{2}$
$\displaystyle=2\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}+2L_{y}^{2}\alpha_{k}^{2}\lVert
h_{f}^{k}\rVert^{2}$ $\displaystyle=2\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}+2L_{y}^{2}\alpha_{k}^{2}\lVert
h_{f}^{k}-\bar{h}_{f}^{k}+\bar{h}_{f}^{k}\rVert^{2},$

where (d) is by Lemma 2. Take expectation conditioned on
$\mathcal{F}_{k}^{{}^{\prime}}$:

$\displaystyle\mathbb{E}[\lVert
y^{k+1}-y^{*}(x^{k+1})\rVert^{2}|\mathcal{F}_{k}^{{}^{\prime}}]$
$\displaystyle\leq 2\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}+2L_{y}^{2}\alpha_{k}^{2}\lVert\bar{h}_{f}^{k}\rVert^{2}+2L_{y}^{2}\alpha_{k}^{2}\tilde{\sigma}_{f}^{2}$
$\displaystyle\leq 2\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}+2L_{y}^{2}\alpha_{b,k}^{2}\lVert\bar{h}_{f}^{k}\rVert^{2}+2L_{y}^{2}\alpha_{b,k}^{2}\tilde{\sigma}_{f}^{2}.$

Then, take total expectations:

$\displaystyle\mathbb{E}[\lVert y^{k+1}-y^{*}(x^{k+1})\rVert^{2}]\leq
2\mathbb{E}[\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}]+2L_{y}^{2}\alpha_{b,k}^{2}\mathbb{E}[\lVert\bar{h}_{f}^{k}\rVert^{2}]+2L_{y}^{2}\alpha_{b,k}^{2}\tilde{\sigma}_{f}^{2}.$
(26)

Now, based on the definition of $V^{k}$ and combining (25) and (26):

$\displaystyle E[V^{k+1}]$
$\displaystyle\leq\mathbb{E}[F(x^{k})-F(x^{*})]-\frac{\alpha_{l,k}}{2}\mathbb{E}[\lVert\nabla
F(x^{k})\rVert^{2}]$
$\displaystyle\quad\quad-\frac{\mathbb{E}[\lVert\bar{h}_{f}^{k}\rVert^{2}]}{2}(\alpha_{l,k}-L_{F}\alpha_{b,k}^{2}-4L_{y}^{2}\alpha_{b,k}^{2})+(\alpha_{b,k}L_{f}^{2}+2)\mathbb{E}[\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}]+$
$\displaystyle\quad\quad\alpha_{b,k}B^{2}+(2L_{y}^{2}\alpha_{b,k}^{2}+\frac{L_{F}\alpha_{b,k}^{2}}{2})\tilde{\sigma}_{f}^{2}$
$\displaystyle\stackrel{{\scriptstyle(e)}}{{\leq}}\mathbb{E}[F(x^{k})-F(x^{*})]-\frac{\alpha_{l,k}}{2}\mathbb{E}[\lVert\nabla
F(x^{k})\rVert^{2}]+(\alpha_{b,k}L_{f}^{2}+2)\mathbb{E}[\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}]+$
$\displaystyle\quad\quad\alpha_{b,k}B^{2}+(2L_{y}^{2}\alpha_{b,k}^{2}+\frac{L_{F}\alpha_{b,k}^{2}}{2})\tilde{\sigma}_{f}^{2}$
$\displaystyle=\mathbb{E}[V^{k}]-\frac{\alpha_{l,k}}{2}\mathbb{E}[\lVert\nabla
F(x^{k})\rVert^{2}]+(\alpha_{b,k}L_{f}^{2}+2)\mathbb{E}[\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}]+$
$\displaystyle\quad\quad\alpha_{b,k}B^{2}+(2L_{y}^{2}\alpha_{b,k}^{2}+\frac{L_{F}\alpha_{b,k}^{2}}{2})\tilde{\sigma}_{f}^{2}-\mathbb{E}[\lVert
y^{k}-y^{*}(x^{k})\rVert^{2}],$

where (e) is because
$\frac{1}{L_{F}+4L_{y}^{2}}\geq\frac{\alpha_{b,0}^{2}}{\alpha_{l,0}}$, which
guarantees that
$\alpha_{l,k}=\frac{\alpha_{l,0}}{\sqrt{k+1}}\geq(L_{F}+4L_{y}^{2})\alpha_{b,k}=(L_{F}+4L_{y}^{2})\frac{\alpha_{b,0}^{2}}{k+1}$.
∎

Lemma 5 and Lemma 6 give two alternatives for bounding the term $\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}$. Lemma 5 is based on the assumption
$\mathbb{E}_{\psi}[g(x,y^{*}(x);\psi)-g(x,y^{*}_{x,\psi};\psi)]\leq\sigma_{g}^{2}$,
$\forall x$. The proof for its one-variable analogous assumption is given in
Loizou et al. [29]. Here we follow a similar approach for the two-variable
function $g(x,y)$. Lemma 5 will be used in the proof of Theorem 3.

######  Lemma 5.

Suppose Assumptions 3, 5, and the bounded optimal function value assumption
(15) hold. Further assume that each sampled function $g(x,y;\psi)$ is convex,
then step size of the form 11 achieves the following:

$\displaystyle\mathbb{E}[\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}]\leq(1-\mu_{g}C_{k})^{T}\mathbb{E}[\lVert
y^{k}-y^{*}(x^{k})\rVert^{2}]+2\beta_{b,k}T\sigma_{g}^{2},$

where $C_{k}=\min\\{\frac{1}{2pL_{g}},\beta_{b,k}\\}$.

###### Proof.

We denote $h_{g}^{k,t}=\nabla_{y}g(x^{k},y^{k,t};\psi)$, and
$\mathcal{F}^{{}^{\prime}}_{k,t}$ be a filtration up to and including $x^{k}$
and $y^{k,t}$. We have,

$\displaystyle\lVert y^{k,t+1}-y^{*}(x^{k})\rVert^{2}$ $\displaystyle=\lVert
y^{k,t}-\beta_{k,t}h_{g}^{k,t}-y^{*}(x^{k})\rVert^{2}$ $\displaystyle=\lVert
y^{k,t}-y^{*}(x^{k})\rVert^{2}-2\beta_{k,t}\langle
y^{k,t}-y^{*}(x^{k}),h_{g}^{k,t}\rangle+\beta_{k,t}^{2}\lVert
h_{g}^{k,t}\rVert^{2}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\lVert
y^{k,t}-y^{*}(x^{k})\rVert^{2}-2\beta_{k,t}\langle
y^{k,t}-y^{*}(x^{k}),h_{g}^{k,t}\rangle+\frac{\beta_{k,t}}{p}[g(x^{k},y^{k,t};\psi)-g(x^{k},y^{*}_{x^{k},\psi};\psi)]$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}\lVert
y^{k,t}-y^{*}(x^{k})\rVert^{2}-2\beta_{k,t}\langle
y^{k,t}-y^{*}(x^{k}),h_{g}^{k,t}\rangle+2\beta_{k,t}[g(x^{k},y^{k,t};\psi)-g(x^{k},y^{*}_{x^{k},\psi};\psi)]$
$\displaystyle=\lVert y^{k,t}-y^{*}(x^{k})\rVert^{2}-2\beta_{k,t}\langle
y^{k,t}-y^{*}(x^{k}),h_{g}^{k,t}\rangle+$ $\displaystyle\quad\quad
2\beta_{k,t}[g(x^{k},y^{k,t};\psi)-g(x^{k},y^{*}(x^{k});\psi)+g(x^{k},y^{*}(x^{k});\psi)-g(x^{k},y^{*}_{x^{k},\psi};\psi)]$
$\displaystyle=\lVert y^{k,t}-y^{*}(x^{k})\rVert^{2}+2\beta_{k,t}[-\langle
y^{k,t}-y^{*}(x^{k}),h_{g}^{k,t}\rangle+g(x^{k},y^{k,t};\psi)-g(x^{k},y^{*}(x^{k});\psi)]+$
$\displaystyle\quad\quad
2\beta_{k,t}[g(x^{k},y^{*}(x^{k});\psi)-g(x^{k},y^{*}_{x^{k},\psi};\psi)]$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}\lVert
y^{k,t}-y^{*}(x^{k})\rVert^{2}+2C_{k}[-\langle
y^{k,t}-y^{*}(x^{k}),h_{g}^{k,t}\rangle+g(x^{k},y^{k,t};\psi)-g(x^{k},y^{*}(x^{k});\psi)]+$
$\displaystyle\quad\quad
2\beta_{b,k}[g(x^{k},y^{*}(x^{k});\psi)-g(x^{k},y^{*}_{x^{k},\psi};\psi)],$

where (a) is by Lemma 1, (b) is by choosing $p\geq\frac{1}{2}$, and, (c) is by
individual convexity of $g(x,y;\psi)$ such that $-\langle
y^{k,t}-y^{*}(x^{k}),h_{g}^{k,t}\rangle+g(x^{k},y^{k,t};\psi)-g(x^{k},y^{*}(x^{k});\psi)\leq
0$ and recalling that
$C_{k}=\min\\{\frac{1}{2pL_{g}},\beta_{b,k}\\}\leq\beta_{k,t}$ by Lemma 1.
Take expectation conditioned on $\mathcal{F}^{{}^{\prime}}_{k,t}$ and note
that
$\mathbb{E}[h_{g}^{k,t}|\mathcal{F}^{{}^{\prime}}_{k,t}]=\nabla_{y}g(x^{k},y^{k+1})$,
$\mathbb{E}[g(x^{k},y^{k,t};\psi)|\mathcal{F}^{{}^{\prime}}_{k,t}]=g(x^{k},y^{k,t})$,
and $\mathbb{E}[g(x^{k},y^{*}(x^{k});\psi)]=g(x^{k},y^{*}(x^{k}))$:

$\displaystyle\mathbb{E}[\lVert
y^{k,t+1}-y^{*}(x^{k})\rVert^{2}|\mathcal{F}^{{}^{\prime}}_{k,t}]$
$\displaystyle\leq\lVert y^{k,t}-y^{*}(x^{k})\rVert^{2}+2C_{k}[-\langle
y^{k,t}-y^{*}(x^{k}),\nabla_{y}g(x^{k},y^{k+1})\rangle+$
$\displaystyle\quad\quad
g(x^{k},y^{k,t})-g(x^{k},y^{*}(x^{k}))]+2\beta_{b,k}\sigma_{g}^{2}.$ (27)

Now, based on the strong convexity of $g$ w.r.t. $y$,

$\displaystyle-\langle
y^{k,t}-y^{*}(x^{k}),\nabla_{y}g(x^{k},y^{k+1})\rangle+g(x^{k},y^{k,t})-g(x^{k},y^{*}(x^{k}))\leq\frac{\mu_{g}}{2}\lVert
y^{k,t}-y^{*}(x^{k})\rVert^{2},$

we can further obtain (by taking total expectations of (27) and using strong-
convexity):

$\displaystyle\mathbb{E}[\lVert y^{k,t+1}-y^{*}(x^{k})\rVert^{2}]$
$\displaystyle\leq(1-\mu_{g}C_{k})\mathbb{E}[\lVert
y^{k,t}-y^{*}(x^{k})\rVert^{2}]+2\beta_{b,k}\sigma_{g}^{2}.$

Solve this recursively from $t=0$ to $t=T-1$ (recall $T$ is the total number
of lower-level steps, $y^{k,0}=y^{k}$ and $y^{k+1}=y^{k,T}$):

$\displaystyle\mathbb{E}[\lVert y^{k+1}-y^{*}(x^{k})\rVert^{2}]$
$\displaystyle\leq(1-\mu_{g}C_{k})^{T}\mathbb{E}[\lVert
y^{k}-y^{*}(x^{k})\rVert^{2}]+2\beta_{b,k}\sigma_{g}^{2}\sum_{j=0}^{T-1}(1-\mu_{g}C_{k})^{j}$
$\displaystyle\leq(1-\mu_{g}C_{k})^{T}\mathbb{E}[\lVert
y^{k}-y^{*}(x^{k})\rVert^{2}]+2\beta_{b,k}T\sigma_{g}^{2}.$

∎

Lemma 6 is based on the standard bounded variance assumption
$\mathbb{E}_{\psi}[\lVert\nabla_{y}g(x,y^{*}(x);\psi)-\nabla_{y}g(x,y^{*}(x))\rVert^{2}]\leq\sigma_{g}^{2}$,
$\forall x,y$ in the bi-level optimization literature [19, 4]. Lemma 6 will be
used in the proof of Theorem 4.

######  Lemma 6.

Suppose Assumptions 3, 5 and the bounded variance assumption (16) hold.
Suppose that
$p\geq\max\\{\frac{\mu_{g}}{\mu_{g}+L_{g}},\frac{\mu_{g}+L_{g}}{4L_{g}}\\}$,
$\beta_{b,0}\leq\min\\{\frac{2}{\mu_{g}+L_{g}},\frac{\mu_{g}+L_{g}}{2\mu_{g}L_{g}},\frac{1}{2pLg-\frac{2\mu_{g}L_{g}}{\mu_{g}+L_{g}}}\\}$.
Then step size of the form 11 achieves the following:

$\displaystyle\mathbb{E}[\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}]\leq(\frac{\beta_{b,k}}{C_{k}}-\frac{2\mu_{g}L_{g}}{\mu_{g}+L_{g}}\beta_{b,k})^{T}\mathbb{E}[\lVert
y^{k}-y^{*}(x^{k})\rVert^{2}]+T\beta_{b,k}^{2}\sigma_{g}^{2},$

where $C_{k}=\min\\{\frac{1}{2pL_{g}},\beta_{b,k}\\}$.

###### Proof.

Similar to the proof of Lemma 5, we can start with

$\displaystyle\lVert y^{k,t+1}-y^{*}(x^{k})\rVert^{2}=\lVert
y^{k,t}-y^{*}(x^{k})\rVert^{2}-2\beta_{k,t}\langle
y^{k,t}-y^{*}(x^{k}),h_{g}^{k,t}\rangle+\beta_{k,t}^{2}\lVert
h_{g}^{k,t}\rVert^{2}.$

Divide both sides by $\beta_{k,t}$

$\displaystyle\frac{\lVert
y^{k,t+1}-y^{*}(x^{k})\rVert^{2}}{\beta_{k,t}}=\frac{\lVert
y^{k,t}-y^{*}(x^{k})\rVert^{2}}{\beta_{k,t}}-2\langle
y^{k,t}-y^{*}(x^{k}),h_{g}^{k,t}\rangle+\beta_{k,t}\lVert
h_{g}^{k,t}\rVert^{2}.$

Then use the facts that $\beta_{k,t}\leq\beta_{b,k}$ and $\beta_{k,t}\geq
C_{k}$,

$\displaystyle\frac{\lVert
y^{k,t+1}-y^{*}(x^{k})\rVert^{2}}{\beta_{b,k}}\leq\frac{\lVert
y^{k,t}-y^{*}(x^{k})\rVert^{2}}{C_{k}}-2\langle
y^{k,t}-y^{*}(x^{k}),h_{g}^{k,t}\rangle+\beta_{b,k}\lVert
h_{g}^{k,t}\rVert^{2}\,.$

Next, take expectation conditioned on the Filtration
$\mathcal{F}^{{}^{\prime}}_{k,t}$ up to and including $x^{k}$ and $y^{k,t}$

$\displaystyle\frac{1}{\beta_{b,k}}\mathbb{E}[\lVert
y^{k,t+1}-y^{*}(x^{k})\rVert^{2}|\mathcal{F}^{{}^{\prime}}_{k,t}]$
$\displaystyle\leq\frac{1}{C_{k}}\mathbb{E}[\lVert
y^{k,t}-y^{*}(x^{k})\rVert^{2}|\mathcal{F}^{{}^{\prime}}_{k,t}]-2\langle
y^{k,t}-y^{*}(x^{k}),\nabla g(x^{k},y^{k,t})\rangle$
$\displaystyle\qquad\qquad+\beta_{b,k}\mathbb{E}[\lVert
h_{g}^{k,t}\rVert^{2}|\mathcal{F}^{{}^{\prime}}_{k,t}]$
$\displaystyle=\frac{1}{C_{k}}\mathbb{E}[\lVert
y^{k,t}-y^{*}(x^{k})\rVert^{2}|\mathcal{F}^{{}^{\prime}}_{k,t}]-2\langle
y^{k,t}-y^{*}(x^{k}),\nabla g(x^{k},y^{k,t})\rangle$
$\displaystyle\qquad\qquad+\beta_{b,k}\mathbb{E}[\lVert h_{g}^{k,t}-\nabla
g(x^{k},y^{k,t})+\nabla
g(x^{k},y^{k,t})\rVert^{2}|\mathcal{F}^{{}^{\prime}}_{k,t}]$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\frac{1}{C_{k}}\mathbb{E}[\lVert
y^{k,t}-y^{*}(x^{k})\rVert^{2}|\mathcal{F}^{{}^{\prime}}_{k,t}]-2\langle
y^{k,t}-y^{*}(x^{k}),\nabla g(x^{k},y^{k,t})\rangle$
$\displaystyle\qquad\qquad+\beta_{b,k}\sigma_{g}^{2}+\beta_{b,k}\lVert\nabla
g(x^{k},y^{k,t})\rVert^{2},$

where (a) is by the bounded variance assumption (16). Based on strong-
convexity of $g$ [34, Theorem 2.1.11], we have

$\displaystyle\frac{1}{\beta_{b,k}}\mathbb{E}[\lVert
y^{k,t+1}-y^{*}(x^{k})\rVert^{2}|\mathcal{F}^{{}^{\prime}}_{k,t}]$
$\displaystyle\leq\frac{1}{C_{k}}\lVert
y^{k,t}-y^{*}(x^{k})\rVert^{2}-\frac{2\mu_{g}L_{g}}{\mu_{g}+L_{g}}\lVert
y^{k,t}-y^{*}(x^{k})\rVert^{2}$
$\displaystyle\qquad\qquad-\frac{2}{\mu_{g}+L_{g}}\lVert\nabla
g(x^{k},y^{k,t})\rVert^{2}+\beta_{b,k}\sigma_{g}^{2}+\beta_{b,k}\lVert\nabla
g(x^{k},y^{k,t})\rVert^{2}$

Multiply by $\beta_{b,k}$ in both sides, group terms, and take total
expectations to reach

$\displaystyle\mathbb{E}[\lVert y^{k,t+1}-y^{*}(x^{k})\rVert^{2}]$
$\displaystyle\leq(\frac{\beta_{b,k}}{C_{k}}-\frac{2\mu_{g}L_{g}}{\mu_{g}+L_{g}}\beta_{b,k})\mathbb{E}[\lVert
y^{k,t}-y^{*}(x^{k})\rVert^{2}]+\beta_{b,k}^{2}\sigma_{g}^{2}$
$\displaystyle\qquad\qquad+\beta_{b,k}(\beta_{b,k}-\frac{2}{\mu_{g}+L_{g}})\lVert\nabla
g(x^{k},y^{k,t})\rVert^{2}$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}(\frac{\beta_{b,k}}{C_{k}}-\frac{2\mu_{g}L_{g}}{\mu_{g}+L_{g}}\beta_{b,k})\mathbb{E}[\lVert
y^{k,t}-y^{*}(x^{k})\rVert^{2}]+\beta_{b,k}^{2}\sigma_{g}^{2},$

where in (b) we have chosen
$\beta_{b,k}\leq\beta_{b,0}\leq\frac{2}{\mu_{g}+L_{g}},\forall k$. Solving
this recursion similar to Lemma 5, we obtain:

$\displaystyle\mathbb{E}[\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}]\leq(\frac{\beta_{b,k}}{C_{k}}-\frac{2\mu_{g}L_{g}}{\mu_{g}+L_{g}}\beta_{b,k})^{T}\mathbb{E}[\lVert
y^{k}-y^{*}(x^{k})\rVert^{2}]+T\beta_{b,k}^{2}\sigma_{g}^{2}.$ (28)

In (28), we require
$\frac{\beta_{b,k}}{C_{k}}-\frac{2\mu_{g}L_{g}}{\mu_{g}+L_{g}}\beta_{b,k}\geq
0$ (recall $C_{k}=\min\\{\frac{1}{2pL_{g}}$,
$\beta_{b,k}\\}\leq\beta_{k,t}\leq\beta_{b,k}$), this is equivalent to
$C_{k}\leq\frac{\mu_{g}+L_{g}}{2\mu_{g}L_{g}}$. In case of
$C_{k}=\frac{1}{2pL_{g}}$, we choose $p\geq\frac{\mu_{g}}{\mu_{g}+L_{g}}$ (to
avoid contradictions, we also choose $p\geq\frac{\mu_{g}+L_{g}}{4L_{g}}$). In
case of $C_{k}=\beta_{b,k}$, we choose
$\beta_{b,k}\leq\beta_{b,0}\leq\frac{\mu_{g}+L_{g}}{2\mu_{g}L_{g}}$.
We also require
$\frac{\beta_{b,k}}{C_{k}}-\frac{2\mu_{g}L_{g}}{\mu_{g}+L_{g}}\beta_{b,k}\leq
1$. In case of $C_{k}=\beta_{b,k}$, this is equivalent to
$\frac{2\mu_{g}L_{g}}{\mu_{g}+L_{g}}\beta_{b,k}\geq 0$, which is satisfied by
all $\beta_{b,k}$. In case of $C_{k}=\frac{1}{2pL_{g}}$, we choose
$\beta_{b,k}\leq\beta_{b,0}\leq\frac{1}{2pL_{g}-\frac{2\mu_{g}L_{g}}{\mu_{g}+L_{g}}}$.
Puting everything together, we have
$p\geq\max\\{\frac{\mu_{g}}{\mu_{g}+L_{g}},\frac{\mu_{g}+L_{g}}{4L_{g}}\\}$,
and
$\beta_{b,0}\leq\min\\{\frac{2}{\mu_{g}+L_{g}},\frac{\mu_{g}+L_{g}}{2\mu_{g}L_{g}},\frac{1}{2pLg-\frac{2\mu_{g}L_{g}}{\mu_{g}+L_{g}}}\\}$.
∎

###  A.2 Single-level Convex Proofs

####  A.2.1 Proof of Theorem 1

The proof of Theorem 1 is similar to the standard proof of decaying-step SGD
(GD) that can be found in e.g. Beck [1]. Here, we give the proof for
completeness.

$\displaystyle\lVert x^{k+1}-x^{*}\rVert^{2}$ $\displaystyle=\lVert
x^{k}-x^{*}\rVert^{2}-2\gamma_{k}\langle x^{k}-x^{*},\nabla
f_{i}(x^{k})\rangle+\gamma_{k}^{2}\lVert\nabla f_{i}(x^{k})\rVert^{2}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\lVert
x^{k}-x^{*}\rVert^{2}-2\gamma_{k}(f_{i}(x^{k})-f_{i}(x^{*}))+\gamma_{k}^{2}\lVert\nabla
f_{i}(x^{k})\rVert^{2}$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}\lVert
x^{k}-x^{*}\rVert^{2}-2\gamma_{k}(f_{i}(x^{k})-f_{i}(x^{*}))+\gamma_{b,k}^{2}G^{2},$

where (a) is by individual convexity of $f_{i}$, and (b) is by Assumption 2.
Take conditional expectation and assume that $\gamma_{k}$ is independent of
sample $k$

$\displaystyle\mathbb{E}[\lVert x^{k+1}-x^{*}\rVert^{2}|x^{k}]$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}\mathbb{\lVert}x^{k}-x^{*}\rVert^{2}-2E[\gamma_{k}|x^{k}](f(x^{k})-f(x^{*}))+\gamma_{b,k}^{2}G^{2}$
$\displaystyle\stackrel{{\scriptstyle(d)}}{{\leq}}\mathbb{\lVert}x^{k}-x^{*}\rVert^{2}-2\gamma_{l,k}(f(x^{k})-f(x^{*}))+\gamma_{b,k}^{2}G^{2},$

where (c) is by independence of $\gamma_{k}$ and $x^{k}$, and (d) is because
$\gamma_{l,k}\leq\gamma_{k},\forall k$. Take total expectation and rearrange

$\displaystyle 2\gamma_{l,k}\mathbb{E}[f(x^{k})-f(x^{*})]\leq\mathbb{E}[\lVert
x^{k}-x^{*}\rVert^{2}]-E[\lVert
x^{k+1}-x^{*}\rVert^{2}]+\gamma_{b,k}^{2}G^{2}$

Using the fact that $\gamma_{l,K-1}\leq\gamma_{l,k}\quad\forall k\in[K-1]$ and
set $\gamma_{b,k}=\frac{\gamma_{0}}{\sqrt{k+1}}$, we obtain

$\displaystyle
2\gamma_{l,K-1}\mathbb{E}[f(x^{k})-f(x^{*})]\leq\mathbb{E}[\lVert
x^{k}-x^{*}\rVert^{2}]-E[\lVert
x^{k+1}-x^{*}\rVert^{2}]+\frac{\gamma_{0}^{2}}{k+1}G^{2}$

Summing over $k=0$ to $k=K-1$ we obtain

$\displaystyle
2\gamma_{l,K-1}\frac{1}{K}\sum_{k=0}^{K-1}\mathbb{E}[f(x^{k})-f(x^{*})]$
$\displaystyle\leq\frac{\lVert x^{0}-x^{*}\rVert^{2}-\mathbb{E}[\lVert
x^{K}-x^{*}\rVert^{2}]}{K}+\gamma_{0}^{2}G^{2}\frac{1}{K}\sum_{k=0}^{K-1}\frac{1}{k+1}$
$\displaystyle\leq\frac{\lVert
x^{0}-x^{*}\rVert^{2}}{K}+\frac{\gamma_{0}^{2}G^{2}\log(K)}{K}$

Define $\bar{x}=\frac{1}{K}\sum_{k=0}^{K-1}x^{k}$, apply Jensen’s inequality
and rearrange

$\displaystyle\mathbb{E}[f(\bar{x}^{K})-f(x^{*})]\leq\frac{\lVert
x^{0}-x^{*}\rVert^{2}}{2\gamma_{l,K-1}K}+\frac{\gamma_{0}^{2}G^{2}\log(K)}{2\gamma_{l,K-1}K}.$

###  A.3 Proof of Theorem 2

The approach of Theorem 2 is similar to [16, Theorem 3.2]. The crucial
difference is that the step size in [16, Theorem 3.2] is constant
($\gamma_{k}=\frac{1}{L_{\max}}$) for $k\leq 4\lceil\kappa\rceil$, whereas for
envelope-type step size it is of the form:

$\displaystyle\gamma_{k}=\min\\{\max\\{\gamma_{l,k},\tilde{\gamma}_{k}\\},\gamma_{b,k}\\},\quad\gamma_{l,k}=\min\\{w,\gamma_{b,k}\\},$

where $\bar{\gamma}_{k}$ can be (e.g. in the case of SPSB ):

$\displaystyle\bar{\gamma}_{k}=\min\\{\frac{f_{i_{k}}(x^{k})-l_{i_{k}}^{*}}{c_{k}\lVert\nabla
f_{i_{k}}(x^{k})\rVert^{2}},\frac{\gamma_{b,0}}{k+1}\\},\quad
w=\frac{1}{2c_{k}L_{\max}}\,.$

The proof of Theorem 2 suggests that the step size can be either
$\frac{f_{i_{k}}(x^{k})-l_{i_{k}}^{*}}{c_{k}\lVert\nabla
f_{i_{k}}(x^{k})\rVert^{2}}$ or $\frac{\gamma_{b,0}}{k+1}$ depending on their
magnitudes for $k\leq k_{0}-1$
($k_{0}=\max\\{1,\lceil\gamma_{0}/\omega\rceil-1\\}$). After $k_{0}$
iterations, the step size is $\gamma_{k}=\frac{\gamma_{b,0}}{k+1}$. This
finding is numerically confirmed by the experimental results in Section B.

To proceed with the proof, we have:

$\displaystyle\lVert x^{k+1}-x^{*}\rVert^{2}$ $\displaystyle=\lVert
x^{k}-\gamma_{k}\nabla f_{i}(x^{k})-x^{*}\rVert^{2}$ $\displaystyle=\lVert
x^{k}-x^{*}\rVert^{2}-2\gamma_{k}\langle\nabla
f_{i}(x^{k}),x^{k}-x^{*}\rangle+\gamma_{k}^{2}\lVert\nabla
f_{i}(x^{k})\rVert^{2}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\lVert
x^{k}-x^{*}\rVert^{2}-2\gamma_{k}\langle\nabla
f_{i}(x^{k}),x^{k}-x^{*}\rangle+\gamma_{b,k}^{2}\lVert\nabla
f_{i}(x^{k})\rVert^{2},$

where (a) is because $\gamma_{k}\leq\gamma_{b,k},\forall k$. Take conditional
expectations

$\displaystyle\mathbb{E}[\lVert x^{k+1}-x^{*}\rVert^{2}|x^{k}]$
$\displaystyle\leq\lVert
x^{k}-x^{*}\rVert^{2}-2\mathbb{E}[\gamma_{k}|x^{k}]\langle\nabla
f(x^{k}),x^{k}-x^{*}\rangle+\gamma_{b,k}^{2}\mathbb{E}[\lVert\nabla
f_{i}(x^{k})\rVert^{2}|x^{k}]$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}\lVert
x^{k}-x^{*}\rVert^{2}-\mu\mathbb{E}[\gamma_{k}|x^{k}]\lVert
x^{k}-x^{*}\rVert^{2}-2\mathbb{E}[\gamma_{k}|x^{k}][f(x^{k})-f(x^{*})]+\gamma_{b,k}^{2}G^{2}$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}\lVert
x^{k}-x^{*}\rVert^{2}-\mu\gamma_{l,k}\lVert
x^{k}-x^{*}\rVert^{2}-2\gamma_{l,k}[f(x^{k})-f(x^{*})]+\gamma_{b,k}^{2}G^{2}$

where (b) is by bounded gradients assumption and strong convexity of $f$. Take
total expectation and rearrange

$\displaystyle
2\mathbb{E}[f(x^{k})-f(x^{*})]\leq\frac{(1-\mu\gamma_{l,k})\mathbb{E}[\lVert
x^{k}-x^{*}\rVert^{2}]-\mathbb{E}[\lVert
x^{k+1}-x^{*}\rVert^{2}]}{\gamma_{l,k}}+\frac{\gamma_{b,k}^{2}G^{2}}{\gamma_{l,k}}$

Choose $k_{0}=\max\\{1,\lceil\gamma_{0}/\omega\rceil-1\\}$, then for $\forall
k$ s.t. $k\geq k_{0}$, we have
$\gamma_{l,k}=\min\\{\omega,\gamma_{b,k}\\}=\gamma_{b,k}=\frac{\gamma_{0}}{k+1}$,
which means $\gamma_{k}=\frac{\gamma_{0}}{k+1}$ after $k_{0}$ steps. Within
the first $k_{0}$ steps, the step size is
$\gamma_{k}=\min\\{\omega,\gamma_{b,k}\\}$. Hence, for $k\geq k_{0}$ we have

$\displaystyle
2\mathbb{E}[f(x^{k})-f(x^{*})]\leq\frac{(1-\mu\gamma_{l,k})\mathbb{E}[\lVert
x^{k}-x^{*}\rVert^{2}]-\mathbb{E}[\lVert
x^{k+1}-x^{*}\rVert^{2}]}{\gamma_{l,k}}+\frac{\gamma_{0}G^{2}}{k+1}$

Now, sum from $k=k_{0}$ to $K-1$

$\displaystyle
2\sum_{k=k_{0}}^{K-1}\mathbb{E}[f(x^{k})-f(x^{*})]\leq\sum_{k=k_{0}}^{K-1}\frac{(1-\mu\gamma_{l,k})\mathbb{E}[\lVert
x^{k}-x^{*}\rVert^{2}]-\mathbb{E}[\lVert
x^{k+1}-x^{*}\rVert^{2}]}{\gamma_{l,k}}+\sum_{k=k_{0}}^{K-1}\frac{\gamma_{0}G^{2}}{k+1}$
(29)

For the first term in (29), call it $A$, we expand it as

$\displaystyle A$ $\displaystyle=\sum_{k=k_{0}+1}^{K-1}\mathbb{E}[\lVert
x^{k}-x^{*}\rVert^{2}](\frac{1}{\gamma_{l,k}}-\frac{1}{\gamma_{l,k-1}}-\mu)+\mathbb{E}[\lVert
x_{k_{0}}-x^{*}\rVert^{2}](\frac{1}{\gamma_{l,k_{0}}}-\mu)-\frac{\mathbb{E}[\lVert
x_{K}-x^{*}\rVert^{2}]}{\gamma_{l,K-1}}$
$\displaystyle\leq\sum_{k=k_{0}+1}^{K-1}\mathbb{E}[\lVert
x^{k}-x^{*}\rVert^{2}](\frac{k+1}{\gamma_{0}}-\frac{k}{\gamma_{0}}-\mu)+\mathbb{E}[\lVert
x_{k_{0}}-x^{*}\rVert^{2}](\frac{k_{0}+1}{\gamma_{0}}-\mu)$
$\displaystyle\stackrel{{\scriptstyle(d)}}{{\leq}}\mathbb{E}[\lVert
x_{k_{0}}-x^{*}\rVert^{2}]\mu k_{0}.$

where (d) is because $\gamma_{0}\geq\frac{1}{\mu}$, we have
$\frac{k+1}{\gamma_{0}}-\frac{k}{\gamma_{0}}-\mu\leq 0,\forall k$ and
$\frac{k_{0}+1}{\gamma_{0}}-\mu\leq\mu k_{0}$. For the second term in (29),
call it B, we have

$\displaystyle
B=\sum_{k=k_{0}}^{K-1}\frac{\gamma_{0}G^{2}}{k+1}\leq\gamma_{0}G^{2}\int_{k_{0}}^{K-1}\frac{1}{x+1}dx=\gamma_{0}G^{2}[\log(K)-\log(k_{0}+1)]\leq\gamma_{0}G^{2}\log
K$

Putting $A$ and $B$ together, we obtain

$\displaystyle\frac{1}{K-k_{0}}\sum_{k=k_{0}}^{K-1}\mathbb{E}[f(x^{k})-f(x^{*})]\leq\frac{\mathbb{E}[\lVert
x^{k_{0}}-x^{*}\rVert^{2}]\mu
k_{0}}{2(K-k_{0})}+\frac{\gamma_{0}G^{2}\log(K)}{2(K-k_{0})}$ (30)

Within the first $k_{0}-1$ iterations, similarly to the above, we have

$\displaystyle\mathbb{E}[\lVert x^{k+1}-x^{*}\rVert]$
$\displaystyle\leq(1-\mu\gamma_{l,k})\mathbb{E}[\lVert
x^{k}-x^{*}\rVert^{2}]-2\gamma_{l,k}\mathbb{E}[f(x^{k})-f(x^{*})]+\gamma_{b,k}^{2}G^{2}$
$\displaystyle\leq(1-\mu\gamma_{l,k})\mathbb{E}[\lVert
x^{k}-x^{*}\rVert^{2}]+\gamma_{b,k}^{2}G^{2}\,.$

For the first $k\leq k_{0}-1$ iterations, $\gamma_{l,k}=\omega$ where
$\omega\mu<1$; thus, we obtain the following

$\displaystyle\mathbb{E}[\lVert x^{k+1}-x^{*}\rVert]\leq(1-\mu\omega)E[\lVert
x^{k}-x^{*}\rVert^{2}]+\gamma_{b,k}^{2}G^{2}.$

Solve this recursively,

$\displaystyle\mathbb{E}[\lVert x^{k_{0}}-x^{*}\rVert^{2}]$
$\displaystyle\leq(1-\mu\omega)^{k_{0}}\lVert
x^{0}-x^{*}\rVert^{2}+\sum_{k=0}^{k_{0}-1}(1-\mu\omega)^{k_{0}-k-1}\frac{\gamma_{0}^{2}G^{2}}{(k+1)^{2}}$
$\displaystyle\leq(1-\mu\omega)^{k_{0}}\lVert
x^{0}-x^{*}\rVert^{2}+\sum_{k=0}^{k_{0}-1}\frac{\gamma_{0}^{2}G^{2}}{(k+1)^{2}}$
$\displaystyle\leq(1-\mu\omega)^{k_{0}}\lVert
x^{0}-x^{*}\rVert^{2}+\gamma_{0}^{2}G^{2}\int_{0}^{k_{0}-1}\frac{1}{(1+x)^{2}}dx$
$\displaystyle=(1-\mu\omega)^{k_{0}}\lVert
x^{0}-x^{*}\rVert^{2}+\gamma_{0}^{2}G^{2}(1-\frac{1}{k_{0}})$
$\displaystyle\leq(1-\mu\omega)^{k_{0}}\lVert
x^{0}-x^{*}\rVert^{2}+\gamma_{0}^{2}G^{2}.$ (31)

Putting this into (30), we obtain

$\displaystyle\frac{1}{K-k_{0}}\sum_{k=k_{0}}^{K-1}\mathbb{E}[f(x^{k})-f(x^{*})]\leq\frac{\mu
k_{0}}{2(K-k_{0})}\\{\exp(-k_{0}\mu\omega)\lVert
x_{0}-x^{*}\rVert^{2}+\gamma_{0}^{2}G^{2}\\}+\frac{\gamma_{0}G^{2}\log
K}{2(K-k_{0})}.$

Define $\bar{x}_{K}=\frac{1}{K-k_{0}}\sum_{k=k_{0}}^{K-1}x^{k}$, then by
Jensen’s inequality we have

$\displaystyle\mathbb{E}[f(\bar{x}_{K})-f(x^{*})]\leq\frac{\mu
k_{0}}{2(K-k_{0})}\\{\exp(-k_{0}\mu\omega)\lVert
x_{0}-x^{*}\rVert^{2}+\gamma_{0}^{2}G^{2}\\}+\frac{\gamma_{0}G^{2}\log
K}{2(K-k_{0})},$

where $k_{0}=\max\\{1,\lceil\gamma_{0}/\omega\rceil-1\\}$ and
$\gamma_{0}\geq\frac{1}{\mu}$.

###  A.4 Bi-level Proofs

####  A.4.1 Proof of Theorem 3

Start with Lemma 4:

$\displaystyle E[V^{k+1}]$
$\displaystyle\leq\mathbb{E}[V^{k}]-\frac{\alpha_{l,k}}{2}\mathbb{E}[\lVert\nabla
F(x^{k})\rVert^{2}]+(\alpha_{b,k}L_{f}^{2}+2)\mathbb{E}[\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}]+$
$\displaystyle\quad\quad\alpha_{b,k}B^{2}+(2L_{y}^{2}\alpha_{b,k}^{2}+\frac{L_{F}\alpha_{b,k}^{2}}{2})\tilde{\sigma}_{f}^{2}-\mathbb{E}[\lVert
y^{k}-y^{*}(x^{k})\rVert^{2}].$

We substitute the result of Lemma 5 for the expression $\mathbb{E}[\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}]$,

$\displaystyle E[V^{k+1}]$
$\displaystyle\leq\mathbb{E}[V^{k}]-\frac{\alpha_{l,k}}{2}\mathbb{E}[\lVert\nabla
F(x^{k})\rVert^{2}]+[(\alpha_{b,k}L_{f}^{2}+2)(1-\mu_{g}C_{k})^{T}-1]\mathbb{E}[\lVert
y^{k}-y^{*}(x^{k})\rVert^{2}]+$ $\displaystyle\quad\quad
2\alpha_{b,k}\beta_{b,k}TL_{f}^{2}\sigma_{g}^{2}+4\beta_{b,k}T\sigma_{g}^{2}+\alpha_{b,k}B^{2}+(2L_{y}^{2}+\frac{L_{F}}{2})\alpha_{b,k}^{2}\tilde{\sigma}_{f}^{2}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\mathbb{E}[V^{k}]-\frac{\alpha_{l,k}}{2}\mathbb{E}[\lVert\nabla
F(x^{k})\rVert^{2}]+[(\alpha_{b,0}L_{f}^{2}+2)(1-\mu_{g}C_{K-1})^{T}-1]\mathbb{E}[\lVert
y^{k}-y^{*}(x^{k})\rVert^{2}]+$ $\displaystyle\quad\quad
2\alpha_{b,k}\beta_{b,k}TL_{f}^{2}\sigma_{g}^{2}+4\beta_{b,k}T\sigma_{g}^{2}+\alpha_{b,k}B^{2}+(2L_{y}^{2}+\frac{L_{F}}{2})\alpha_{b,k}^{2}\tilde{\sigma}_{f}^{2}$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}E[V^{k}]-\frac{\alpha_{l,k}}{2}E[\lVert\nabla
F(x^{k})\rVert^{2}]+2\alpha_{b,k}\beta_{b,k}TL_{f}^{2}\sigma_{g}^{2}+4\beta_{b,k}T\sigma_{g}^{2}+\alpha_{b,k}B^{2}+(2L_{y}^{2}+\frac{L_{F}}{2})\alpha_{b,k}^{2}\tilde{\sigma}_{f}^{2},$

where (a) is by $\alpha_{b,0}\geq\alpha_{b,k},C_{K-1}\leq C_{k},\forall k$,
hence
$(\alpha_{b,0}L_{f}^{2})(1-\mu_{g}C_{K-1})^{T}\geq(\alpha_{b,k}L_{f}^{2})(1-\mu_{g}C_{k})^{T}$
(recall that $C_{k}=\min\\{\frac{1}{2pL_{g}},\beta_{b,k}\\}$ in Lemma 1); (b)
is by $T\geq\frac{\log[\alpha_{b,0}L_{f}^{2}+2]}{-\log(1-\mu C_{K-1})}$, which
implies that $(\alpha_{b,0}L_{f}^{2}+2)(1-\mu_{g}C_{K-1})^{T}\leq 1$. Now,
rearrange and use the fact that $\alpha_{l,K-1}\leq\alpha_{l,k}$,

$\displaystyle\alpha_{l,K-1}\mathbb{E}[\lVert\nabla F(x^{k})\rVert^{2}]\leq
2\mathbb{E}[V^{k}]-2\mathbb{E}[V^{k+1}]+4\alpha_{b,k}\beta_{b,k}TL_{f}^{2}\sigma_{g}^{2}+8\beta_{b,k}T\sigma_{g}^{2}+2\alpha_{b,k}B^{2}+(4L_{y}^{2}+L_{F})\alpha_{b,k}^{2}\tilde{\sigma}_{f}^{2}.$

Then sum over $k=0$ to $K-1$:

$\displaystyle\frac{1}{K}\sum_{k=0}^{K-1}\mathbb{E}[\lVert\nabla
F(x^{k})\rVert^{2}]$
$\displaystyle\leq\frac{2V^{0}}{\alpha_{l,K-1}K}+\frac{4TL_{f}^{2}\sigma_{g}^{2}}{\alpha_{l,K-1}K}\sum_{k=0}^{K-1}\alpha_{b,k}\beta_{b,k}+\frac{8T\sigma_{g}^{2}}{\alpha_{l,K-1}K}\sum_{k=0}^{K-1}\beta_{b,k}+\frac{2B^{2}}{\alpha_{l,K-1}K}\sum_{k=0}^{K-1}\alpha_{b,k}+$
$\displaystyle\quad\quad\frac{(4L_{y}^{2}+L_{F})\tilde{\sigma}_{f}^{2}}{\alpha_{l,K-1}K}\sum_{k=0}^{K-1}\alpha_{b,k}^{2}$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{=}}\frac{2V^{0}}{\alpha_{l,K-1}K}+\frac{4TL_{f}^{2}\sigma_{g}^{2}\alpha_{b,0}\beta_{b,0}}{\alpha_{l,K-1}K}\sum_{k=0}^{K-1}\frac{1}{(k+1)^{1.5}}+\frac{8T\sigma_{g}^{2}\beta_{b,0}}{\alpha_{l,K-1}K}\sum_{k=0}^{K-1}\frac{1}{k+1}+$
$\displaystyle\quad\quad+\frac{2B^{2}\alpha_{b,0}}{\alpha_{l,K-1}K}\sum_{k=0}^{K-1}\frac{1}{\sqrt{k+1}}+\frac{(4L_{y}^{2}+L_{F})\tilde{\sigma}_{f}^{2}\alpha_{b,0}^{2}}{\alpha_{l,K-1}K}\sum_{k=0}^{K-1}\frac{1}{k+1}$
$\displaystyle\stackrel{{\scriptstyle(d)}}{{\leq}}\frac{2V^{0}}{\alpha_{l,0}\sqrt{K}}+\frac{8TL_{f}^{2}\sigma_{g}^{2}\alpha_{b,0}\beta_{b,0}}{\alpha_{l,0}\sqrt{K}}+\frac{8T\sigma_{g}^{2}\beta_{b,0}\log(K)}{\alpha_{l,0}\sqrt{K}}+$
$\displaystyle\quad\quad\frac{2B^{2}\alpha_{b,0}}{\alpha_{l,0}}+\frac{(4L_{y}^{2}+L_{F})\tilde{\sigma}_{f}^{2}\alpha_{b,0}^{2}\log(K)}{\alpha_{l,0}\sqrt{K}},$

where we substituted $\alpha_{b,k}=\frac{\alpha_{b,0}}{\sqrt{k+1}}$ and
$\beta_{b,k}=\frac{\beta_{b,0}}{k+1}$ in (c); (d) is based on
$\sum_{k=0}^{K-1}\frac{1}{(k+1)^{1.5}}\leq 2$,
$\sum_{k=0}^{K-1}\frac{1}{k+1}\leq\log(K)$, and
$\sum_{k=0}^{K-1}\frac{1}{\sqrt{k+1}}\leq\sqrt{K}$. Recall that in Lemma 2, we
have $L_{f}\sim\mathcal{O}(\kappa^{2})$, $L_{y}\sim\mathcal{O}(\kappa)$, and
$L_{F}\sim\mathcal{O}(\kappa^{3})$. Also recall that
$\alpha_{b,k}=\frac{\alpha_{b,0}}{\sqrt{k+1}}$ and
$\alpha_{l,k}=\frac{\alpha_{l,0}}{\sqrt{k+1}}$, hence we choose
$\alpha_{l,0}\sim\alpha_{b,0}\sim\mathcal{O}(\kappa^{-3})$,
$T\sim\mathcal{O}(\kappa)$, and $\beta_{b,0}\sim\mathcal{O}(\kappa^{-2})$.
Then we can obtain $\frac{1}{K}\sum_{k=0}^{K-1}\mathbb{E}[\lVert\nabla
F(x^{k})\rVert^{2}]\leq\mathcal{O}(\frac{\kappa^{3}}{\sqrt{K}}+\frac{\kappa^{2}\log(K)}{\sqrt{K}})$.

####  A.4.2 Theorem 4 and its proof

######  Theorem 4.

Suppose $f$ and $g$ satisfy Assumptions 3, 4, and 5, and, learning-rate upper
bounds $\alpha_{b,k}=\frac{\alpha_{b,0}}{\sqrt{k+1}}$ and
$\alpha_{l,k}=\frac{\alpha_{l,0}}{\sqrt{k+1}}$ with $\alpha_{b,0}$ and
$\alpha_{l,0}$ satisfying
$\frac{1}{L_{F}+4L_{y}^{2}}\geq\frac{\alpha_{b,0}^{2}}{\alpha_{l,0}}$ and
$\alpha_{l,0}\leq\alpha_{b,0}$. Further assume that $\alpha_{k}$ is
independent of the stochastic hypergradient $h_{f}^{k}$. Then, under the
Bounded-variance assumption in (16) with
$p\geq\max\\{\frac{\mu_{g}}{\mu_{g}+L_{g}},\frac{\mu_{g}+L_{g}}{4L_{g}}\\}$,
$\beta_{b,0}\leq\min\\{\frac{2}{\mu_{g}+L_{g}},\frac{\mu_{g}+L_{g}}{2\mu_{g}L_{g}},\frac{1}{2pLg-\frac{2\mu_{g}L_{g}}{\mu_{g}+L_{g}}}\\}$,
and
$T\geq\frac{\log(\frac{3}{2}\alpha_{b,0}L_{f}^{2}+2)}{\min\\{-\log(1-\frac{2\mu_{g}L_{g}}{\mu_{g}+L_{g}}\beta_{b,K-1}),-\log((2pL_{g}-\frac{2\mu_{g}L_{g}}{\mu_{g}+L_{g}})\beta_{b,0})\\}}$,
BiSPS achieves the following rate:

$\displaystyle\frac{1}{K}\sum_{k=0}^{K-1}\mathbb{E}[\lVert\nabla
F(x^{k})\rVert^{2}]$
$\displaystyle\leq\mathcal{\tilde{O}}(\frac{\kappa^{3}}{\sqrt{K}}+\frac{\kappa^{2}\log
K}{\sqrt{K}}).$

###### Proof.

Start with Lemma 4:

$\displaystyle E[V^{k+1}]$
$\displaystyle\leq\mathbb{E}[V^{k}]-\frac{\alpha_{l,k}}{2}\mathbb{E}[\lVert\nabla
F(x^{k})\rVert^{2}]+(\alpha_{b,k}L_{f}^{2}+2)\mathbb{E}[\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}]+$
$\displaystyle\quad\quad\alpha_{b,k}B^{2}+(2L_{y}^{2}\alpha_{b,k}^{2}+\frac{L_{F}\alpha_{b,k}^{2}}{2})\tilde{\sigma}_{f}^{2}-\mathbb{E}[\lVert
y^{k}-y^{*}(x^{k})\rVert^{2}].$

We substitute the result of Lemma 6 for the expression $\mathbb{E}[\lVert
y^{k+1}-y^{*}(x^{k})\rVert^{2}]$,

$\displaystyle E[V^{k+1}]$
$\displaystyle\leq\mathbb{E}[V^{k}]-\frac{\alpha_{l,k}}{2}\mathbb{E}[\lVert\nabla
F(x^{k})\rVert^{2}]+[(\alpha_{b,k}L_{f}^{2}+2)(\frac{\gamma_{b,k}}{C_{k}}-\frac{2\mu_{g}L_{g}}{\mu_{g}+L_{g}}\beta_{b,k})^{T}-1]\mathbb{E}[\lVert
y^{k}-y^{*}(x^{k})\rVert^{2}]+$
$\displaystyle\quad\quad+L_{f}^{2}T\sigma_{g}^{2}\alpha_{b,k}\beta_{b,k}^{2}+2T\sigma_{g}^{2}\beta_{b,k}^{2}+\alpha_{b,k}B^{2}+[2L_{y}^{2}+\frac{L_{F}}{2}]\alpha_{b,k}^{2}\tilde{\sigma}_{f}^{2}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}E[V^{k}]-\frac{\alpha_{l,k}}{2}\mathbb{E}[\lVert\nabla
F(x^{k})\rVert^{2}]+L_{f}^{2}T\sigma_{g}^{2}\alpha_{b,k}\beta_{b,k}^{2}+2T\sigma_{g}^{2}\beta_{b,k}^{2}+\alpha_{b,k}B^{2}+[2L_{y}^{2}+\frac{L_{F}}{2}]\alpha_{b,k}^{2}\tilde{\sigma}_{f}^{2},$

where in (a) we have chosen
$T\geq\max\\{\frac{\log(\alpha_{b,0}L_{f}^{2}+2)}{-\log(1-\frac{2\mu_{g}L_{g}}{\mu_{g}+L_{g}}\beta_{b,K-1})},\frac{\log(\alpha_{b,0}^{2}L_{f}^{2}+2)}{-\log(2pL_{g}-\frac{2\mu_{g}L_{g}}{\mu_{g}+L_{g}})\beta_{b,0}}\\}$.
This ensures that
$T\geq\frac{\log(\alpha_{b,0}L_{f}^{2}+2)}{-\log(\frac{\beta_{b,k}}{C_{k}}-\frac{2\mu_{g}L_{g}}{\mu_{g}+L_{g}}\beta_{b,k})},\forall
k$, which guarantees that
$(\alpha_{b,k}L_{f}^{2}+2)(\frac{\beta_{b,k}}{C_{k}}-\frac{2\mu_{g}L_{g}}{\mu_{g}+L_{g}}\beta_{b,k})^{T}-1\leq
0$. Now, rearrange and use the fact that $\alpha_{l,K-1}\leq\alpha_{l,k}$,

$\displaystyle\alpha_{l,K-1}\mathbb{E}[\lVert\nabla F(x^{k})\rVert^{2}]\leq
2\mathbb{E}[V^{k}]-2\mathbb{E}[V^{k+1}]+2\alpha_{b,k}\beta_{b,k}^{2}TL_{f}^{2}\sigma_{g}^{2}+4\beta_{b,k}^{2}T\sigma_{g}^{2}+2\alpha_{b,k}B^{2}+(4L_{y}^{2}+L_{F})\alpha_{b,k}^{2}\tilde{\sigma}_{f}^{2}.$

Sum over $k=0$ to $k=K-1$:

$\displaystyle\frac{1}{K}\sum_{k=0}^{K-1}\mathbb{E}[\lVert\nabla
F(x^{k})\rVert^{2}]$
$\displaystyle\leq\frac{2V^{0}}{\alpha_{l,K-1}K}+\frac{2TL_{f}^{2}\sigma_{g}^{2}}{\alpha_{l,K-1}K}\sum_{k=0}^{K-1}\alpha_{b,k}\beta_{b,k}^{2}+\frac{4T\sigma_{g}^{2}}{\alpha_{l,K-1}K}\sum_{k=0}^{K-1}\beta_{b,k}^{2}+$
$\displaystyle\quad\quad\frac{2B^{2}}{\alpha_{l,K-1}K}\sum_{k=0}^{K-1}\alpha_{b,k}+\frac{(4L_{y}^{2}+L_{F})\tilde{\sigma}_{f}^{2}}{\alpha_{l,K-1}K}\sum_{k=0}^{K-1}\alpha_{b,k}^{2}$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}\frac{2V^{0}}{\alpha_{l,K-1}K}+\frac{2TL_{f}^{2}\sigma_{g}^{2}\alpha_{b,0}\beta_{b,0}^{2}}{\alpha_{l,K-1}K}\sum_{k=0}^{K-1}\frac{1}{(k+1)^{3/2}}+\frac{4T\sigma_{g}^{2}\beta_{b,0}^{2}}{\alpha_{l,K-1}K}\sum_{k=0}^{K-1}\frac{1}{k+1}+$
$\displaystyle\quad\quad\frac{2B^{2}\alpha_{b,0}}{\alpha_{l,K-1}K}\sum_{k=0}^{K-1}\frac{1}{(k+1)^{0.5}}+\frac{(4L_{y}^{2}+L_{F})\tilde{\sigma}_{f}^{2}\alpha_{b,0}^{2}}{\alpha_{l,K-1}K}\sum_{k=0}^{K-1}\frac{1}{k+1}$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}\frac{2V^{0}}{\alpha_{l,K-1}K}+\frac{2TL_{f}^{2}\sigma_{g}^{2}\alpha_{b,0}\beta_{b,0}^{2}}{\alpha_{b,0}\sqrt{K}}+\frac{4T\sigma_{g}^{2}\beta_{b,0}^{2}\log(K)}{\alpha_{b,0}\sqrt{K}}+$
$\displaystyle\quad\quad+\frac{2B^{2}\alpha_{b,0}}{\alpha_{l,0}}+\frac{(4L_{y}^{2}+L_{F})\tilde{\sigma}_{f}^{2}\alpha_{b,0}^{2}\log(K)}{\alpha_{l,0}\sqrt{K}},$

where in (b) we have substituted
$\alpha_{b,k}=\frac{\alpha_{b,0}}{\sqrt{k+1}}$ and
$\beta_{b,k}=\frac{\beta_{b,0}}{\sqrt{k+1}}$, and (c) is by
$\sum_{k=0}^{K-1}\frac{1}{(k+1)^{3/2}}\leq 2$,
$\sum_{k=0}^{K-1}\frac{1}{k+1}\leq\log(K)$, and
$\sum_{k=0}^{K-1}\frac{1}{(k+1)^{1/2}}\leq\sqrt{K}$. Similar to the proof of
Theorem 3, we choose
$\alpha_{l,0}\sim\alpha_{b,0}\sim\mathcal{O}(\kappa^{-3})$,
$T\sim\mathcal{O}(\kappa)$, and $\beta_{b,0}\sim\mathcal{O}(\kappa^{-1})$ to
obtain $\frac{1}{K}\sum_{k=0}^{K-1}\mathbb{E}[\lVert\nabla
F(x^{k})\rVert^{2}]\leq\mathcal{O}(\frac{\kappa^{3}}{\sqrt{K}}+\frac{\kappa^{2}\log(K)}{\sqrt{K}})$.
∎

##  Appendix B Additional Experiment Results

This section is organized as follows. First, we discuss synthetic quadratics
experiments. Second, we provide more details on the sensitivity of the
algorithms to the choices of $\delta$ in (14), on the reset procedure, and on
the search cost of BiSLS. Third, we compare the empirical performance of
1-sample vs 2-samples implementations of our algorithms for single-level
convex and bi-level optimization. Some additional results for hyper-
representation learning and data distillation experiments are also presented.
We run $5$ independent runs for all our experiments.

###  B.1 Synthetic Quadratics

The experiments on quadratic functions are adapted from Loizou et al. [29].
The training objective is as follows:

$f(x)=\frac{1}{2}\,(x-x_{1}^{*})^{T}H_{1}(x-x_{1}^{*})\,+\,\frac{1}{2}\,(x-x_{2}^{*})^{T}H_{2}(x-x_{2}^{*}),$

where $H_{i}$ ($i=1,2$) are positive definite. The optimal solutions
$x_{i}^{*}$ ($i=1,2$) are generated randomly from a standard normal
distribution. Specifically, $H_{i}$ is defined as follows:

$H_{i}=O^{T}\cdot\text{Diag}(\log(1+\lambda_{i}))\cdot O,\quad i=1,2,$

where $O$ and $\lambda_{i}$ are taken from the spectral decomposition of
$P^{T}P$, and $P$ is generated from the standard normal distribution. Figure 8
shows the convergence of various algorithms with different starting points.
Interestingly, both SPSB with either 1 sample or 2 samples (1 sample for
computing the gradient and the other for computing the step size) converge to
the optimal solution (labelled with a star).

Figure 8: Iterate trajectories of different starting points for the synthetic
quadratic experiments.

###  B.2 Sensitivity of $\delta$, reset, and search cost

In this section, we discuss the effects of $\delta$ in (14), $\eta$ in
Algorithm 2, and comment on the search cost of the options in the Reset
Algorithm 2. Recall that the line-search condition (8) assumes that we can
find a largest $\gamma_{k}\leq\gamma_{b,0}$ to satisfy it. However, in
practice, we apply a backtracking procedure, i.e.
$\gamma_{k}=\gamma_{k}*w,0<w<1$, until $\gamma_{k}$ satisfies (8). Therefore,
the found learning rate is not guaranteed to be the largest. Nonetheless, we
assume that $\gamma_{k}$ is the largest to simplify our analysis given above
(similar arguments apply to line-search at both upper and lower-level in the
bi-level optimization). The experiments in this section are based on hyper-
representation learning [42]. In this case, the objective of the induced
bilevel-optimization problem can be written as:

$\displaystyle\min_{w}F(w)=\frac{1}{2D_{X_{1}}}\lVert\tilde{f}(X_{1};w)c^{*}(w)-Y_{1}\rVert^{2}$
$\displaystyle s.t.\quad
c^{*}(w)=\operatorname*{argmin}_{c}\frac{1}{2D_{X_{2}}}\lVert\tilde{f}(X_{2};w)c-Y_{2}\rVert^{2}+\frac{\lambda}{2}\lVert
c\rVert^{2}\,,$

where $(X_{1},Y_{1})$ and $(X_{2},Y_{2})$ are validation and training data
sets with sizes $D_{X_{1}}$ and $D_{X_{2}}$, respectively;
$\tilde{f}(\cdot;w)$ are the embedding layers of the model parameterized by
$w$; and, $c$ is the classification layer. Moreover, we use conjugate gradient
methods (CG) [17, 18] to solve the linear system when computing the
hypergradient for hyper-representation learning experiments.

##### Reset

While Algorithm 2 (reset) can be applied to both upper and lower-level
problems, we focus our discussions here on the upper-level learning rate
($\alpha_{k}$). This is because we empirically find it to be more critical for
the convergence performance (see Figure 6(a)). As shown in Algorithm 2, reset
has $3$ options. Options $1$, $2$, and $3$ search starting from
$\alpha_{b,0}$, $\alpha_{k-1}$, and $\eta\alpha_{k-1}$, at iteration $k$
respectively. Option $1$ has the highest search cost as it always starts from
the same initial upper bound ($\alpha_{b,0}$). Option $2$ ensures the
monotonicity of the learning rate due to
$\alpha_{k}\leq\alpha_{b,k}=\alpha_{k-1}$. Option $3$ chooses the search
starting point at iteration $k$ ($\alpha_{b,k}$) by multiplying the previous
learning rate ($\alpha_{k-1}$) by a factor $\eta\geq 1$. As in the single-
level convex case where monotonicity in the step size can potentially lead to
slow convergence (see Figure 3), we again observe that monotonicity in the
upper learning rate (i.e. option 2) leads to poorer performance when compared
against options 1 or 3 as shown in Figure 9. Finally, we compare the
performance of different $\eta$s in option $3$ (note that $\eta=1$ in option 3
is equivalent to option 2). We observe in Figure 10 that different $\eta$s
perform equally well. This shows the robustness of our algorithm to the choice
of $\eta$. As mentioned previously, the choice of option $3$ over option $1$
are due to $2$ reasons: (a) reduced search cost; (b) provides an overall non-
increasing and non-monotonic trend of upper bound $\alpha_{b,k}$. We discuss
search cost of different $\eta$s in option 3 below.

##### Search Cost

We investigate the line-search cost based on options $1$, $2$, and $3$ in
reset (Algorithm 2) for the upper-level problem. For option 1, the search cost
in terms of number of search rounds (i.e. number of evaluations using (14))
per iteration for BiSLS-SGD and BiSLS-Adam is $89\pm 15$ and $115\pm 16$,
respectively. The search cost of BiSLS-Adam is higher than that of BiSLS-SGD
because the feasible learning rate range for Adam is typically smaller than
SGD at the upper-level. Based on the results in Figure 11, we observe the use
of option 3 in reset can significantly reduce the search cost for both BiSLS-
Adam and BiSLS-SGD. For example, choosing $\eta=2$ in option 3 results in an
average upper-level search cost of only $\sim 9$ rounds per iteration, which
is much smaller than that of option 1. As we have shown in Figure 9, options 1
and 3 have nearly the same performance. Therefore, option 3 is an efficient
algorithm that maintains good performance while reducing computation cost.
Option 2 has the lowest search cost ($\sim 4$ rounds per iteration). However,
its performance is not as good as option 1 or 3 as observed in Figure 9.
Moreover, the average lower-level search cost is only $1$ round per iteration
when option 1 is used (see Figure 11(b)).

##### Sensitivity on $\delta$

As mentioned in Sec 2, due to the stochastic error in hypergradient
computation, further complicated by the approximation error of $y^{*}(x)$ (see
(14)), a learning rate is not guaranteed to be found in the bi-level case.
Specifically, this is in contrast to the single-level convex problems. To
avoid this, we introduce in (14) a $\delta$ slack to give some tolerance to
such errors. Here, we give a thorough investigation of the effects of $\delta$
on performance. We vary its magnitude across $6$ orders for both reset options
$1$ and $3$ (see Algorithm 2 and discussions on reset above). We observe that
despite a large difference on the magnitudes of $\delta$, they all share very
similar performance for both BiSLS-SGD and BiSLS-Adam: see Figures 12 and 13.
We summarize the key fins in this section as follows: ① The option 3 in reset
has good empirical performance (outperforms option 2) and is an effective way
to reduce search cost (Figure 9, 11); ② BiSLS is highly robust to different
choices of $\eta$ in option 3 and $\delta$ in (14) (Figure 10, 12, 13).

(a) Different Opt (BiSLS-SGD)

(b) Different Opt (BiSLS-Adam)

Figure 9: Validation loss against iterations with search options $1$, $2$, and
$3$ for the upper-level learning rate. The results for BiSLS-SGD and BiSLS-
Adam are in (a) and (b), respectively. For the lower-level search, we fix it
option 1 with $\beta_{b,0}=100$. Results are based on hyper-representation
learning.

(a) Different $\eta$s (BiSLS-SGD)

(b) Different $\eta$s (BiSLS-Adam)

Figure 10: Validation loss against iterations for different $\eta$s based on
reset option 3. Results for BiSLS-SGD are given in (a) and for BiSLS-Adam are
given in (b). Note that $\eta=1$ in reset option 3 is equivalent to reset
option 2. For the lower-level search, we fix it option 1 with
$\beta_{b,0}=100$. Results are based on hyper-representation learning.

(a) Upper-level search cost

(b) Lower-level search cost

Figure 11: (a) Upper-level search cost measured in terms of average number of
search rounds per iteration against different $\eta$s in reset (Algorithm 2).
(b) Lower-level search cost measured in the same way as upper-level against
different lower-level search starting points ($\beta_{b,0}$). The lower-level
search is done with option 1 (see above for discussions about these options).

(a) Different $\delta$s/Opt-1 (BiSLS-SGD)

(b) Different $\delta$s/Opt-1 (BiSLS-Adam)

Figure 12: Validation loss against iterations for different $\delta$s based on
reset option 1. Results for BiSLS-SGD (a) and for BiSLS-Adam (b). For the
lower-level search, we fix it option 1 with $\beta_{b,0}=100$. Results are
based on hyper-representation learning.

(a) Different $\delta$s/Opt-3 (BiSLS-SGD)

(b) Different $\delta$s/Opt-3 (BiSLS-Adam)

Figure 13: Validation loss against iterations for different $\delta$s based on
reset option 3 ($\eta=10$). Results for BiSLS-SGD (a) and for BiSLS-Adam (b).
For the lower-level search, we fix it option 1 with $\beta_{b,0}=100$. Results
are based on hyper-representation learning.

###  B.3 Data distillation objective and additional results

We let $\mathcal{L}_{S}(w)$ denote the loss evaluated on dataset $S$ with
model weights $w$. The objective of _data distillation_ can be expressed as a
BO problem as follows:

$D^{*}\,=\,\operatorname*{argmin}_{D}\mathcal{L}_{\tilde{V}}(w^{*}(D))\quad\text{s.t.}\quad
w^{*}(D)\,=\,\operatorname*{argmin}_{w}\mathcal{L}_{D}(w),$

where $\tilde{V}$ is of the same size as $D$ and subsampled from the entire
(original) dataset $V$. The solution $D^{*}$ is the distilled data, e.g. $9$
MNIST digits each corresponding to a different label. In figure 14(a), we show
the performance of BiSPS for different values of $\alpha_{b,0}$ in comparison
with BiSLS-SGD, and observe that BiSLS-SGD has better performance. In 14(b),
we show the results when we increase the number of lower-level iterations (T)
from 20 to 50. As observed for $T=20$ (in Figure 7), BiSLS-SGD here also
outperforms a fine-tuned Adam or SGD.

(a) Different $\alpha_{b,0}$s (BiSPS)

(b) T = 50, Identity

Figure 14: Validation loss against iterations. (a) Comparison between BiSPS
with different $\alpha_{b,0}$s and BiSLS-SGD. (b) Comparison between BiSLS-SGD
to fine-tuned Adam/SGD ($\beta_{k}$ fixed at $10^{-4}$). Inverse Hessian in
(2) treated as the Identity [30] when computing the hypergradient. Recall that
T is the total number of lower-level iterations and we have shown the results
for $T=20$ in Figure 7(b).

###  B.4 1-sample or 2-samples versions of algorithms for convex and bi-level
optimization

We provide additional results to compare the performance of 1-sample and
2-samples (one for computing the gradient and the other for computing the step
size) versions of our algorithms for SPSB and BiSPS used for single-level and
bi-level optimization, respectively. In the single-level case (Figure 15), we
observe that 2-samples SPSB performs just as well as 1-sample SPSB .
Interestingly, we observe that their step sizes also follow a similar pattern.
That is: an initial increase followed by a regime where
$\gamma_{k}=\frac{f_{i_{k}}(x^{k})-l_{i_{k}}^{*}}{c\lVert\nabla
f_{i_{k}}(x^{k})\rVert^{2}}$ is frequently used, and eventually changes to
decaying-step SGD. This seems to also match with Theorem 2 where a transition
point for SPSB
($k_{0}=\max\\{1,\lceil\gamma_{0}/\omega\rceil-1\\},w=\frac{1}{2cL_{\max}},\gamma_{0}\geq\frac{1}{\mu}$)
is predicted. At the same time, we also note that (perhaps, unsurprisingly)
the 2-samples version seems to have a slightly more oscillatory behavior than
the 1-sample version as shown in Figure 15. SLSB with either 1-sample or
2-samples also result in a similar performance and step size. Overall, despite
the requirements of Theorems 1 and 2 for a 2-samples assumption, the empirical
performance of 1-sample and 2-samples for either SPSB or SLSB appears to be
very similar. Moving on to the bi-level case, recall that Theorems 3 and 4
require the 2-samples assumption (i.e., $\alpha_{k}$ independent of
$h_{f}^{k}$) for the upper-level learning rate. We empirically verify this
assumption with both hyper-representation learning and data distillation
experiments. For hyper-representation learning experiments in Figure 16, BiSPS
with either 1-sample or 2-samples for different values of $\alpha_{b,0}$ show
similar performance. In fact, for $\alpha_{b,0}=0.1$ we even observe that the
2-samples variant outperforms the 1-sample BiSPS. For data distillation
experiments in Figure 17, the performances of 1-sample and 2-samples BiSPS are
similar to each other when $\alpha_{b,0}=10.0$ or $\alpha_{b,0}=50.0$. In
general, the performance difference between 1-sample and 2-samples in the
single-level or bi-level settings is small.

(a) Train loss

(b) Step size

Figure 15: Binary linear classification on w8a dataset using logistic loss
[3]. Train loss (left) and step size (right) against iterations. We choose
$\gamma_{b,0}=1000$ for all algorithms. The upper bound for either SPSB or
SLSB decays as $\gamma_{b,k}=\frac{\gamma_{b,0}}{\sqrt{k+1}}$. For decaying-
step SGD, the learning rate schedule is $\frac{\gamma_{b,0}}{\sqrt{k+1}}$.
Figure 16: Comparison between BiSPS (2-samples), BiSPS (1-sample) and
decaying-step SGD. Experiments are based on hyper-representation learning. For
either version of BiSPS, the lower-level learning rate ($\beta_{k}$) is fixed
at $10$. The hypergradient is computed using conjugate gradient [17].

Figure 17: Comparison between BiSPS (2-samples), BiSPS (1-sample) and
decaying-step SGD. Experiments are based on data distillation. For either
version of BiSPS, the lower-level learning rate ($\beta_{k}$) is fixed at
$10^{-4}$. The Inverse Hessian in (2) is treated as the Identity when
computing the hypergradient [30].