text:
\section{Introduction}
\label{Introduction}
(Stochastic) bilevel optimization is a widely confronted
problem in machine learning with various applications such as
meta-learning \citep{finn2017model, bertinetto2018meta,rajeswaran2019meta}, hyper-parameter optimization \citep{franceschi2018bilevel,shaban2019truncated,baydin2017online,bergstra2011algorithms,luketina2016scalable}, reinforcement learning \citep{hong2020two}, and few-shot learning \citep{koch2015siamese,santoro2016meta,vinyals2016matching}.
The basic form of this problem can be  defined as follows
\begin{equation}\label{pb1}
  \begin{aligned}
  \min _{\mathbf{x} \in \mathbb{R}^{d_{1 }}}  R(\mathbf{x})&= F\left(\mathbf{x}, \mathbf{y}^{*}(\mathbf{x})\right):=\mathbb{E}_{\xi}\left[f\left(\mathbf{x}, \mathbf{y}^{*}(\mathbf{x}) ; \xi\right)\right] \\
  \text { s.t. }  \mathbf{y}^{*}(\mathbf{x})&=\arg \min _{\mathbf{y} \in \mathbb{R}^{d_{2}}}\left\{G(\mathbf{x}, \mathbf{y}):=\mathbb{E}_{\zeta}[g(\mathbf{x}, \mathbf{y} ; \zeta)]\right\},
  \end{aligned}
  \end{equation}
where $f$ : $\mathbb{R}^{d_1} \times \mathbb{R}^{d_{2 }} \rightarrow \mathbb{R} $ and $g$ : $\mathbb{R}^{d_1} \times \mathbb{R}^{d_{2 }} \rightarrow \mathbb{R} $ are two continuously differentiable loss functions with respect to $\mathbf{x}$ and  $\mathbf{y}$. Problem (\ref{pb1}) has an optimization hierarchy of two levels, where the outer-level objective function $f$ depends on the minimizer of the inner-level objective function $g$.

 Due to its importance,

 the above bilevel optimization problem has received considerable attention in recent years. A natural way to solve problem (\ref{pb1}) is to apply alternating stochastic gradient updates with approximating $\nabla_{\mathbf{y}} g(\mathbf{x}, \mathbf{y})$ and $\nabla f(\mathbf{x} , \mathbf{y})$, respectively. Briefly speaking, previous efforts mainly examined two types of methods to perceive an approximate solution that is close to the optimum $y^* (\mathbf{x})$. One is to utilize the single-timescale strategy \citep{chen2021single, guo2021randomized, khanduri2021near,hu2022improved}, where the updates for $\mathbf{y}$ and $\mathbf{x}$ are carried out simultaneously. The other one is to apply the two-timescale strategy \citep{ghadimi2018approximation, ji2021bilevel,hong2020two,pedregosa2016hyperparameter}, where the  update of $\mathbf{y}$ is repeated multiple times to achieve a more accurate approximation before conducting the  update of $\mathbf{x}$.

 While there is a long list of work on bilevel optimization,
most of the existing work only focuses on either analyzing its convergence behaviors \citep{ghadimi2018approximation, hong2020two, ji2021bilevel} or improving its convergence rate, based on the convexity and the smoothness properties of  $f(\cdot, \cdot)$ and/or $g(\cdot,\cdot)$ \citep{liu2020generic, li2020improved}.

 Contrarily,
 only little effort is devoted to

 understanding the generalization behavior of the problem. To the best of our knowledge, only one recent work considers
 generalization

 \citep{bao2021stability}, which presents the first expected stability bound

 for the bilevel optimization problem.  Unfortunately, the problem is still far from well-understood and   there are several undesirable issues in this work: (1) Their study is only for the hyper-parameter optimization problem, which is only a special instance of bilevel optimization. Thus, the stability for general bilevel optimization is still unknown; (2) Their result is only for the uniform stability (which could be deduced from argument stability with certain conditions, see Definition \ref{def:stab} for details), leaving the analysis of other stronger definitions of algorithmic stability open; (3) \cite{bao2021stability} only considers a general setting for inner function. However, there are some technical flaws in their analysis (see Remark \ref{remark_bao} for detailed discussions) and
 their bound is quite loose in some cases. For example, in the case where the inner function is strongly convex, it is known that even in the classical stochastic convex optimization, strongly convex functions have better stability than the convex ones \citep{hardt2016train}. Thus, it is possible to provide improved stability bounds in some cases via different and novel analysis; (4) Their work only considers a special class of optimization algorithms, which allows the outer level parameters to be updated continuously but needs to reinitialize the inner level parameters before each iteration in the inner level. However, such algorithms  are not commonly used in practice due to their inefficiency;  (5) Their experiments take only hyper-parameter optimization into consideration and neglect other applications in the bilevel optimization instances.

To address all the aforementioned issues, we give in this paper a thorough analysis on the generalization behaviors of first-order (gradient-based) methods for general bilevel optimization problem. We
employ the recent advances of algorithmic stability to investigate the generalization behaviors
in different settings. Specifically,  our main contributions can be summarized as follows:
\begin{itemize}
    \item Firstly, we establish a fundamental connection between generalization error and different notations  of algorithmic stability (argument stability and uniform stability) for
     bilevel optimization problem in both  expectation and high probability forms. Specifically, we

    show that the high probability form of the generalization error bound can be improved from

   $\bigO(\sqrt{n})$ to $\bigO(\log n)$ compared with the result in \citet{bao2021stability}.

    \item Next, we present the first stability bounds for gradient-based methods with either single-timescale or two-timescale update strategy under different standard settings. In detail, we consider the settings of
    strongly-convex-strongly-convex (SC-SC), convex-convex (C-C), and nonconvex-nonconvex (NC-NC), and further extend our analysis to a particular nonconvex-strongly-convex (NC-SC) setting that is widely appeared in practice. Especially, for the SC-SC and C-C cases with single-timescale update strategy we significantly improve the generalization bounds compared with \citep{bao2021stability} (see Table \ref{tab:main results} for details).  To the best of our knowledge, our work is the first thorough generalization analysis for general  bilevel optimization problem.
    \item Thirdly, we provide the first generalization bounds for the case where

    both the outer and inner level parameters are subject to continuous (iterative) changes. Compared to the previous work \citep{bao2021stability}, our work does not need the reinitialization step before each iteration in the inner level

    and hence our algorithm can carry over the last updated inner level parameters, which is more general and practical.
    \item Finally, we conduct empirical studies to corroborate our theories
    via meta-learning and hyperparameter optimization, which are two  applications of bilevel optimization.
\end{itemize}
Due to space limitations, all the proofs and additional experiments are included in Appendix.
\begin{table*}[t]
\captionsetup{font=scriptsize}
\centering
\begin{center}
\footnotesize
\scriptsize
\begin{sc}
\scalebox{0.77}{
\begin{tabular}{ccccccccc}
\toprule
\multirow{2}{*}{Reference} & \multicolumn{4}{c}{stability bounds in various settings} & \multicolumn{2}{c}{stability}                         & \multicolumn{2}{c}{consistency}         \\ \cmidrule(lr){2-5}  \cmidrule(lr){6-7}  \cmidrule(lr){8-9}
                           & SC-SC         & C-C        & NC-NC        & NC-SC        & argument                  & uniform                   & Outer            & inner                \\ \hline
UD \citep{bao2021stability}                         & \multicolumn{4}{c}{$\bigO(K^{\kappa_1}/m_1)$} & & \checkmark & \checkmark &                           \\
SSGD (This work)                  & $\bigO(1/m_1)$ & $\bigO(C^{K/2} /m_1)$ & $\bigO(K^{\kappa_2}/m_1)$ & $\bigO(K^{\kappa_3}/m_1)$ &\checkmark &       \checkmark                    & \checkmark & \checkmark \\
TSGD (This work)               & $\bigO((\kappa_4 T)^K/m_1)$ & $\bigO((T)^K/m_1)$ & $\bigO( T^{1-\kappa_4} K^{\kappa_5}/m_1)$ & $\bigO( T^{1-\kappa_6} K^{\kappa_6}/m_1)$           & \checkmark & \checkmark & \checkmark & \checkmark \\
\bottomrule
\end{tabular}}
\caption{
Comparison of related work with ours.
$\kappa_i$: a constant such that $\kappa_i\in(0,1)$ for all $i$ above; $T$: inner iterations; $K$: outer iterations; $m_1$: size of outer dataset. UD indicates the unrolled differentiation (Algorithm \ref{alg:UD}). SSGD and TSGD stand for Algorithm \ref{alg:SSGD} and Algorithm \ref{alg:TSGD}, the single-timescale and two-timescale methods, via stochastic gradient descent. We provide a comparison of TSGD and UD in the Appendix \ref{comparision between ud and tsgd}. Consistency indicates whether it is allowed to update inner-level or outer-level parameters continuously.
}
\label{tab:main results}
\end{sc}
\end{center}
\vspace{-0.5cm}
\end{table*}

\subsection{Related work}
\label{related_work}

Research at the interface between generalization and the bilevel problem can be roughly classified into two categories. The first one includes all the research on bilevel optimization. In recent decades, extensive studies have been done on this topic, which suggests that bilevel optimization has a wide range of applications in machine learning such as
hyper-parameter optimization \citep{franceschi2018bilevel,lorraine2018stochastic,okuno2021lp}, meta learning \citep{bertinetto2018meta,rajeswaran2019meta,soh2020meta} and reinforcement learning \citep{yang2018convergent,tschiatschek2019learner}. Most of the existing work studies the problem  from an optimization perspective. For example,
\citet{ghadimi2018approximation,ji2021bilevel} provide the convergence rate analysis based on the nonconvex-strongly-convex assumption for the two functions $f(\cdot, \cdot)$ and $g(\cdot,\cdot)$. \citep{grazzi2020iteration} considers the iteration complexity for hypergradient computation. \citep{liu2020generic, li2020improved} present an asymptotic analysis for the convex-strongly-convex setting. Perhaps the most related one to ours
from the generalization standpoint ({\em i.e.,} the expectation of population risk and empirical risk) is \citet{bao2021stability}. \textbf{Table \ref{tab:main results}} is the comparison between our work and  the results in \citep{bao2021stability}.  We can see that
our analysis is more comprehensive as we consider various settings. Firstly, in this paper we consider the general bilevel problem defined in Problem (\ref{pb1}), while \citet{bao2021stability} only concentrates  on the hyper-parameter optimization, which is just an instance of  bilevel optimization problem (to be detailed in Section \ref{experiment} and Appendix \ref{additional experiment}).
Secondly, our work adopts a stronger version of stability called {\em argument stability}, which implies the previously used {\em uniform stability} if the function is sufficiently smooth. Furthermore, our work does not need to reinitialize the inner-level parameters and allows them to carry over their last updated parameters at each time updating the inner level. This indicates that
$\mathbf{y}$ in the inner level is updated iteratively  and depends on the current parameter of $\mathbf{x}$, which is more common and efficient in practice.

The second category includes all the work on stability analysis. There is a long list of research on
stability and generalization \citep{bousquet2002stability, mukherjee2006learning, shalev2010learnability}. \citet{bousquet2002stability} first introduces the notion of \textit{uniform stability} and establishes the first framework of stability analysis. \citet{hardt2016train} later extends the stability analysis  to iterative algorithms based on stochastic gradient methods for the vanilla stochastic optimization. After that, there are  subsequent studies on generalization analysis for various problems via algorithmic stability, such as minmax problems \citep{lei2021stability, farnia2021train,zhang2021generalization} and pairwise learning \citep{yang2021simple, lei2020sharper, xue2021differentially,huai2020pairwise}. However, it is notable that due to the additional stochastic function in the constraint in the bilevel optimization, all the previous techniques and results cannot be applied to our problem. Although the generalization analysis of minmax optimization is somewhat similar to ours, it involves only one objective function $f$ and a single level in algorithms for typical minmax optimization problems, while in the bilevel optimization algorithms there is an inner level and an outer level, which is considerably more challenging.

\section{Preliminaries}
\label{preliminaries}
\subsection{Definitions and Assumptions}
In the following, we give some necessary definitions and assumptions that are widely used in bilevel optimization \citep{ghadimi2018approximation,ji2021bilevel,khanduri2021near} and generalization analysis \citep{hardt2016train,lei2021stability}.
\begin{definition}[Joint Lipschitz Continuity]
    A function $f(\mathbf{x},\mathbf{y})$ is jointly $L$-Lipschitz over $\mathbb{R}^{d_1} \times \mathbb{R}^{d_2}$, if for all $\mathbf{x} \in \mathbb{R}^{d_1}, \mathbf{y} \in \mathbb{R}^{d_2}$, the following holds,   $|f(\mathbf{x},\mathbf{y})-f(\mathbf{x}^{\prime},\mathbf{y}^{\prime})| \leq L\sqrt{\left\|\mathbf{x}-\mathbf{x}^{\prime}\right\|_{2}^{2}+\left\|\mathbf{y}-\mathbf{y}^{\prime}\right\|_{2}^{2}}$.
\end{definition}

\begin{definition}[Smoothness]
     A function $f$ is $l$-smooth over a set $S$ if for all $u, w \in S$ the following is true, $\|\nabla f(u)-\nabla f(w)\| \leq l\|u-w\|$.
\end{definition}

\begin{definition}[Strong Convexity]
     A function $f$ is $\mu$-strongly-convex over a set $S$, if for all $u, w \in S$, the following holds, $f(u)+\langle\nabla f(u), w-u\rangle+\frac{\mu}{2}\|w-u\|^{2} \leq f(w)$.
\end{definition}

\begin{assumption}[Inner-level Function Assumption]\label{asu1}
  We assume the inner stochastic function  $g(\mathbf{x},\mathbf{y})$ in (\ref{pb1}) satisfies the following:\\
    \subasu $g(\mathbf{x},\mathbf{y})$ is jointly $L_g$-Lipschitz for any $\mathbf{x}\in \mathbb{R}^{d_1}$ and $\mathbf{y} \in \mathbb{R}^{d_2}$.\\
    \subasu $g(\mathbf{x},\mathbf{y})$ is continuously differentiable and  $l_g$-smooth for any $(\mathbf{x},\mathbf{y}) \in \mathbb{R}^{d_1} \times \mathbb{R}^{d_2}$.

\end{assumption}

\begin{assumption}[Outer-level Function Assumption]\label{asu2}
  We assume the outer stochastic function  $f(\mathbf{x},\mathbf{y})$ in (\ref{pb1}) satisfies the following:\\
    \subasu $f(\mathbf{x},\mathbf{y})$ is jointly $L_f$-Lipschitz for any $ \mathbf{x} \in \mathbb{R}^{d_1}$ and $\mathbf{y} \in \mathbb{R}^{d_2}$.\\
    \subasu $f(\mathbf{x},\mathbf{y})$ is continuously differentiable and $l_f$-smooth for any $(\mathbf{x},\mathbf{y}) \in \mathbb{R}^{d_1} \times \mathbb{R}^{d_2}$.

\end{assumption}

\subsection{Problem Formulation}\label{sec:problem}
Given two distributions $\mathbb{D}_1$ and $\mathbb{D}_2$, in the (stochastic) optimization problem we aim to find the minimizer of Problem (\ref{pb1}). However, since the distributions are often unknown, in practice we only have two finite-size datasets  $D_{m_1}= \{\xi_i\mid i=1,...,m_1\}\sim \mathbb{D}_1^{m_1}$ and $D_{m_2}= \{\zeta_i\mid i=1,...,m_2\} \sim \mathbb{D}_2^{m_2}$, where each $\xi_i$ and $\zeta_i$ are i.i.d. sampled from $\mathbb{D}_1$ and $\mathbb{D}_2$, respectively. Based on these datasets, we will design some (randomized) algorithm ${\bf A}$ with output ${\bf A}(D_{m_1}, D_{m_2})= (\mathbf{x}, \mathbf{y})\in \mathbb{R}^{d_1}\times \mathbb{R}^{d_2}$. Our goal is to investigate the generalization behavior of such output. Note that although there are two stochastic functions in the bilevel optimization problem, we only care about the generalization of the outer-level one since it is the one that we prefer to minimize.

Below we define the generalization error to measure the generalization behavior.
Given  distribution $\mathbb{D}_1$  and  a finite data $D_{m_1}\sim \mathbb{D}_1^{m_1} $, the population risk function $R(\mathbf{x},\mathbf{y},\mathbb{D}_1)$
of $\mathbf{x}$, $\mathbf{y}$ on $\mathbb{D}_1$ is defined as $R(\mathbf{x},\mathbf{y},\mathbb{D}_1) := \mathbb{E}_{{\xi} \sim \mathbb{D}_1}\left[f\left(\mathbf{x}, \mathbf{y}(\mathbf{x}) ; \xi\right)\right]$, and its empirical risk function  on $D_{m_1}$ is
 $R_s(\mathbf{x},\mathbf{y}, {D}_{m_1}) = \frac{1}{m_1} \sum_{i=1}^{m_1}\left[f\left(\mathbf{x}, \mathbf{y}(\mathbf{x}) ; \xi_i\right)\right]$. Moreover, for a fixed hyperparameter $\mathbf{x} \in \mathbb{R}^{d_{1 }} $ and  $\mathbf{y}(\mathbf{x}) \in \mathbb{R}^{d_{2}}$ ( note that $\mathbf{y}(\mathbf{x})$ might be dependent on $\mathbf{x}$),
  we define the difference between the population risk and the empirical risk over $(\mathbf{x}, \mathbf{y}(\mathbf{x}))$ as the {\em bilevel generalization error of $(\mathbf{x}, \mathbf{y}(\mathbf{x}))$}:
  $\mathbb{E}_s[R(\mathbf{x},\mathbf{y})-R_s(\mathbf{x},\mathbf{y})]$, where $\mathbb{E}_s$ denotes the expectation of $D_{m_1}\sim \mathbb{D}_1^{m_1}$. When there is no ambiguity, we simplify thereafter the notations as follows:  $R(\mathbf{x},\mathbf{y},\mathbb{D}_1) =R(\mathbf{x},\mathbf{y}) $ and $R_s(\mathbf{x},\mathbf{y},D_{m_1}) = R_s(\mathbf{x},\mathbf{y})$.

  Our goal is thus to analyze the {\em bilevel generalization error of the output of algorithm ${\bf A}(D_{m_1}, D_{m_2})$} based on $D_{m_1}$ and  $D_{m_2}$. Since the generalized error depends on the algorithm itself, in the following we will introduce the algorithms to be considered in this paper.

Most of the existing algorithms adopt the following idea:
first approximate $\mathbf{y}^*$ on $D_{m_2}$ for a given parameter $\mathbf{x}$ in the inner level and then seek the hyperparameter $\mathbf{x}^*(D_{m_1},D_{m_2})$ with corresponding hypothesis $\mathbf{y}^*(\mathbf{x}^*(D_{m_1},D_{m_2}),D_{m_2})$ by the below estimation:
\begin{equation}
\label{eq:approximation}
\begin{aligned}
     \hat{\mathbf{x}}(D_{m_1},D_{m_2}) & \approx \arg \min_{\mathbf{x}} R_s(\mathbf{x},\hat{\mathbf{y}}(\mathbf{x},D_{m_2}),D_{m_1}),\\
      \text{where} \quad \hat{\mathbf{y}}(\mathbf{x},D_{m_2}) & \approx \arg \min_{\mathbf{y}}G_s(\mathbf{x},\mathbf{y},D_{m_2}),
\end{aligned}
\end{equation}
      where $G_s(\mathbf{x},\mathbf{y},D_{m_2})$ is the empirical risk of $G(\mathbf{x},\mathbf{y})$ over $D_{m_2}$, {\em i.e.,}  $G(\mathbf{x},\mathbf{y},D_{m_2}) = \frac{1}{m_2} \sum_{i=1}^{m_2}g\left(\mathbf{x}, \mathbf{y}(\mathbf{x}) ; \zeta_i\right)$. Most of the current gradient-based (first-order) algorithms for approximating (\ref{eq:approximation}) can be categorized into two classes:  single-timescale methods and two-timescale methods. The single-timescale method performs the updates for $\mathbf{y}$ and $\mathbf{x}$ simultaneously via stochastic gradient descent (SGD), while the two-timescale method updates $\mathbf{y}$ multiple times  before updating $\mathbf{x}$ (via stochastic gradient descent). As there are numerous approaches for both classes (see Related Work section for details), in this paper we will analyze the generalization behaviors for the most classical and standard one in each class, {\em i.e.,} single-timescale SGD (SSGD; Algorithm \ref{alg:SSGD}) and two-timescale SGD (TSGD; Algorithm \ref{alg:TSGD}). There is a long list of work \citep{chen2021single}, \citep{ghadimi2018approximation,ji2021bilevel}

      based on either SSGD or TSGD.

\begin{figure}[]
\begin{minipage}{0.5\textwidth}
\vspace{-33pt}
\begin{algorithm}[H]
  \begin{algorithmic}[1]
  \caption{Single-timescale SGD (SSGD)}
  \label{alg:SSGD}
    \STATE {\bfseries Input:} number of iterations $K$, step sizes $\alpha_x$, $\alpha_y$, initialization $\mathbf{x}_0$,$\mathbf{y}_0$, Datasets $D_{m_1}$ and $D_{m_2}$
    \STATE {\bfseries Output:} $\mathbf{x}_K$, $\mathbf{y}_K$
    \FOR {$k=0$ {\bfseries to} $K-1$}
     \STATE Uniformly sample $i\in [m_2]$, $j\in [m_1]$
      \STATE $\mathbf{y}_{k+1} = \mathbf{y}_k - \alpha_y \nabla_{\mathbf{y}} g(\mathbf{x}_k, \mathbf{y}_k; \zeta_i) $
      \STATE $\mathbf{x}_{k+1} =\mathbf{x}_k -\alpha_{x} \nabla f(\mathbf{x}_k, \mathbf{y}_k;\xi_j) $
    \ENDFOR
    \STATE {\bf return} $\mathbf{x}_K$ and $\mathbf{y}_K$
  \end{algorithmic}
\end{algorithm}
\end{minipage}
\begin{minipage}{0.50\textwidth}
\begin{algorithm}[H]
  \begin{algorithmic}[1]
  \caption{Two-timescale SGD (TSGD) }
  \label{alg:TSGD}
    \STATE {\bfseries Input:} number of iterations $K$, step sizes $\alpha_x$, $\alpha_y$, initialization $\mathbf{x}_0$,$\mathbf{y}_0$
    \STATE {\bfseries Output:} $\mathbf{x}_K$, $\mathbf{y}_K$
    \FOR {$k=0$ {\bfseries to} $K-1$}
      \STATE $\mathbf{y}_{k}^{0} \gets \mathbf{y}_{k-1}^{T}$
      \FOR{$t=0$ {\bfseries to} $T-1$}
      \STATE Uniformly sample $i\in [m_2]$
      \STATE 	$\mathbf{y}_{k}^{t+1} = \mathbf{y}_k^{t} - \alpha_y \nabla_{\mathbf{y}} g(\mathbf{x}_k^t, \mathbf{y}_k^t; \zeta_i) $
      \ENDFOR
      \STATE Uniformly sample $j\in [m_1]$
      \STATE $\mathbf{x}_{k+1} =\mathbf{x}_k -\alpha_{x} \nabla f(\mathbf{x}_k, \mathbf{y}_k^T;\xi_j) $
    \ENDFOR
    \STATE {\bf return} $\mathbf{x}_K$, $\mathbf{y}_K^{T}$
  \end{algorithmic}
\end{algorithm}
\end{minipage}
\end{figure}

\section{Generalization and Stability for Bilevel Optimization}
\label{generalization and stability}
Algorithmic stability is one of the classical approaches to analyzing  the generalization bound for algorithms. Roughly speaking, the algorithmic stability of (randomized) algorithm ${\bf A}$ measures how the output of algorithm  ${\bf A}$ changes if we change one data sample in the input dataset.
While there are various notions of stability, most of the existing work on analyzing the stability of stochastic optimization, pairwise learning and minimax optimization focuses on the  uniform-stability \citep{bousquet2002stability} and the argument-stability \citep{liu2017algorithmic,lei2020fine}. Thus, we also adopt these two notions of stability for the bilevel optimization problem. Briefly speaking, uniform-stability focuses on the resulting change in population risk function, while the argument-stability considers the resulting change in arguments, {\em i.e.,} the output of the algorithm.
\begin{definition}[Algorithmic Stability]\label{def:stab}
    Let $\mathbf{A}: \mathbb{D}_1^{m_1} \times \mathbb{D}_2^{m_2} \mapsto \mathbb{R}^{d_1} \times \mathbb{R}^{d_2}$ be a randomized algorithm.
    \begin{itemize}
    \item[(a)] $\mathbf{A}$ is $\beta$-uniformly-stable if for all datasets $D_{m_1}, D^{\prime}_{m_1} \sim \mathbb{D}_1^{m_1}$ and $ D_{m_2} \sim \mathbb{D}_2^{m_2}$ such that $D_{m_1}$ and $D^{\prime}_{m_1}$ differ in at most one sample, we have the following for any $\xi \sim \mathbb{D}_1$:
    \begin{equation*}
        \mathbb{E}_{\mathbf{A}}[|f(\mathbf{A}(D_{m_1},D_{m_2}),\xi) - f(\mathbf{A}(D^{\prime}_{m_1},D_{m_2}),\xi) |] \leq \beta.
    \end{equation*}
    $\mathbf{A}$ is $\beta$-uniformly-stable with probability at least $1-\delta$ if we have the following for any $\xi \sim \mathbb{D}_1$ with probability at least $1-\delta$:
     \begin{equation*}
         \left| f(\mathbf{A}(D_{m_1},D_{m_2}),\xi) - f(\mathbf{A}(D^{\prime}_{m_1},D_{m_2}),\xi) \right|  \leq \beta.
    \end{equation*}
    \item[(b)] $\mathbf{A}$ is $\beta$-argument-stable in expectation if for all datasets $D_{m_1}, D^{\prime}_{m_1} \sim \mathbb{D}_1^{m_1}$ and $ D_{m_2} \sim \mathbb{D}_2^{m_2}$ such that $D_{m_1}$ and $D^{\prime}_{m_1}$ differ in at most one sample, we have:
    \begin{equation*}
         \mathbb{E}_{\mathbf{A}}[\|\mathbf{A}(D_{m_1},D_{m_2}) - \mathbf{A}(D^{\prime}_{m_1},D_{m_2}) \|_{2} ] \leq \beta.
    \end{equation*}
    \end{itemize}
\end{definition}
Note that the definition of uniform stability in expectation is the same as the definition in \citep{bao2021stability}. Thus, our other definitions can be considered as extensions of the previous stability for bilevel optimization. In the following, we present Theorem \ref{thm1_connection} as our first result, which shows a crucial relationship between generalization error and algorithmic stability for an algorithm ${\bf A}$.

\begin{theorem}\label{thm1_connection}
    Let $\mathbf{A}: \xi^{m_1} \times \zeta^{m_2} \mapsto \mathbb{R}^{d_1} \times \mathbb{R}^{d_2}$ be a randomized BO algorithm.
    \begin{itemize}
    \item[(a)] If $\mathbf{A}$ is $\beta$-uniform-stable in expectation, then the following holds for $D_{m_1} \sim \mathbb{D}_1^{m_1}$, $D_{m_2} \sim \mathbb{D}_2^{m_2}$:

    \begin{equation*}
        \mathbb{E}_{\mathbf{A},D_{m_1}}[R(\mathbf{A}(D_{m_1},D_{m_2}))-R_{s}(\mathbf{A}(D_{m_1},D_{m_2}))] \leq \beta.
    \end{equation*}

    \item[(b)] If $\mathbf{A}$ is $\beta$-argument-stable in expectation and Assumption \ref{asu2} holds, then the following holds for $D_{m_1} \sim \mathbb{D}_1^{m_1}$, $D_{m_2} \sim \mathbb{D}_2^{m_2}$:
    \begin{equation*}
        \mathbb{E}_{\mathbf{A},D_{m_1}}[R(\mathbf{A}(D_{m_1},D_{m_2}))-R_{s}(\mathbf{A}(D_{m_1},D_{m_2}))] \leq L_f\beta.
    \end{equation*}

    \item[(c)] Assume that $\left|f(\mathbf{x},\mathbf{y};\xi)\right| \leq M$ for some $M \geq 0$. If $\mathbf{A}$ is $\beta$-uniform-stable almost surely, then for $D_{m_1} \sim \mathbb{D}_1^{m_1}$, $D_{m_2} \sim \mathbb{D}_2^{m_2}$, the following holds with probability $1-\delta$:
    \begin{equation*}
        \begin{aligned}
            &\left|R(\mathbf{A}(D_{m_1},D_{m_2}))-R_{s}(\mathbf{A}(D_{m_1},D_{m_2})) \right|\\
            & \leq 2\beta+e\left(\frac{4 M}{\sqrt{m_{1}}} \sqrt{\log \frac{e}{\delta}}+12 \sqrt{2} \beta \lceil  \log _{2} m_1  \rceil \sqrt{\log \frac{e}{\delta}}\right).
        \end{aligned}
    \end{equation*}

    \end{itemize}
\end{theorem}
\begin{remark}
The above theorem suggests that the generalization error can be controlled by several notions of algorithmic stability. Part (a) and Part (b) show that the expectation of generalization error can be bounded by uniform stability and argument stability with the Lipschitz constant, respectively; Part (c) indicates that the generalization error for the algorithm is no more than
$\bigO(\beta\log(m_1) +1/\sqrt{m_1})$ with probability $1-\delta$. Compared with the existing work \citep{bao2021stability}, Theorem \ref{thm1_connection} considers argument stability additionally, which is a stronger notion of stability than uniform stability (since uniform stability can be deduced from argument stability with the condition that the function is sufficiently smooth).
Moreover, we use the McDiarmid's inequality and the equivalence of tails and moments for the random variable with a mixture of sub-gaussian and sub-exponential tails (Lemma 1 in \citet{bousquet2020sharper}), which provide a significantly improved high probability bound   in Part (c) ({\em i.e.,} improving from $\bigO(\beta\sqrt{m_1})$ in \citet{bao2021stability} to $\bigO(\beta \log m_1)$).
\end{remark}

\section{Stability Analysis for Bilevel Optimization Algorithms }
\label{main results}
Motivated by Theorem \ref{thm1_connection}, we can see that to analyze the generalization behaviors for any algorithm, it is sufficient to analyze its stability. As mentioned in the previous Section \ref{sec:problem}, we will consider the stability of SSGD and TSGD. For simplicity we let  SC-SC denote the case where $f$ and $g$ both  are strongly convex functions. C-C, NC-NC, and NC-SC are also denoted in a similar manner with ''C'' representing convex function and ''NC'' representing nonconvex function.

\subsection{Stability Bounds for Single-timescale SGD}
As we can see from Algorithm \ref{alg:SSGD}, SSGD updates $\mathbf{y}$ and $\mathbf{x}$ simultaneously. In the following we develop stability bounds for this algorithm in different settings.

\begin{theorem}\label{ssgd_stab}
    Suppose that Assumptions \ref{asu1} and \ref{asu2} hold and Algorithm $\mathbf{A}$ is SSGD with $K$ iterations:
    \begin{itemize}
        \item [(a)] Assume that Problem (\ref{pb1}) is SC-SC with strongly convexity parameters $\mu_f$ and $\mu_g$. Let $\alpha_x = \alpha_y$ (see  Lemma \ref{lemma:expansive} for details) be the step sizes. Denote $l = \max\{ l_f, l_g\}$. Then, $\mathbf{A}$ is $\beta$-argument-stable in expectation, where $$ \beta \leq \bigO\left(  \big(L_f^{2}+L_g^{2}\big)^{\frac{1}{2}}\left(m_1 \left(\mu_f+\mu_{g}- (\alpha_{x} l)^2/2 + 0.25 \right) \right)^{-1} \right).$$

        \item [(b)] Assume that Problem (\ref{pb1}) is C-C. Let $\alpha_x$, $\alpha_y$ be the step sizes. Then,  $\mathbf{A}$ is $\beta$-argument-stable in expectation, where
        $$\beta \leq \bigO \left( m_1^{-1}\sqrt{\left(\alpha_{x} L_f \right)^{2}+\left(\alpha_{y} L_g\right)^{2}}\left(2+2\max \left\{\left(\alpha_{x} l_f\right)^{2}, \left(\alpha_{y} l_g\right)^{2} \right\} \right)^{K/ 2} \right). $$

        \item [(c)] Assume that Problem (\ref{pb1}) is NC-NC. Let the step sizes satisfy $\max  \left\{ \alpha_x,\alpha_y  \right\} \leq c/k$ for some constant $c \geq 0$ and $l = \max \left\{ l_f, l_g\right\}$. Then, $\mathbf{A}$ is $\beta$-argument-stable in expectation, where
        $$\beta \leq \bigO \Big( (m_1  c l)^{-1}\left(2 c L_f \sqrt{l_f^2+ l_g^2}\right)^{\frac{1}{ c l+1}} \cdot K^{\frac{ c l}{ c l+1}} \Big),$$
    \end{itemize}
    where $l_f$, $l_g$ and $L_f$, $L_g$ are smoothness constants and Lipschitz constants for $f$, $g$, respectively.
\end{theorem}

\begin{remark}
Note that the above stability bounds are independent of the specific form of the objective function $f(\cdot,\cdot)$ and the exact form of the sample distribution $\mathbb{D}_1$, which are more reliant on the properties of the loss functions and sample size $m_1$, and the stability bounds in the C-C and NC-NC cases
are related to the number of iterations additionally. Specifically, Part(a) establishes a stability bound of $\bigO(1/m_1)$ in the
SC-SC
setting and Part(b) considers a C-C
case with a stability bound  $\bigO( {C^{K/2}}/{m_1})$ related to the number of iterations and the data size, where $C$ is a constant.
The NC-NC
case is discussed in Part(c) which provides a stability bound of $\bigO(K^{\frac{cl}{cl+1}}/m_1)$, where $c$ is a constant to control the step size and $l$ is the larger smoothness number  of $l_f$ and $l_g$. The conclusions here match the existing results in minmax problems \citep{lei2021stability,farnia2021train}.
\end{remark}

\subsection{Stability Bounds for two-timescale SGD}
Compared with the above SSGD,
Two-timescale SGD (TSGD; Algorithm \ref{alg:TSGD}) always achieves more accurate approximate solutions by updating $\mathbf{y}$ multiple times before updating $\mathbf{x}$. In this section, we extend our analysis from SSGD to TSGD.
Particularly, compared with the results in \citet{bao2021stability}, we provide stability bounds in Theorem \ref{tsgd_stab} for the case where the inner level parameter ($\mathbf{y}$) is updated iteratively ({\em i.e.,} consistency).
We further explore in Theorem \ref{sc_nc} a particular NC-SC setting, which is commonly appeared in bilevel optimization applications such as meta learning and hyperparameter optimization.

\begin{theorem}
\label{tsgd_stab}
    Suppose that Assumptions \ref{asu1} and \ref{asu2} hold and $|g(\cdot,\cdot)| \leq 1$. Let $\mathbf{A}$ be the TSGD  algorithm with $K$ outer-iterations and $T$ inner-iterations. Then we have
    \begin{itemize}
         \item [(a)] Assume that Problem (\ref{pb1}) is SC-SC. Let  $\alpha_x = \alpha_y \leq \frac{1}{l_g}$ be the step sizes. Then,  $\mathbf{A}$ is $\beta$-argument-stable in expectation, where
         $$\beta \leq \bigO \left( {m_{1}}^{-1}\sqrt{L_f^2\alpha_x^2 + \Big(\frac{2T}{\alpha_y(2-\alpha_y l_g)}\Big)^2 }
        \left( 1- \alpha_x \mu_g + \alpha_x l T \right)^{K} \right).$$

        \item [(b)] Assume that Problem (\ref{pb1}) is C-C. Let $\alpha l = \max \{\alpha_x l_f, \alpha_y l_g\}$ and  $\alpha_x$, $\alpha_y \leq \frac{1}{l_g}$ be the step sizes. Then, $\mathbf{A}$ is $\beta$-argument-stable in expectation, where
        $$\beta \leq \bigO \left( m_{1}^{-1}\sqrt{L_f^2\alpha_x^2 + \Big(\frac{2T}{\alpha_y(2-\alpha_y l_g)}\Big)^2 }
        \left(2+2( \alpha l T)^2\right)^{K/ 2}\right). $$

        \item [(c)] Assume that Problem (\ref{pb1}) is NC-NC. Let the step sizes satisfy $\max  \left\{ \alpha_x,\alpha_y  \right\} \leq c/k$ for some constant $c \geq 0$ and $l = \max \left\{ l_f, l_g\right\}$. Then, $\mathbf{A}$ is $\beta$-argument-stable in expectation, where
        $$\beta \leq \bigO \left( ( m_1 T c l )^{-1}\left(2 c L_f \sqrt{l_f^2+T^2 l_g^2}\right)^{\frac{1}{T c l+1}} \cdot K^{\frac{T c l}{T c l+1}} \right).$$
    \end{itemize}

\end{theorem}

\begin{remark}
Compared with the previous results for SSGD, the stability bounds of TSGD depend on the number of iterations in the outer level loop, the number of iterations in the inner level loop, and the data size in the outer level loop. If the step sizes are sufficiently small,  we can see that the bounds in Theorem \ref{tsgd_stab} are asymptotically the same as the bounds of SSGD in Theorem \ref{ssgd_stab}. Thus, Theorem \ref{tsgd_stab} can be considered as a generalization of the previous one.  However, in comparison to Theorem \ref{alg:SSGD}, we can also see that our bounds are exponential in the inner loop iteration $T$. This is caused by the accumulation of the sub-sampled gradients in the inner level, which makes the   outer loop update  overly expansive. The exponential dependence on $T$ also reveals the difference
with the existing work of stability analysis for other problems, such as simple SGD and minmax problems.
\end{remark}

\begin{remark}
\label{remark_bao}
Comparing our results with the ones in
\citep{bao2021stability}, we have the following observations. 1) They only established the uniform stability bound for the Unrolled Differentiation algorithm \ref{alg:UD}, where the algorithm is reinitialized at each time entering the inner level loop, indicating that it takes into account the changes to only one parameter in the outer level loop, while our algorithm considers the update for both parameters.
2) Its proof needs to assume that the update of $y$ in the inner level after the reinitialization will not be affected by  the value specified for $\mathbf{x}$. However, this assumption is quite uncommon and is probably the reason that they do not need to make any assumption on the inner level objective function. 3) The term $\ell(\lambda_t^{\prime}, \hat{\theta}\left(\lambda_t^{\prime}, S^{t r}\right), z)$ in their paper, which corresponds to the inner loss function $g\left( \mathbf{x}_k^{\prime}, \mathbf{y}_k^{\prime}\left( \mathbf{x}_k^{\prime}, D_{m_2}\right), \eta  \right)$ in our paper, includes three arguments ({\em i.e.,} $\lambda_t^{\prime}$, $\hat{\theta}\left(\lambda_t^{\prime}, S^{t r}\right)$ and $z$), where $\hat{\theta}\left(\lambda_t^{\prime}, S^{t r}\right)$ is a function of
$\lambda_t^{\prime}$. This implies that  $\hat{\theta}\left(\lambda_t^{\prime}, S^{t r}\right)$ will change in the inner level updates as a result of the change of $\lambda_t^{\prime}$, and thus cannot be treated as an argument independent of $\lambda_t^{\prime}$, which is misused and causes technical flaws in their proof subsequently.
In contrast, our work allows the inner level parameters to be updated consistently ({\em i.e.,} carrying over the value in the last update), \textit{instead of being reinitialized at each time entering the inner level loop}. Specifically, we allow $y_k^{T}$ to be employed at the beginning of the $(k+1)$-th outer level iteration, rather than $y_0$. This enables us to obtain different stability bounds for different inner level objective functions.
\end{remark}

In the following, we extend our analysis to a particular NC-SC
setting that is frequently encountered in real-world applications and optimization analysis.

\begin{theorem}
\label{sc_nc}
    Suppose that Assumptions \ref{asu1} and \ref{asu2} hold, $0 \leq f(\cdot,\cdot) \leq 1$ and Problem (\ref{pb1}) is NC-SC.  Let $\mathbf{A}$ be the TSGD Algorithm with $K$ outer-iterations and $T$ inner-iterations with $\max  \left\{ \alpha_x,\alpha_y  \right\} \leq c/k$ for constant $c \geq 0$. Denote $l = \max \left\{ l_f, l_g\right\}$. Then, $\mathbf{A}$ is $\beta$-uniform-stable in expectation, where
    $$\beta \leq \bigO \left( \frac{\left(2 c L_f \sqrt{l_f^2+l_g^2 T^2}\right)^{\frac{1}{c (Tl+l-\mu_g)+1}} \cdot K^{\frac{c (Tl+l-\mu_g)}{c (Tl+l-\mu_g)+1}}\left(Tl+l-\mu_g +2/c \right)}{m_1 (Tl+l-\mu_g)} \right).$$

\end{theorem}

\begin{remark}
Compared with our previous analysis, we now sketch the technique differences in our analysis. We consider the bound of the term $(\delta_{\mathbf{x},k},\delta_{\mathbf{y},k}) ^T=(\|\mathbf{x}_k - \mathbf{x}_k^{\prime} \|, \| \mathbf{y}_k - \mathbf{y}_k^{\prime}\|)^T$, while we employ $\delta_k = \sqrt{\left\|\mathbf{x}_k-\mathbf{x}_k^{\prime}\right\|_{2}^{2}+\left\|\mathbf{y}_k-\mathbf{y_k}^{\prime}\right\|_{2}^{2}}$ in the previous analysis, where $(\mathbf{x}_k,\mathbf{y}_k)$, $(\mathbf{x}_k^{\prime},\mathbf{y}_k^{\prime})$ are the outputs of TSGD after $k$ iterations for $D_{m_1}$ and $D_{m_1}^{\prime}$ respectively with $D_{m_1}$ and $D_{m_1}^{\prime}$ differing  in one sample. In the NC-SC
setting, we show that $(\delta_{\mathbf{x},k+1},\delta_{\mathbf{y},k+1}) ^T \leq ((1+\alpha_x l)\delta_{\mathbf{x},k},(1+\alpha_x Tl)\delta_{\mathbf{y},k})^T$ ($\leq$ means the entry-wise inequality), which means our term can be controlled. Then,  we take the expectation of it to derive our uniform stability bound. To achieve the generalization error over continuously changing parameters,
it is imperative to take into account the growth of $(\delta_{\mathbf{x},k},\delta_{\mathbf{y},k})$ instead of $\delta_{\mathbf{x},k}$ in \citep{bao2021stability}. Appendix \ref{thm_tsgd} provides more details.
\end{remark}

Thus, based on our previous results, we now provide the first generalization bounds in the NC-NC setting for both SSGD and TSGD.

\begin{corollary}
\label{ssgd_gen}
    Assume that the problem is NC-NC, $|f(\cdot, \cdot ; \xi)| \leq 1$ for all $\xi$, and  Assumptions \ref{asu1} and \ref{asu2} hold. Denote $l= \max \{l_f,l_g\}$ with $\max \{\alpha_x, \alpha_y\} \leq c/k$ for constant $c \geq 0$. Then, the generalization error of SSGD \ref{alg:SSGD} with $K$ iterations is bounded by
    $\bigO( K^{\frac{cl}{cl+1}}/{m_1} )$.
\end{corollary}

\begin{corollary}
\label{tsgd_gen}
    Assume that the problem is NC-NC, $|f(\cdot, \cdot ; \xi)| \leq 1$ for all $\xi$, and Assumptions \ref{asu1} and \ref{asu2} hold. Let $l= \max \{l_f,l_g\}$ with $\max \{\alpha_x, \alpha_y\} \leq c/k$. Then, the generalization error of TSGD \ref{alg:TSGD} with $K$ outer iterations and $T$ inner iterations is bounded by
    $\bigO( T^{\frac{1}{T c l+1}} K^{1-\frac{1}{T c l+1}}/({m_1}))$.
\end{corollary}

\begin{remark}
By Theorem \ref{thm1_connection}, \ref{ssgd_stab}, \ref{tsgd_stab}], we can derive the above corollaries on generalization error from  stability bounds. Corollary \ref{ssgd_gen} and Corollary \ref{tsgd_gen} show that extremely high number of iterations ($K$ for SSGD and $K$,$T$ for TSGD) will drastically reduce the stability of these algorithms and increase the generalization error, which will make these algorithms increase the risk of overfitting. We will also verify it in the following experiments.
\end{remark}
\section{Experiments}
\label{experiment}
In this section, we empirically validate our previous theoretical results on real world datasets. Two experiments, including meta-learning and hyperparameter optimization, are conducted via Algorithm \ref{alg:TSGD} TSGD (note that when $T=1$, TSGD is just SSGD). Due to the space limitation, we just present the meta learning experiment here, leaving the hyperparameter optimization experiment and other details in the Appendix \ref{additional experiment}.
\subsection{Meta Learning}
\label{exp:meta learning}
Consider the few-shot meta-learning problem with $M$ tasks $\{ \mathcal{T}_i, i=1, ..., M \}$ sampled from distribution $P_{\mathcal{T}}$.
We aim to learn a model that can rapidly adapt to different tasks.
Firstly, the embedding model $\phi$ is shared by all tasks to learn embedded features. Secondly, the task-specific parameter $w_i$ is to adapt the shared embedding to its own sub-problem. Thus, the overall problem of meta-learning can be formulated as follow:
\begin{subequations}
\begin{align}
  &\min _\phi \mathcal{L}_{\mathcal{D}}\left(\phi, \bar{w}^*\right)=\mathbb{E}_{\xi \in \mathcal{D}_i^{\mathrm{te}},\mathcal{T}_i }\left[\mathcal{L}\left(\phi, w_i^* ; \xi\right) \right], \label{eq:meta_learner}\\
   &\text { s.t. } \bar{w}^*=\underset{\bar{w}}{\arg \min } \Big[\mathcal{L}_{\mathcal{D}^{\mathrm{tr}}}(\phi, \bar{w})=\mathbb{E}_{\mathcal{T}_i }\left[ \mathcal{L}_{\mathcal{D}_i^{\mathrm{tr}}}\left(\phi, w_i\right)\Big]\right]. \label{eq:base_learner}
\end{align}
\end{subequations}
\vspace{-0.15in}

where $\mathcal{D}^{\mathrm{tr}}_i$ and $\mathcal{D}^{\mathrm{te}}_i$ are the training and testing datasets for task $\mathcal{T}_i$. Each $w_i$ is computed from one or more gradient descent updates from $\bar{w}$ on the corresponding task (rapid adaptation), {\em i.e.,} $w_i=\bar{w}-\alpha \nabla_{\bar{w}} \mathcal{L}_{\mathcal{D}_{tr}}(\phi, w_i)$.
In the inner level, the base learner optimizes the series of ${w_i}$ for each tasks (Equation \ref{eq:base_learner}). In the outer level, the meta-learner optimizes the embedding model $\phi$ using the minimizers $w^*_i$ learned from the inner level and computes the loss from the testing dataset (Equation \ref{eq:meta_learner}).

\noindent {\bf Settings and Implementation}
We evaluate the behavior of the 5-way-1-shot task on the Omnilot dataset \citep{lake2015human}, {\em i.e.,} it aims to classify 5 unseen classes from only 1 labeled sample. It contains 1623 different handwritten characters from 50 different alphabets. The image is in greyscale with a size $28\times28$. We follow similar settings in \cite{ji2021bilevel}. A five-layer fully-connected network is constructed, where the task-specific parameter $w_i$ corresponds to the last layer of the network and the shared embedding model $\phi$ corresponds to all preceding layers. Thus, we train the two sets of layers separately in the outer and inner level of optimization. We build our model and establish our training using the software library {\fontfamily{cmss}\selectfont learn2learn} \citep{Arnold2020}. We follow the official train-validation-test partition and train $\phi$, $w_i$ using the training set. The size of each layer in the network is $784 \rightarrow 256 \rightarrow 128 \rightarrow 64 \rightarrow 64 \rightarrow 5$. We set the number of tasks for training and testing set to 2000 and the batch size of tasks to 32. The learning rate of $\phi$ and $w_i$ are 0.002 and 0.01, respectively. Results are evaluated based on the average of 5 trial runs with different random seeds.

\noindent {\bf Results Evaluation}
Figure \ref{fig:metalearning} presents the learning curves on training set, testing set and the generalization error with different values of inner iterations $T$ and outer iterations $K$. Generalization error is estimated by the difference between training and testing loss. On one hand, it can be seen that the model easily overfits on the testing set as $K$ increases drastically (Figure \ref{fig:metatestloss}) and the effect of $T$ is very limited. On the other hand, with an appropriate value of $K$, smaller $T$ (\textit{i}.\textit{e}) will result in underfitting on the testing loss ($T=1$ in the Figure \ref{fig:metagenloss} causes highest generalization error due to the underfitting training process). The trend of generalization error in terms of $K$ and $T$ matches with our analysis in Theorem \ref{sc_nc}.

\begin{figure}[]
     \centering
     \begin{subfigure}[1]{0.32\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/meta_train_loss_2500.pdf}
         \caption{Meta-training Loss}
         \label{fig:metatrainloss}
     \end{subfigure}
     \hfill
     \begin{subfigure}[2]{0.32\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/meta_test_loss_2500.pdf}
         \caption{Meta-testing Loss}
         \label{fig:metatestloss}
     \end{subfigure}
     \hfill
     \begin{subfigure}[3]{0.32\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/meta_gap_loss_2500.pdf}
         \caption{Generalization Error}
         \label{fig:metagenloss}
     \end{subfigure}
        \caption{Results of meta learning with various values of $T$ and $K$}
        \label{fig:metalearning}
\end{figure}

\section{Conclusion}
\label{conclusion}
We give a thorough analysis on the generalization of first-order (gradient-based) methods for the bilevel optimization framework. In particular, we establish a quantitative connection between generalization and algorithmic stability and provide the first generalization bounds of the continuous updates for inner parameters and outer parameters in multiple settings. Our experiments suggest that inappropriate iterations will cause underfitting and overfitting easily. The tendency of generalization error also validates our theoretical results.

From the discussion in previous sections, we only discussed the first-order method, while there exists a number of estimating second-order and momentum-based approaches to solve the bilevel optimization problem. Dealing with the approximation of hypergradient in generalization analysis is another direction for future work.

\section{Introduction}
\label{Introduction}
(Stochastic) bilevel optimization is a widely confronted
problem in machine learning with various applications such as
meta-learning \citep{finn2017model, bertinetto2018meta,rajeswaran2019meta}, hyper-parameter optimization \citep{franceschi2018bilevel,shaban2019truncated,baydin2017online,bergstra2011algorithms,luketina2016scalable}, reinforcement learning \citep{hong2020two}, and few-shot learning \citep{koch2015siamese,santoro2016meta,vinyals2016matching}.
The basic form of this problem can be  defined as follows
\begin{equation}\label{pb1}
  \begin{aligned}
  \min _{\mathbf{x} \in \mathbb{R}^{d_{1 }}}  R(\mathbf{x})&= F\left(\mathbf{x}, \mathbf{y}^{*}(\mathbf{x})\right):=\mathbb{E}_{\xi}\left[f\left(\mathbf{x}, \mathbf{y}^{*}(\mathbf{x}) ; \xi\right)\right] \\
  \text { s.t. }  \mathbf{y}^{*}(\mathbf{x})&=\arg \min _{\mathbf{y} \in \mathbb{R}^{d_{2}}}\left\{G(\mathbf{x}, \mathbf{y}):=\mathbb{E}_{\zeta}[g(\mathbf{x}, \mathbf{y} ; \zeta)]\right\},
  \end{aligned}
  \end{equation}
where $f$ : $\mathbb{R}^{d_1} \times \mathbb{R}^{d_{2 }} \rightarrow \mathbb{R} $ and $g$ : $\mathbb{R}^{d_1} \times \mathbb{R}^{d_{2 }} \rightarrow \mathbb{R} $ are two continuously differentiable loss functions with respect to $\mathbf{x}$ and  $\mathbf{y}$. Problem (\ref{pb1}) has an optimization hierarchy of two levels, where the outer-level objective function $f$ depends on the minimizer of the inner-level objective function $g$.

 Due to its importance,

 the above bilevel optimization problem has received considerable attention in recent years. A natural way to solve problem (\ref{pb1}) is to apply alternating stochastic gradient updates with approximating $\nabla_{\mathbf{y}} g(\mathbf{x}, \mathbf{y})$ and $\nabla f(\mathbf{x} , \mathbf{y})$, respectively. Briefly speaking, previous efforts mainly examined two types of methods to perceive an approximate solution that is close to the optimum $y^* (\mathbf{x})$. One is to utilize the single-timescale strategy \citep{chen2021single, guo2021randomized, khanduri2021near,hu2022improved}, where the updates for $\mathbf{y}$ and $\mathbf{x}$ are carried out simultaneously. The other one is to apply the two-timescale strategy \citep{ghadimi2018approximation, ji2021bilevel,hong2020two,pedregosa2016hyperparameter}, where the  update of $\mathbf{y}$ is repeated multiple times to achieve a more accurate approximation before conducting the  update of $\mathbf{x}$.

 While there is a long list of work on bilevel optimization,
most of the existing work only focuses on either analyzing its convergence behaviors \citep{ghadimi2018approximation, hong2020two, ji2021bilevel} or improving its convergence rate, based on the convexity and the smoothness properties of  $f(\cdot, \cdot)$ and/or $g(\cdot,\cdot)$ \citep{liu2020generic, li2020improved}.

 Contrarily,
 only little effort is devoted to

 understanding the generalization behavior of the problem. To the best of our knowledge, only one recent work considers
 generalization

 \citep{bao2021stability}, which presents the first expected stability bound

 for the bilevel optimization problem.  Unfortunately, the problem is still far from well-understood and   there are several undesirable issues in this work: (1) Their study is only for the hyper-parameter optimization problem, which is only a special instance of bilevel optimization. Thus, the stability for general bilevel optimization is still unknown; (2) Their result is only for the uniform stability (which could be deduced from argument stability with certain conditions, see Definition \ref{def:stab} for details), leaving the analysis of other stronger definitions of algorithmic stability open; (3) \cite{bao2021stability} only considers a general setting for inner function. However, there are some technical flaws in their analysis (see Remark \ref{remark_bao} for detailed discussions) and
 their bound is quite loose in some cases. For example, in the case where the inner function is strongly convex, it is known that even in the classical stochastic convex optimization, strongly convex functions have better stability than the convex ones \citep{hardt2016train}. Thus, it is possible to provide improved stability bounds in some cases via different and novel analysis; (4) Their work only considers a special class of optimization algorithms, which allows the outer level parameters to be updated continuously but needs to reinitialize the inner level parameters before each iteration in the inner level. However, such algorithms  are not commonly used in practice due to their inefficiency;  (5) Their experiments take only hyper-parameter optimization into consideration and neglect other applications in the bilevel optimization instances.

To address all the aforementioned issues, we give in this paper a thorough analysis on the generalization behaviors of first-order (gradient-based) methods for general bilevel optimization problem. We
employ the recent advances of algorithmic stability to investigate the generalization behaviors
in different settings. Specifically,  our main contributions can be summarized as follows:
\begin{itemize}
    \item Firstly, we establish a fundamental connection between generalization error and different notations  of algorithmic stability (argument stability and uniform stability) for
     bilevel optimization problem in both  expectation and high probability forms. Specifically, we

    show that the high probability form of the generalization error bound can be improved from

   $\bigO(\sqrt{n})$ to $\bigO(\log n)$ compared with the result in \citet{bao2021stability}.

    \item Next, we present the first stability bounds for gradient-based methods with either single-timescale or two-timescale update strategy under different standard settings. In detail, we consider the settings of
    strongly-convex-strongly-convex (SC-SC), convex-convex (C-C), and nonconvex-nonconvex (NC-NC), and further extend our analysis to a particular nonconvex-strongly-convex (NC-SC) setting that is widely appeared in practice. Especially, for the SC-SC and C-C cases with single-timescale update strategy we significantly improve the generalization bounds compared with \citep{bao2021stability} (see Table \ref{tab:main results} for details).  To the best of our knowledge, our work is the first thorough generalization analysis for general  bilevel optimization problem.
    \item Thirdly, we provide the first generalization bounds for the case where

    both the outer and inner level parameters are subject to continuous (iterative) changes. Compared to the previous work \citep{bao2021stability}, our work does not need the reinitialization step before each iteration in the inner level

    and hence our algorithm can carry over the last updated inner level parameters, which is more general and practical.
    \item Finally, we conduct empirical studies to corroborate our theories
    via meta-learning and hyperparameter optimization, which are two  applications of bilevel optimization.
\end{itemize}
Due to space limitations, all the proofs and additional experiments are included in Appendix.
\begin{table*}[t]
\captionsetup{font=scriptsize}
\centering
\begin{center}
\footnotesize
\scriptsize
\begin{sc}
\scalebox{0.77}{
\begin{tabular}{ccccccccc}
\toprule
\multirow{2}{*}{Reference} & \multicolumn{4}{c}{stability bounds in various settings} & \multicolumn{2}{c}{stability}                         & \multicolumn{2}{c}{consistency}         \\ \cmidrule(lr){2-5}  \cmidrule(lr){6-7}  \cmidrule(lr){8-9}
                           & SC-SC         & C-C        & NC-NC        & NC-SC        & argument                  & uniform                   & Outer            & inner                \\ \hline
UD \citep{bao2021stability}                         & \multicolumn{4}{c}{$\bigO(K^{\kappa_1}/m_1)$} & & \checkmark & \checkmark &                           \\
SSGD (This work)                  & $\bigO(1/m_1)$ & $\bigO(C^{K/2} /m_1)$ & $\bigO(K^{\kappa_2}/m_1)$ & $\bigO(K^{\kappa_3}/m_1)$ &\checkmark &       \checkmark                    & \checkmark & \checkmark \\
TSGD (This work)               & $\bigO((\kappa_4 T)^K/m_1)$ & $\bigO((T)^K/m_1)$ & $\bigO( T^{1-\kappa_4} K^{\kappa_5}/m_1)$ & $\bigO( T^{1-\kappa_6} K^{\kappa_6}/m_1)$           & \checkmark & \checkmark & \checkmark & \checkmark \\
\bottomrule
\end{tabular}}
\caption{
Comparison of related work with ours.
$\kappa_i$: a constant such that $\kappa_i\in(0,1)$ for all $i$ above; $T$: inner iterations; $K$: outer iterations; $m_1$: size of outer dataset. UD indicates the unrolled differentiation (Algorithm \ref{alg:UD}). SSGD and TSGD stand for Algorithm \ref{alg:SSGD} and Algorithm \ref{alg:TSGD}, the single-timescale and two-timescale methods, via stochastic gradient descent. We provide a comparison of TSGD and UD in the Appendix \ref{comparision between ud and tsgd}. Consistency indicates whether it is allowed to update inner-level or outer-level parameters continuously.
}
\label{tab:main results}
\end{sc}
\end{center}
\vspace{-0.5cm}
\end{table*}

\subsection{Related work}
\label{related_work}

Research at the interface between generalization and the bilevel problem can be roughly classified into two categories. The first one includes all the research on bilevel optimization. In recent decades, extensive studies have been done on this topic, which suggests that bilevel optimization has a wide range of applications in machine learning such as
hyper-parameter optimization \citep{franceschi2018bilevel,lorraine2018stochastic,okuno2021lp}, meta learning \citep{bertinetto2018meta,rajeswaran2019meta,soh2020meta} and reinforcement learning \citep{yang2018convergent,tschiatschek2019learner}. Most of the existing work studies the problem  from an optimization perspective. For example,
\citet{ghadimi2018approximation,ji2021bilevel} provide the convergence rate analysis based on the nonconvex-strongly-convex assumption for the two functions $f(\cdot, \cdot)$ and $g(\cdot,\cdot)$. \citep{grazzi2020iteration} considers the iteration complexity for hypergradient computation. \citep{liu2020generic, li2020improved} present an asymptotic analysis for the convex-strongly-convex setting. Perhaps the most related one to ours
from the generalization standpoint ({\em i.e.,} the expectation of population risk and empirical risk) is \citet{bao2021stability}. \textbf{Table \ref{tab:main results}} is the comparison between our work and  the results in \citep{bao2021stability}.  We can see that
our analysis is more comprehensive as we consider various settings. Firstly, in this paper we consider the general bilevel problem defined in Problem (\ref{pb1}), while \citet{bao2021stability} only concentrates  on the hyper-parameter optimization, which is just an instance of  bilevel optimization problem (to be detailed in Section \ref{experiment} and Appendix \ref{additional experiment}).
Secondly, our work adopts a stronger version of stability called {\em argument stability}, which implies the previously used {\em uniform stability} if the function is sufficiently smooth. Furthermore, our work does not need to reinitialize the inner-level parameters and allows them to carry over their last updated parameters at each time updating the inner level. This indicates that
$\mathbf{y}$ in the inner level is updated iteratively  and depends on the current parameter of $\mathbf{x}$, which is more common and efficient in practice.

The second category includes all the work on stability analysis. There is a long list of research on
stability and generalization \citep{bousquet2002stability, mukherjee2006learning, shalev2010learnability}. \citet{bousquet2002stability} first introduces the notion of \textit{uniform stability} and establishes the first framework of stability analysis. \citet{hardt2016train} later extends the stability analysis  to iterative algorithms based on stochastic gradient methods for the vanilla stochastic optimization. After that, there are  subsequent studies on generalization analysis for various problems via algorithmic stability, such as minmax problems \citep{lei2021stability, farnia2021train,zhang2021generalization} and pairwise learning \citep{yang2021simple, lei2020sharper, xue2021differentially,huai2020pairwise}. However, it is notable that due to the additional stochastic function in the constraint in the bilevel optimization, all the previous techniques and results cannot be applied to our problem. Although the generalization analysis of minmax optimization is somewhat similar to ours, it involves only one objective function $f$ and a single level in algorithms for typical minmax optimization problems, while in the bilevel optimization algorithms there is an inner level and an outer level, which is considerably more challenging.

\section{Preliminaries}
\label{preliminaries}
\subsection{Definitions and Assumptions}
In the following, we give some necessary definitions and assumptions that are widely used in bilevel optimization \citep{ghadimi2018approximation,ji2021bilevel,khanduri2021near} and generalization analysis \citep{hardt2016train,lei2021stability}.
\begin{definition}[Joint Lipschitz Continuity]
    A function $f(\mathbf{x},\mathbf{y})$ is jointly $L$-Lipschitz over $\mathbb{R}^{d_1} \times \mathbb{R}^{d_2}$, if for all $\mathbf{x} \in \mathbb{R}^{d_1}, \mathbf{y} \in \mathbb{R}^{d_2}$, the following holds,   $|f(\mathbf{x},\mathbf{y})-f(\mathbf{x}^{\prime},\mathbf{y}^{\prime})| \leq L\sqrt{\left\|\mathbf{x}-\mathbf{x}^{\prime}\right\|_{2}^{2}+\left\|\mathbf{y}-\mathbf{y}^{\prime}\right\|_{2}^{2}}$.
\end{definition}

\begin{definition}[Smoothness]
     A function $f$ is $l$-smooth over a set $S$ if for all $u, w \in S$ the following is true, $\|\nabla f(u)-\nabla f(w)\| \leq l\|u-w\|$.
\end{definition}

\begin{definition}[Strong Convexity]
     A function $f$ is $\mu$-strongly-convex over a set $S$, if for all $u, w \in S$, the following holds, $f(u)+\langle\nabla f(u), w-u\rangle+\frac{\mu}{2}\|w-u\|^{2} \leq f(w)$.
\end{definition}

\begin{assumption}[Inner-level Function Assumption]\label{asu1}
  We assume the inner stochastic function  $g(\mathbf{x},\mathbf{y})$ in (\ref{pb1}) satisfies the following:\\
    \subasu $g(\mathbf{x},\mathbf{y})$ is jointly $L_g$-Lipschitz for any $\mathbf{x}\in \mathbb{R}^{d_1}$ and $\mathbf{y} \in \mathbb{R}^{d_2}$.\\
    \subasu $g(\mathbf{x},\mathbf{y})$ is continuously differentiable and  $l_g$-smooth for any $(\mathbf{x},\mathbf{y}) \in \mathbb{R}^{d_1} \times \mathbb{R}^{d_2}$.

\end{assumption}

\begin{assumption}[Outer-level Function Assumption]\label{asu2}
  We assume the outer stochastic function  $f(\mathbf{x},\mathbf{y})$ in (\ref{pb1}) satisfies the following:\\
    \subasu $f(\mathbf{x},\mathbf{y})$ is jointly $L_f$-Lipschitz for any $ \mathbf{x} \in \mathbb{R}^{d_1}$ and $\mathbf{y} \in \mathbb{R}^{d_2}$.\\
    \subasu $f(\mathbf{x},\mathbf{y})$ is continuously differentiable and $l_f$-smooth for any $(\mathbf{x},\mathbf{y}) \in \mathbb{R}^{d_1} \times \mathbb{R}^{d_2}$.

\end{assumption}

\subsection{Problem Formulation}\label{sec:problem}
Given two distributions $\mathbb{D}_1$ and $\mathbb{D}_2$, in the (stochastic) optimization problem we aim to find the minimizer of Problem (\ref{pb1}). However, since the distributions are often unknown, in practice we only have two finite-size datasets  $D_{m_1}= \{\xi_i\mid i=1,...,m_1\}\sim \mathbb{D}_1^{m_1}$ and $D_{m_2}= \{\zeta_i\mid i=1,...,m_2\} \sim \mathbb{D}_2^{m_2}$, where each $\xi_i$ and $\zeta_i$ are i.i.d. sampled from $\mathbb{D}_1$ and $\mathbb{D}_2$, respectively. Based on these datasets, we will design some (randomized) algorithm ${\bf A}$ with output ${\bf A}(D_{m_1}, D_{m_2})= (\mathbf{x}, \mathbf{y})\in \mathbb{R}^{d_1}\times \mathbb{R}^{d_2}$. Our goal is to investigate the generalization behavior of such output. Note that although there are two stochastic functions in the bilevel optimization problem, we only care about the generalization of the outer-level one since it is the one that we prefer to minimize.

Below we define the generalization error to measure the generalization behavior.
Given  distribution $\mathbb{D}_1$  and  a finite data $D_{m_1}\sim \mathbb{D}_1^{m_1} $, the population risk function $R(\mathbf{x},\mathbf{y},\mathbb{D}_1)$
of $\mathbf{x}$, $\mathbf{y}$ on $\mathbb{D}_1$ is defined as $R(\mathbf{x},\mathbf{y},\mathbb{D}_1) := \mathbb{E}_{{\xi} \sim \mathbb{D}_1}\left[f\left(\mathbf{x}, \mathbf{y}(\mathbf{x}) ; \xi\right)\right]$, and its empirical risk function  on $D_{m_1}$ is
 $R_s(\mathbf{x},\mathbf{y}, {D}_{m_1}) = \frac{1}{m_1} \sum_{i=1}^{m_1}\left[f\left(\mathbf{x}, \mathbf{y}(\mathbf{x}) ; \xi_i\right)\right]$. Moreover, for a fixed hyperparameter $\mathbf{x} \in \mathbb{R}^{d_{1 }} $ and  $\mathbf{y}(\mathbf{x}) \in \mathbb{R}^{d_{2}}$ ( note that $\mathbf{y}(\mathbf{x})$ might be dependent on $\mathbf{x}$),
  we define the difference between the population risk and the empirical risk over $(\mathbf{x}, \mathbf{y}(\mathbf{x}))$ as the {\em bilevel generalization error of $(\mathbf{x}, \mathbf{y}(\mathbf{x}))$}:
  $\mathbb{E}_s[R(\mathbf{x},\mathbf{y})-R_s(\mathbf{x},\mathbf{y})]$, where $\mathbb{E}_s$ denotes the expectation of $D_{m_1}\sim \mathbb{D}_1^{m_1}$. When there is no ambiguity, we simplify thereafter the notations as follows:  $R(\mathbf{x},\mathbf{y},\mathbb{D}_1) =R(\mathbf{x},\mathbf{y}) $ and $R_s(\mathbf{x},\mathbf{y},D_{m_1}) = R_s(\mathbf{x},\mathbf{y})$.

  Our goal is thus to analyze the {\em bilevel generalization error of the output of algorithm ${\bf A}(D_{m_1}, D_{m_2})$} based on $D_{m_1}$ and  $D_{m_2}$. Since the generalized error depends on the algorithm itself, in the following we will introduce the algorithms to be considered in this paper.

Most of the existing algorithms adopt the following idea:
first approximate $\mathbf{y}^*$ on $D_{m_2}$ for a given parameter $\mathbf{x}$ in the inner level and then seek the hyperparameter $\mathbf{x}^*(D_{m_1},D_{m_2})$ with corresponding hypothesis $\mathbf{y}^*(\mathbf{x}^*(D_{m_1},D_{m_2}),D_{m_2})$ by the below estimation:
\begin{equation}
\label{eq:approximation}
\begin{aligned}
     \hat{\mathbf{x}}(D_{m_1},D_{m_2}) & \approx \arg \min_{\mathbf{x}} R_s(\mathbf{x},\hat{\mathbf{y}}(\mathbf{x},D_{m_2}),D_{m_1}),\\
      \text{where} \quad \hat{\mathbf{y}}(\mathbf{x},D_{m_2}) & \approx \arg \min_{\mathbf{y}}G_s(\mathbf{x},\mathbf{y},D_{m_2}),
\end{aligned}
\end{equation}
      where $G_s(\mathbf{x},\mathbf{y},D_{m_2})$ is the empirical risk of $G(\mathbf{x},\mathbf{y})$ over $D_{m_2}$, {\em i.e.,}  $G(\mathbf{x},\mathbf{y},D_{m_2}) = \frac{1}{m_2} \sum_{i=1}^{m_2}g\left(\mathbf{x}, \mathbf{y}(\mathbf{x}) ; \zeta_i\right)$. Most of the current gradient-based (first-order) algorithms for approximating (\ref{eq:approximation}) can be categorized into two classes:  single-timescale methods and two-timescale methods. The single-timescale method performs the updates for $\mathbf{y}$ and $\mathbf{x}$ simultaneously via stochastic gradient descent (SGD), while the two-timescale method updates $\mathbf{y}$ multiple times  before updating $\mathbf{x}$ (via stochastic gradient descent). As there are numerous approaches for both classes (see Related Work section for details), in this paper we will analyze the generalization behaviors for the most classical and standard one in each class, {\em i.e.,} single-timescale SGD (SSGD; Algorithm \ref{alg:SSGD}) and two-timescale SGD (TSGD; Algorithm \ref{alg:TSGD}). There is a long list of work \citep{chen2021single}, \citep{ghadimi2018approximation,ji2021bilevel}

      based on either SSGD or TSGD.

\begin{figure}[]
\begin{minipage}{0.5\textwidth}
\vspace{-33pt}
\begin{algorithm}[H]
  \begin{algorithmic}[1]
  \caption{Single-timescale SGD (SSGD)}
  \label{alg:SSGD}
    \STATE {\bfseries Input:} number of iterations $K$, step sizes $\alpha_x$, $\alpha_y$, initialization $\mathbf{x}_0$,$\mathbf{y}_0$, Datasets $D_{m_1}$ and $D_{m_2}$
    \STATE {\bfseries Output:} $\mathbf{x}_K$, $\mathbf{y}_K$
    \FOR {$k=0$ {\bfseries to} $K-1$}
     \STATE Uniformly sample $i\in [m_2]$, $j\in [m_1]$
      \STATE $\mathbf{y}_{k+1} = \mathbf{y}_k - \alpha_y \nabla_{\mathbf{y}} g(\mathbf{x}_k, \mathbf{y}_k; \zeta_i) $
      \STATE $\mathbf{x}_{k+1} =\mathbf{x}_k -\alpha_{x} \nabla f(\mathbf{x}_k, \mathbf{y}_k;\xi_j) $
    \ENDFOR
    \STATE {\bf return} $\mathbf{x}_K$ and $\mathbf{y}_K$
  \end{algorithmic}
\end{algorithm}
\end{minipage}
\begin{minipage}{0.50\textwidth}
\begin{algorithm}[H]
  \begin{algorithmic}[1]
  \caption{Two-timescale SGD (TSGD) }
  \label{alg:TSGD}
    \STATE {\bfseries Input:} number of iterations $K$, step sizes $\alpha_x$, $\alpha_y$, initialization $\mathbf{x}_0$,$\mathbf{y}_0$
    \STATE {\bfseries Output:} $\mathbf{x}_K$, $\mathbf{y}_K$
    \FOR {$k=0$ {\bfseries to} $K-1$}
      \STATE $\mathbf{y}_{k}^{0} \gets \mathbf{y}_{k-1}^{T}$
      \FOR{$t=0$ {\bfseries to} $T-1$}
      \STATE Uniformly sample $i\in [m_2]$
      \STATE 	$\mathbf{y}_{k}^{t+1} = \mathbf{y}_k^{t} - \alpha_y \nabla_{\mathbf{y}} g(\mathbf{x}_k^t, \mathbf{y}_k^t; \zeta_i) $
      \ENDFOR
      \STATE Uniformly sample $j\in [m_1]$
      \STATE $\mathbf{x}_{k+1} =\mathbf{x}_k -\alpha_{x} \nabla f(\mathbf{x}_k, \mathbf{y}_k^T;\xi_j) $
    \ENDFOR
    \STATE {\bf return} $\mathbf{x}_K$, $\mathbf{y}_K^{T}$
  \end{algorithmic}
\end{algorithm}
\end{minipage}
\end{figure}

\section{Generalization and Stability for Bilevel Optimization}
\label{generalization and stability}
Algorithmic stability is one of the classical approaches to analyzing  the generalization bound for algorithms. Roughly speaking, the algorithmic stability of (randomized) algorithm ${\bf A}$ measures how the output of algorithm  ${\bf A}$ changes if we change one data sample in the input dataset.
While there are various notions of stability, most of the existing work on analyzing the stability of stochastic optimization, pairwise learning and minimax optimization focuses on the  uniform-stability \citep{bousquet2002stability} and the argument-stability \citep{liu2017algorithmic,lei2020fine}. Thus, we also adopt these two notions of stability for the bilevel optimization problem. Briefly speaking, uniform-stability focuses on the resulting change in population risk function, while the argument-stability considers the resulting change in arguments, {\em i.e.,} the output of the algorithm.
\begin{definition}[Algorithmic Stability]\label{def:stab}
    Let $\mathbf{A}: \mathbb{D}_1^{m_1} \times \mathbb{D}_2^{m_2} \mapsto \mathbb{R}^{d_1} \times \mathbb{R}^{d_2}$ be a randomized algorithm.
    \begin{itemize}
    \item[(a)] $\mathbf{A}$ is $\beta$-uniformly-stable if for all datasets $D_{m_1}, D^{\prime}_{m_1} \sim \mathbb{D}_1^{m_1}$ and $ D_{m_2} \sim \mathbb{D}_2^{m_2}$ such that $D_{m_1}$ and $D^{\prime}_{m_1}$ differ in at most one sample, we have the following for any $\xi \sim \mathbb{D}_1$:
    \begin{equation*}
        \mathbb{E}_{\mathbf{A}}[|f(\mathbf{A}(D_{m_1},D_{m_2}),\xi) - f(\mathbf{A}(D^{\prime}_{m_1},D_{m_2}),\xi) |] \leq \beta.
    \end{equation*}
    $\mathbf{A}$ is $\beta$-uniformly-stable with probability at least $1-\delta$ if we have the following for any $\xi \sim \mathbb{D}_1$ with probability at least $1-\delta$:
     \begin{equation*}
         \left| f(\mathbf{A}(D_{m_1},D_{m_2}),\xi) - f(\mathbf{A}(D^{\prime}_{m_1},D_{m_2}),\xi) \right|  \leq \beta.
    \end{equation*}
    \item[(b)] $\mathbf{A}$ is $\beta$-argument-stable in expectation if for all datasets $D_{m_1}, D^{\prime}_{m_1} \sim \mathbb{D}_1^{m_1}$ and $ D_{m_2} \sim \mathbb{D}_2^{m_2}$ such that $D_{m_1}$ and $D^{\prime}_{m_1}$ differ in at most one sample, we have:
    \begin{equation*}
         \mathbb{E}_{\mathbf{A}}[\|\mathbf{A}(D_{m_1},D_{m_2}) - \mathbf{A}(D^{\prime}_{m_1},D_{m_2}) \|_{2} ] \leq \beta.
    \end{equation*}
    \end{itemize}
\end{definition}
Note that the definition of uniform stability in expectation is the same as the definition in \citep{bao2021stability}. Thus, our other definitions can be considered as extensions of the previous stability for bilevel optimization. In the following, we present Theorem \ref{thm1_connection} as our first result, which shows a crucial relationship between generalization error and algorithmic stability for an algorithm ${\bf A}$.

\begin{theorem}\label{thm1_connection}
    Let $\mathbf{A}: \xi^{m_1} \times \zeta^{m_2} \mapsto \mathbb{R}^{d_1} \times \mathbb{R}^{d_2}$ be a randomized BO algorithm.
    \begin{itemize}
    \item[(a)] If $\mathbf{A}$ is $\beta$-uniform-stable in expectation, then the following holds for $D_{m_1} \sim \mathbb{D}_1^{m_1}$, $D_{m_2} \sim \mathbb{D}_2^{m_2}$:

    \begin{equation*}
        \mathbb{E}_{\mathbf{A},D_{m_1}}[R(\mathbf{A}(D_{m_1},D_{m_2}))-R_{s}(\mathbf{A}(D_{m_1},D_{m_2}))] \leq \beta.
    \end{equation*}

    \item[(b)] If $\mathbf{A}$ is $\beta$-argument-stable in expectation and Assumption \ref{asu2} holds, then the following holds for $D_{m_1} \sim \mathbb{D}_1^{m_1}$, $D_{m_2} \sim \mathbb{D}_2^{m_2}$:
    \begin{equation*}
        \mathbb{E}_{\mathbf{A},D_{m_1}}[R(\mathbf{A}(D_{m_1},D_{m_2}))-R_{s}(\mathbf{A}(D_{m_1},D_{m_2}))] \leq L_f\beta.
    \end{equation*}

    \item[(c)] Assume that $\left|f(\mathbf{x},\mathbf{y};\xi)\right| \leq M$ for some $M \geq 0$. If $\mathbf{A}$ is $\beta$-uniform-stable almost surely, then for $D_{m_1} \sim \mathbb{D}_1^{m_1}$, $D_{m_2} \sim \mathbb{D}_2^{m_2}$, the following holds with probability $1-\delta$:
    \begin{equation*}
        \begin{aligned}
            &\left|R(\mathbf{A}(D_{m_1},D_{m_2}))-R_{s}(\mathbf{A}(D_{m_1},D_{m_2})) \right|\\
            & \leq 2\beta+e\left(\frac{4 M}{\sqrt{m_{1}}} \sqrt{\log \frac{e}{\delta}}+12 \sqrt{2} \beta \lceil  \log _{2} m_1  \rceil \sqrt{\log \frac{e}{\delta}}\right).
        \end{aligned}
    \end{equation*}

    \end{itemize}
\end{theorem}
\begin{remark}
The above theorem suggests that the generalization error can be controlled by several notions of algorithmic stability. Part (a) and Part (b) show that the expectation of generalization error can be bounded by uniform stability and argument stability with the Lipschitz constant, respectively; Part (c) indicates that the generalization error for the algorithm is no more than
$\bigO(\beta\log(m_1) +1/\sqrt{m_1})$ with probability $1-\delta$. Compared with the existing work \citep{bao2021stability}, Theorem \ref{thm1_connection} considers argument stability additionally, which is a stronger notion of stability than uniform stability (since uniform stability can be deduced from argument stability with the condition that the function is sufficiently smooth).
Moreover, we use the McDiarmid's inequality and the equivalence of tails and moments for the random variable with a mixture of sub-gaussian and sub-exponential tails (Lemma 1 in \citet{bousquet2020sharper}), which provide a significantly improved high probability bound   in Part (c) ({\em i.e.,} improving from $\bigO(\beta\sqrt{m_1})$ in \citet{bao2021stability} to $\bigO(\beta \log m_1)$).
\end{remark}

\section{Stability Analysis for Bilevel Optimization Algorithms }
\label{main results}
Motivated by Theorem \ref{thm1_connection}, we can see that to analyze the generalization behaviors for any algorithm, it is sufficient to analyze its stability. As mentioned in the previous Section \ref{sec:problem}, we will consider the stability of SSGD and TSGD. For simplicity we let  SC-SC denote the case where $f$ and $g$ both  are strongly convex functions. C-C, NC-NC, and NC-SC are also denoted in a similar manner with ''C'' representing convex function and ''NC'' representing nonconvex function.

\subsection{Stability Bounds for Single-timescale SGD}
As we can see from Algorithm \ref{alg:SSGD}, SSGD updates $\mathbf{y}$ and $\mathbf{x}$ simultaneously. In the following we develop stability bounds for this algorithm in different settings.

\begin{theorem}\label{ssgd_stab}
    Suppose that Assumptions \ref{asu1} and \ref{asu2} hold and Algorithm $\mathbf{A}$ is SSGD with $K$ iterations:
    \begin{itemize}
        \item [(a)] Assume that Problem (\ref{pb1}) is SC-SC with strongly convexity parameters $\mu_f$ and $\mu_g$. Let $\alpha_x = \alpha_y$ (see  Lemma \ref{lemma:expansive} for details) be the step sizes. Denote $l = \max\{ l_f, l_g\}$. Then, $\mathbf{A}$ is $\beta$-argument-stable in expectation, where $$ \beta \leq \bigO\left(  \big(L_f^{2}+L_g^{2}\big)^{\frac{1}{2}}\left(m_1 \left(\mu_f+\mu_{g}- (\alpha_{x} l)^2/2 + 0.25 \right) \right)^{-1} \right).$$

        \item [(b)] Assume that Problem (\ref{pb1}) is C-C. Let $\alpha_x$, $\alpha_y$ be the step sizes. Then,  $\mathbf{A}$ is $\beta$-argument-stable in expectation, where
        $$\beta \leq \bigO \left( m_1^{-1}\sqrt{\left(\alpha_{x} L_f \right)^{2}+\left(\alpha_{y} L_g\right)^{2}}\left(2+2\max \left\{\left(\alpha_{x} l_f\right)^{2}, \left(\alpha_{y} l_g\right)^{2} \right\} \right)^{K/ 2} \right). $$

        \item [(c)] Assume that Problem (\ref{pb1}) is NC-NC. Let the step sizes satisfy $\max  \left\{ \alpha_x,\alpha_y  \right\} \leq c/k$ for some constant $c \geq 0$ and $l = \max \left\{ l_f, l_g\right\}$. Then, $\mathbf{A}$ is $\beta$-argument-stable in expectation, where
        $$\beta \leq \bigO \Big( (m_1  c l)^{-1}\left(2 c L_f \sqrt{l_f^2+ l_g^2}\right)^{\frac{1}{ c l+1}} \cdot K^{\frac{ c l}{ c l+1}} \Big),$$
    \end{itemize}
    where $l_f$, $l_g$ and $L_f$, $L_g$ are smoothness constants and Lipschitz constants for $f$, $g$, respectively.
\end{theorem}

\begin{remark}
Note that the above stability bounds are independent of the specific form of the objective function $f(\cdot,\cdot)$ and the exact form of the sample distribution $\mathbb{D}_1$, which are more reliant on the properties of the loss functions and sample size $m_1$, and the stability bounds in the C-C and NC-NC cases
are related to the number of iterations additionally. Specifically, Part(a) establishes a stability bound of $\bigO(1/m_1)$ in the
SC-SC
setting and Part(b) considers a C-C
case with a stability bound  $\bigO( {C^{K/2}}/{m_1})$ related to the number of iterations and the data size, where $C$ is a constant.
The NC-NC
case is discussed in Part(c) which provides a stability bound of $\bigO(K^{\frac{cl}{cl+1}}/m_1)$, where $c$ is a constant to control the step size and $l$ is the larger smoothness number  of $l_f$ and $l_g$. The conclusions here match the existing results in minmax problems \citep{lei2021stability,farnia2021train}.
\end{remark}

\subsection{Stability Bounds for two-timescale SGD}
Compared with the above SSGD,
Two-timescale SGD (TSGD; Algorithm \ref{alg:TSGD}) always achieves more accurate approximate solutions by updating $\mathbf{y}$ multiple times before updating $\mathbf{x}$. In this section, we extend our analysis from SSGD to TSGD.
Particularly, compared with the results in \citet{bao2021stability}, we provide stability bounds in Theorem \ref{tsgd_stab} for the case where the inner level parameter ($\mathbf{y}$) is updated iteratively ({\em i.e.,} consistency).
We further explore in Theorem \ref{sc_nc} a particular NC-SC setting, which is commonly appeared in bilevel optimization applications such as meta learning and hyperparameter optimization.

\begin{theorem}
\label{tsgd_stab}
    Suppose that Assumptions \ref{asu1} and \ref{asu2} hold and $|g(\cdot,\cdot)| \leq 1$. Let $\mathbf{A}$ be the TSGD  algorithm with $K$ outer-iterations and $T$ inner-iterations. Then we have
    \begin{itemize}
         \item [(a)] Assume that Problem (\ref{pb1}) is SC-SC. Let  $\alpha_x = \alpha_y \leq \frac{1}{l_g}$ be the step sizes. Then,  $\mathbf{A}$ is $\beta$-argument-stable in expectation, where
         $$\beta \leq \bigO \left( {m_{1}}^{-1}\sqrt{L_f^2\alpha_x^2 + \Big(\frac{2T}{\alpha_y(2-\alpha_y l_g)}\Big)^2 }
        \left( 1- \alpha_x \mu_g + \alpha_x l T \right)^{K} \right).$$

        \item [(b)] Assume that Problem (\ref{pb1}) is C-C. Let $\alpha l = \max \{\alpha_x l_f, \alpha_y l_g\}$ and  $\alpha_x$, $\alpha_y \leq \frac{1}{l_g}$ be the step sizes. Then, $\mathbf{A}$ is $\beta$-argument-stable in expectation, where
        $$\beta \leq \bigO \left( m_{1}^{-1}\sqrt{L_f^2\alpha_x^2 + \Big(\frac{2T}{\alpha_y(2-\alpha_y l_g)}\Big)^2 }
        \left(2+2( \alpha l T)^2\right)^{K/ 2}\right). $$

        \item [(c)] Assume that Problem (\ref{pb1}) is NC-NC. Let the step sizes satisfy $\max  \left\{ \alpha_x,\alpha_y  \right\} \leq c/k$ for some constant $c \geq 0$ and $l = \max \left\{ l_f, l_g\right\}$. Then, $\mathbf{A}$ is $\beta$-argument-stable in expectation, where
        $$\beta \leq \bigO \left( ( m_1 T c l )^{-1}\left(2 c L_f \sqrt{l_f^2+T^2 l_g^2}\right)^{\frac{1}{T c l+1}} \cdot K^{\frac{T c l}{T c l+1}} \right).$$
    \end{itemize}

\end{theorem}

\begin{remark}
Compared with the previous results for SSGD, the stability bounds of TSGD depend on the number of iterations in the outer level loop, the number of iterations in the inner level loop, and the data size in the outer level loop. If the step sizes are sufficiently small,  we can see that the bounds in Theorem \ref{tsgd_stab} are asymptotically the same as the bounds of SSGD in Theorem \ref{ssgd_stab}. Thus, Theorem \ref{tsgd_stab} can be considered as a generalization of the previous one.  However, in comparison to Theorem \ref{alg:SSGD}, we can also see that our bounds are exponential in the inner loop iteration $T$. This is caused by the accumulation of the sub-sampled gradients in the inner level, which makes the   outer loop update  overly expansive. The exponential dependence on $T$ also reveals the difference
with the existing work of stability analysis for other problems, such as simple SGD and minmax problems.
\end{remark}

\begin{remark}
\label{remark_bao}
Comparing our results with the ones in
\citep{bao2021stability}, we have the following observations. 1) They only established the uniform stability bound for the Unrolled Differentiation algorithm \ref{alg:UD}, where the algorithm is reinitialized at each time entering the inner level loop, indicating that it takes into account the changes to only one parameter in the outer level loop, while our algorithm considers the update for both parameters.
2) Its proof needs to assume that the update of $y$ in the inner level after the reinitialization will not be affected by  the value specified for $\mathbf{x}$. However, this assumption is quite uncommon and is probably the reason that they do not need to make any assumption on the inner level objective function. 3) The term $\ell(\lambda_t^{\prime}, \hat{\theta}\left(\lambda_t^{\prime}, S^{t r}\right), z)$ in their paper, which corresponds to the inner loss function $g\left( \mathbf{x}_k^{\prime}, \mathbf{y}_k^{\prime}\left( \mathbf{x}_k^{\prime}, D_{m_2}\right), \eta  \right)$ in our paper, includes three arguments ({\em i.e.,} $\lambda_t^{\prime}$, $\hat{\theta}\left(\lambda_t^{\prime}, S^{t r}\right)$ and $z$), where $\hat{\theta}\left(\lambda_t^{\prime}, S^{t r}\right)$ is a function of
$\lambda_t^{\prime}$. This implies that  $\hat{\theta}\left(\lambda_t^{\prime}, S^{t r}\right)$ will change in the inner level updates as a result of the change of $\lambda_t^{\prime}$, and thus cannot be treated as an argument independent of $\lambda_t^{\prime}$, which is misused and causes technical flaws in their proof subsequently.
In contrast, our work allows the inner level parameters to be updated consistently ({\em i.e.,} carrying over the value in the last update), \textit{instead of being reinitialized at each time entering the inner level loop}. Specifically, we allow $y_k^{T}$ to be employed at the beginning of the $(k+1)$-th outer level iteration, rather than $y_0$. This enables us to obtain different stability bounds for different inner level objective functions.
\end{remark}

In the following, we extend our analysis to a particular NC-SC
setting that is frequently encountered in real-world applications and optimization analysis.

\begin{theorem}
\label{sc_nc}
    Suppose that Assumptions \ref{asu1} and \ref{asu2} hold, $0 \leq f(\cdot,\cdot) \leq 1$ and Problem (\ref{pb1}) is NC-SC.  Let $\mathbf{A}$ be the TSGD Algorithm with $K$ outer-iterations and $T$ inner-iterations with $\max  \left\{ \alpha_x,\alpha_y  \right\} \leq c/k$ for constant $c \geq 0$. Denote $l = \max \left\{ l_f, l_g\right\}$. Then, $\mathbf{A}$ is $\beta$-uniform-stable in expectation, where
    $$\beta \leq \bigO \left( \frac{\left(2 c L_f \sqrt{l_f^2+l_g^2 T^2}\right)^{\frac{1}{c (Tl+l-\mu_g)+1}} \cdot K^{\frac{c (Tl+l-\mu_g)}{c (Tl+l-\mu_g)+1}}\left(Tl+l-\mu_g +2/c \right)}{m_1 (Tl+l-\mu_g)} \right).$$

\end{theorem}

\begin{remark}
Compared with our previous analysis, we now sketch the technique differences in our analysis. We consider the bound of the term $(\delta_{\mathbf{x},k},\delta_{\mathbf{y},k}) ^T=(\|\mathbf{x}_k - \mathbf{x}_k^{\prime} \|, \| \mathbf{y}_k - \mathbf{y}_k^{\prime}\|)^T$, while we employ $\delta_k = \sqrt{\left\|\mathbf{x}_k-\mathbf{x}_k^{\prime}\right\|_{2}^{2}+\left\|\mathbf{y}_k-\mathbf{y_k}^{\prime}\right\|_{2}^{2}}$ in the previous analysis, where $(\mathbf{x}_k,\mathbf{y}_k)$, $(\mathbf{x}_k^{\prime},\mathbf{y}_k^{\prime})$ are the outputs of TSGD after $k$ iterations for $D_{m_1}$ and $D_{m_1}^{\prime}$ respectively with $D_{m_1}$ and $D_{m_1}^{\prime}$ differing  in one sample. In the NC-SC
setting, we show that $(\delta_{\mathbf{x},k+1},\delta_{\mathbf{y},k+1}) ^T \leq ((1+\alpha_x l)\delta_{\mathbf{x},k},(1+\alpha_x Tl)\delta_{\mathbf{y},k})^T$ ($\leq$ means the entry-wise inequality), which means our term can be controlled. Then,  we take the expectation of it to derive our uniform stability bound. To achieve the generalization error over continuously changing parameters,
it is imperative to take into account the growth of $(\delta_{\mathbf{x},k},\delta_{\mathbf{y},k})$ instead of $\delta_{\mathbf{x},k}$ in \citep{bao2021stability}. Appendix \ref{thm_tsgd} provides more details.
\end{remark}

Thus, based on our previous results, we now provide the first generalization bounds in the NC-NC setting for both SSGD and TSGD.

\begin{corollary}
\label{ssgd_gen}
    Assume that the problem is NC-NC, $|f(\cdot, \cdot ; \xi)| \leq 1$ for all $\xi$, and  Assumptions \ref{asu1} and \ref{asu2} hold. Denote $l= \max \{l_f,l_g\}$ with $\max \{\alpha_x, \alpha_y\} \leq c/k$ for constant $c \geq 0$. Then, the generalization error of SSGD \ref{alg:SSGD} with $K$ iterations is bounded by
    $\bigO( K^{\frac{cl}{cl+1}}/{m_1} )$.
\end{corollary}

\begin{corollary}
\label{tsgd_gen}
    Assume that the problem is NC-NC, $|f(\cdot, \cdot ; \xi)| \leq 1$ for all $\xi$, and Assumptions \ref{asu1} and \ref{asu2} hold. Let $l= \max \{l_f,l_g\}$ with $\max \{\alpha_x, \alpha_y\} \leq c/k$. Then, the generalization error of TSGD \ref{alg:TSGD} with $K$ outer iterations and $T$ inner iterations is bounded by
    $\bigO( T^{\frac{1}{T c l+1}} K^{1-\frac{1}{T c l+1}}/({m_1}))$.
\end{corollary}

\begin{remark}
By Theorem \ref{thm1_connection}, \ref{ssgd_stab}, \ref{tsgd_stab}], we can derive the above corollaries on generalization error from  stability bounds. Corollary \ref{ssgd_gen} and Corollary \ref{tsgd_gen} show that extremely high number of iterations ($K$ for SSGD and $K$,$T$ for TSGD) will drastically reduce the stability of these algorithms and increase the generalization error, which will make these algorithms increase the risk of overfitting. We will also verify it in the following experiments.
\end{remark}
\section{Experiments}
\label{experiment}
In this section, we empirically validate our previous theoretical results on real world datasets. Two experiments, including meta-learning and hyperparameter optimization, are conducted via Algorithm \ref{alg:TSGD} TSGD (note that when $T=1$, TSGD is just SSGD). Due to the space limitation, we just present the meta learning experiment here, leaving the hyperparameter optimization experiment and other details in the Appendix \ref{additional experiment}.
\subsection{Meta Learning}
\label{exp:meta learning}
Consider the few-shot meta-learning problem with $M$ tasks $\{ \mathcal{T}_i, i=1, ..., M \}$ sampled from distribution $P_{\mathcal{T}}$.
We aim to learn a model that can rapidly adapt to different tasks.
Firstly, the embedding model $\phi$ is shared by all tasks to learn embedded features. Secondly, the task-specific parameter $w_i$ is to adapt the shared embedding to its own sub-problem. Thus, the overall problem of meta-learning can be formulated as follow:
\begin{subequations}
\begin{align}
  &\min _\phi \mathcal{L}_{\mathcal{D}}\left(\phi, \bar{w}^*\right)=\mathbb{E}_{\xi \in \mathcal{D}_i^{\mathrm{te}},\mathcal{T}_i }\left[\mathcal{L}\left(\phi, w_i^* ; \xi\right) \right], \label{eq:meta_learner}\\
   &\text { s.t. } \bar{w}^*=\underset{\bar{w}}{\arg \min } \Big[\mathcal{L}_{\mathcal{D}^{\mathrm{tr}}}(\phi, \bar{w})=\mathbb{E}_{\mathcal{T}_i }\left[ \mathcal{L}_{\mathcal{D}_i^{\mathrm{tr}}}\left(\phi, w_i\right)\Big]\right]. \label{eq:base_learner}
\end{align}
\end{subequations}
\vspace{-0.15in}

where $\mathcal{D}^{\mathrm{tr}}_i$ and $\mathcal{D}^{\mathrm{te}}_i$ are the training and testing datasets for task $\mathcal{T}_i$. Each $w_i$ is computed from one or more gradient descent updates from $\bar{w}$ on the corresponding task (rapid adaptation), {\em i.e.,} $w_i=\bar{w}-\alpha \nabla_{\bar{w}} \mathcal{L}_{\mathcal{D}_{tr}}(\phi, w_i)$.
In the inner level, the base learner optimizes the series of ${w_i}$ for each tasks (Equation \ref{eq:base_learner}). In the outer level, the meta-learner optimizes the embedding model $\phi$ using the minimizers $w^*_i$ learned from the inner level and computes the loss from the testing dataset (Equation \ref{eq:meta_learner}).

\noindent {\bf Settings and Implementation}
We evaluate the behavior of the 5-way-1-shot task on the Omnilot dataset \citep{lake2015human}, {\em i.e.,} it aims to classify 5 unseen classes from only 1 labeled sample. It contains 1623 different handwritten characters from 50 different alphabets. The image is in greyscale with a size $28\times28$. We follow similar settings in \cite{ji2021bilevel}. A five-layer fully-connected network is constructed, where the task-specific parameter $w_i$ corresponds to the last layer of the network and the shared embedding model $\phi$ corresponds to all preceding layers. Thus, we train the two sets of layers separately in the outer and inner level of optimization. We build our model and establish our training using the software library {\fontfamily{cmss}\selectfont learn2learn} \citep{Arnold2020}. We follow the official train-validation-test partition and train $\phi$, $w_i$ using the training set. The size of each layer in the network is $784 \rightarrow 256 \rightarrow 128 \rightarrow 64 \rightarrow 64 \rightarrow 5$. We set the number of tasks for training and testing set to 2000 and the batch size of tasks to 32. The learning rate of $\phi$ and $w_i$ are 0.002 and 0.01, respectively. Results are evaluated based on the average of 5 trial runs with different random seeds.

\noindent {\bf Results Evaluation}
Figure \ref{fig:metalearning} presents the learning curves on training set, testing set and the generalization error with different values of inner iterations $T$ and outer iterations $K$. Generalization error is estimated by the difference between training and testing loss. On one hand, it can be seen that the model easily overfits on the testing set as $K$ increases drastically (Figure \ref{fig:metatestloss}) and the effect of $T$ is very limited. On the other hand, with an appropriate value of $K$, smaller $T$ (\textit{i}.\textit{e}) will result in underfitting on the testing loss ($T=1$ in the Figure \ref{fig:metagenloss} causes highest generalization error due to the underfitting training process). The trend of generalization error in terms of $K$ and $T$ matches with our analysis in Theorem \ref{sc_nc}.

\begin{figure}[]
     \centering
     \begin{subfigure}[1]{0.32\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/meta_train_loss_2500.pdf}
         \caption{Meta-training Loss}
         \label{fig:metatrainloss}
     \end{subfigure}
     \hfill
     \begin{subfigure}[2]{0.32\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/meta_test_loss_2500.pdf}
         \caption{Meta-testing Loss}
         \label{fig:metatestloss}
     \end{subfigure}
     \hfill
     \begin{subfigure}[3]{0.32\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/meta_gap_loss_2500.pdf}
         \caption{Generalization Error}
         \label{fig:metagenloss}
     \end{subfigure}
        \caption{Results of meta learning with various values of $T$ and $K$}
        \label{fig:metalearning}
\end{figure}

\section{Conclusion}
\label{conclusion}
We give a thorough analysis on the generalization of first-order (gradient-based) methods for the bilevel optimization framework. In particular, we establish a quantitative connection between generalization and algorithmic stability and provide the first generalization bounds of the continuous updates for inner parameters and outer parameters in multiple settings. Our experiments suggest that inappropriate iterations will cause underfitting and overfitting easily. The tendency of generalization error also validates our theoretical results.

From the discussion in previous sections, we only discussed the first-order method, while there exists a number of estimating second-order and momentum-based approaches to solve the bilevel optimization problem. Dealing with the approximation of hypergradient in generalization analysis is another direction for future work.

redpajama_set_name: RedPajamaArXiv


Why do we find it so easy to retaliate? What is it about human nature that makes us want to make sure the person who hurts us is hurt even more? Even if we don't act on it, we still wish it and that's every bit as bad in God's eyes. Nobody denies the difficulty in turning the other cheek, but that is exactly what Jesus taught. You may be thinking, I know that somewhere in the Bible it says an eye for an eye. And you would be right, but the Old Testament intent for "an eye for an eye" was that of proportional justice, the idea that no matter the social class, the punishment would fit the crime. But sinful people, as they often do, take the intent of a passage, use it out of context, and misrepresent it. Such was the case in the circumstance of the Old Testament references to "an eye for an eye" approach. It is and never was meant to be a license for payback.
Jesus said, "To one who strikes you on the cheek, offer the other also" (Luke 6:29). No, Jesus wasn't kidding! In fact, in the verses just prior, He said, "Love your enemies, do good to those who hate you, bless those who curse you, pray for those who abuse you" (Luke 6:27–28).
Today, as you go through your day and people strike you on the cheek, some lightly and some not so lightly, resist the temptation to retaliate. To love your enemies is a command. To turn the other cheek is a command. But most importantly is who it is that gave the command. He is the one who expressed love for His enemies and a willingness to be humiliated by turning the other cheek like no other. His name is Jesus.
Thank You, heavenly Father for Your Word. I confess my failure to take Your Word seriously on this issue. Forgive me and help me do what is unnatural in my own strength to do. Help me to please You and love those who strike out at me. After all, You loved me first. Amen.
From Whom Are We Now Seeking Our Independence?
On Thursday, we celebrate the independence of our nation. As the early settlers charted the course of our nation, they had no shame concerning the things of God. The first line in the Mayflower Compact was, "In the name of God, Amen" and when the Pilgrims' signed it in 1620, they acknowledged God's sovereignty and sought to glorify Him. Our founding fathers believed the truths of the Bible and acknowledged their faith in God as they wrote the Constitution. Many prominent universities were founded on the Word of God. God is in our pledge and it is "in Him" we claim to trust on our currency. But given the current state of affairs in our country, can we say that He is in our hearts? As we reflect this July 4th, we might ask, from whom are we now seeking our independence?
"Righteousness exalts a nation, but sin is a reproach to any people" Proverbs is a collection of wise sayings. Proverbs 14:34 teaches that actions aligned with God's standards exalt or "lift up" a nation while sin is a disgrace. The United States has been richly blessed. Throughout our history, we have taken up the cause for much that is right and good. But have we, in our pride, forgotten what David wrote in Psalm 8:4, "What is man that you are mindful of him", that as great as God is, it's amazing He would take notice of and care for us. Have we really come to believe that "man is the measure of all things", that we determine what is right and true? Perhaps we have. Today, many deny His sovereignty. Many of those educational institutions that were founded on the Word of God are now major influencers in denying Him. Our Creator, whom we used to hold close, we now appear to be, in increasing measure, pushing away. It will make you wonder if there is any reason for hope for our nation.
But there is hope! There is always hope because our hope lies, not in our government or our courts, but in a sovereign God. Nothing that happens ever happens outside of His control. We hope, knowing God will preserve His people through times such as these. We hope, knowing that God is building His church and nothing can stop Him.
It's easy to get discouraged, confused and quite frankly, angry by recent events and our nations progressive departure from God's will. But it's also important to remember the greater the darkness the greater the opportunity for the light of the gospel of Christ to shine. The opportunity and privilege we have as Christians is to display what the gospel has done in our own lives, to remember who we were before Christ saved us (an enemy of God), and to let what He has graciously done in our hearts motivate us in our witness. The gospel is not only for those whose sins are different from ours, but for all of our sins. Pray that in a time of increasing darkness the gospel would go forth and trust that it would bring a harvest. Pray that our nation would realign with the will of God and instead of seeking our independence from Him, we would reaffirm our dependence on Him. "If my people who are called by name humble themselves, and pray and seek my face and turn from their wicked ways, then I will hear from heaven and will forgive their sin and heal their land" (2 Chronicles 7:14).
Father, as a nation we have strayed. We have strayed in so many ways on so many issues. Forgive us. Help us to pursue righteousness and let Your glory be our desire. Thank You God for forgiving me for my sin. Thank You for Jesus Christ whose perfect atonement cleanses all sin. Let our nation exalt Him above all. Thank You for Your grace and mercy. Let us be a people who seek to please You alone. Amen!


Centrelink, Australia's social security office will start sending SMS messages to clients, reminding them to attend job interviews. Of course, the next logical progression from this (which the Liberal Party should love) would be to integrate this with mobile phone tracking/triangulation, to keep track of dole recipients; much like those bracelets they put on home-detention inmates, only the client pays for the technology. If they leave the 5km radius of their home or registered appointments, or enter an off-limits area during business hours, their benefits get cut off until they wait in queue, grovel and fill in forms. After all, dole recipients are parasites, the philosophy goes, and the experience of receiving welfare should be made as humiliating as possible. And with the government's new telecommunications intereception bill, giving a wide range of government agencies the power to intercept people's communications without a warrant in the name of Fighting Terrorism, this could be possible.
So what did I do this weekend, you ask? Well, today I spent much of the day installing Debian 3.0 on my Linux box. It looks quite doovy (it has a lot of packages), and it's about time I replaced the patched-over RedHat 6.x setup it was.


text:
@interface GHPageModel : NSObject

@property (nonatomic, strong) NSString * title;
@property (nonatomic, strong) NSString * url;
@property (nonatomic, strong) NSString * category;
@property (nonatomic, strong) NSString * type;
@property (nonatomic, strong) NSNumber * priority;
@property (nonatomic, strong) NSString * identifier;
@property (nonatomic, strong) NSArray  * children;

@end

redpajama_set_name: RedPajamaGithub


Next HIGH TIDE in Pagadian is at , which is in from now.
Next LOW TIDE in Pagadian is at , which is in from now.
The local time in Pagadian is . See the detailed Pagadian tide chart and tide table below.
Tide Times are +08 (UTC +8.0hrs). Last Spring High Tide at Pagadian was on Sun 21 Apr (height:2.42m 7.9ft). Next high Spring Tide at Pagadian will be on Tue 07 May (height:2.34m 7.7ft).
Pagadian Tide Chart. The largest known tidal range at Pagadian is 2.68m 8.8 feet.
The tide chart above shows the height and times of high tide and low tide for Pagadian. The red flashing dot shows the tide time right now. The yellow shading corresponds to daylight hours between sunrise and sunset at Pagadian.
NOTE: Use of this site for Pagadian tide times is subject to our terms.