A New Simple Stochastic Gradient Descent Type Algorithm With Lower Computational Complexity for Bilevel Optimization

Haimei Huo, Risheng Liu, , Zhixun Su Haimei Huo is with the School of Mathematical Sciences, Dalian University of Technology, Dalian 116024, China (e-mail: ab1234@mail.dlut.edu.cn).Risheng Liu(Corresponding author) is with the DUT-RU International School of Information Science and Engineering, Dalian University of Technology, and also with the Key Laboratory for Ubiquitous Network and Service Software of Liaoning Province, Dalian 116024, China (e-mail: rsliu@dlut.edu.cn).Zhixun Su(Corresponding author) is with the School of Mathematical Sciences, Dalian University of Technology, with the State Key Laboratory of Structural Analysis for Industrial Equipment, and also with the Key Laboratory for Computational Mathematics and Data Intelligence of Liaoning Province, Dalian 116024, China (e-mail: zxsu@dlut.edu.cn).

Abstract

Bilevel optimization has been widely used in many machine learning applications such as hyperparameter optimization and meta learning. Recently, many simple stochastic gradient descent(SGD) type algorithms(without using momentum and variance techniques) have been proposed to solve the bilevel optimization problems. However, all the existing simple SGD type algorithms estimate the hypergradient via stochastic estimation of Neumann series. In the paper, we propose to estimate the hypergradient via SGD-based Estimation(i.e., solving the linear system with SGD). By using warm start initialization strategy, a new simple SGD type algorithm SSGD based on SGD-based Estimation is proposed. We provide the convergence rate guarantee for SSGD and show that SSGD outperforms the best known computational complexity achieved by the existing simple SGD type algorithms. Our experiments validate our theoretical results and demonstrate the efficiency of our proposed algorithm SSGD in hyperparameter optimization applications.

Index Terms:

Bilevel optimization, stochastic gradient descent(SGD)-based Estimation, simple SGD type algorithms, warm start.

I Introduction

Bilevel optimization has become a powerful tool for many machine learning applications such as meta learning[1, 2, 3, 4], hyperparameter optimization[5], and reinforcement learning[6, 7]. It has a nested structure which contains two levels of optimization tasks, and the feasible region of the upper level(UL) function is restricted by the optimal solutions of the lower level(LL) problem. Generally, bilevel optimization can be expressed as follows:

		$\displaystyle\min\limits_{\boldsymbol{x}\in\mathbb{R}^{p}}\Phi(\boldsymbol{x}):=f(\boldsymbol{x},\boldsymbol{y}^{*}(\boldsymbol{x}))$		(1)
		$\displaystyle s.t.~{}\boldsymbol{y}^{*}(\boldsymbol{x})\in\mathcal{S}(\boldsymbol{x}):={\underset{\boldsymbol{y}\in\mathbb{R}^{q}}{\operatorname{arg\,min}}\,g(\boldsymbol{x},\boldsymbol{y})}$

where $f(\boldsymbol{x},\boldsymbol{y})$ is the UL function, $g(\boldsymbol{x},\boldsymbol{y})$ is the LL function, and $\mathcal{S}(\boldsymbol{x})$ is the solution set of the LL problem.

Many deterministic algorithms have been proposed to solve the bilevel problem (1). In early works, the bilevel problem is reformulated into single level problem using the first order optimality conditions of the LL problem[8, 9, 10, 11]. However, such type of methods often involve a large number of constraints, and are not applicable in complex machine learning problems. To address the problem, gradient based algorithms have been revisited to solve problem (1). Depending on the way to approximate the hypergradient $\nabla\Phi(\boldsymbol{x})$ , these algorithms can be generally categorized into the iterative differentiation(ITD) based approach[1, 5, 12, 13, 14] and the approximate implicit differentiation(AID) based approach[12, 15, 16, 17].

ITD based approach treats the iterative optimization algorithm for the LL problem as a dynamical system, and approximates the hypergradient by back-propagation(BP) along the dynamical system [12, 18], this kind of hypergradient approximation method termed as BP in the paper. Since computing the gradient along the entire dynamical system is computationally challenging for high-dimensional problems, [5] proposed to truncate the path to compute the gradient. On the other hand, the convergence of ITD based approach for solving problem $(\ref{eq1})$ has also been studied. [1] assumed that the LL function $g(\boldsymbol{x},\cdot)$ is strongly convex, and provided the asymptotic convergence analysis. To address the case that the LL function $g(\boldsymbol{x},\cdot)$ is convex, [19] proposed BDA and the asymptotic convergence analysis is provided. For the nonasymptotic convergence analysis of ITD based approach, see [14].

AID based approach applies implicit function theorem to the first order optimality conditions of the LL problem, and approximates the hypergradient by solving a linear system[15, 16]. There are mainly two methods to solve the linear system: conjugate gradient descent(CG)[12] and approximating the Hessian inverse with Neumann series[17], these two hypergradient approximation methods termed as CG and Neumann series(NS) in the paper. For the nonasymptotic convergence analysis of the AID based approach, see [14, 17].

In fact, in many machine learning applications, the bilevel optimization takes the following form:

		$\displaystyle\min\limits_{\boldsymbol{x}\in\mathbb{R}^{p}}\Phi(\boldsymbol{x}):=f\left(\boldsymbol{x},\boldsymbol{y}^{}(\boldsymbol{x})\right):=\mathbb{E}_{\boldsymbol{\xi}}\left[F(\boldsymbol{x},\boldsymbol{y}^{}(\boldsymbol{x});\boldsymbol{\xi})\right]$		(2)
		$\displaystyle s.t.~{}\boldsymbol{y}^{*}(\boldsymbol{x}):={\underset{\boldsymbol{y}\in\mathbb{R}^{q}}{\operatorname{arg\,min}}\,g(\boldsymbol{x},\boldsymbol{y})}:=\mathbb{E}_{\boldsymbol{\zeta}}[G(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\zeta)}]$

where UL function $f(\boldsymbol{x},\boldsymbol{y})$ and LL function $g(\boldsymbol{x},\boldsymbol{y})$ take the expected values with respect to(w.r.t.) random variables $\boldsymbol{\xi}$ and $\boldsymbol{\zeta}$ , respectively; $g(\boldsymbol{x},\cdot)$ is strongly convex w.r.t. $\boldsymbol{y}$ and $\Phi(\boldsymbol{x})$ is nonconvex w.r.t. $\boldsymbol{x}$ . Notice that if $\boldsymbol{\xi}$ and $\boldsymbol{\zeta}$ are discrete random variables, the values of $\boldsymbol{\xi}$ and $\boldsymbol{\zeta}$ are chosen from the given data sets $\mathcal{D}_{n}:=\{\boldsymbol{\xi}_{1},\ldots,\boldsymbol{\xi}_{n}\}$ and $\mathcal{D}_{m}:=\{\boldsymbol{\zeta}_{1},\ldots,\boldsymbol{\zeta}_{m}\}$ , respectively, and the probability to sample data from $\mathcal{D}_{n}$ and $\mathcal{D}_{m}$ is equal, respectively, then the objective functions in problem (2) can be written as $f(\boldsymbol{x},\boldsymbol{y}):=\frac{1}{n}\sum_{i=1}^{n}F(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\xi}_{i})$ and $g(\boldsymbol{x},\boldsymbol{y}):=\frac{1}{m}\sum_{i=1}^{m}G(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\zeta}_{i})$ .

In order to achieve better efficiency than deterministic algorithms for solving problem (2), many stochastic algorithms have been proposed. Without considering the loop to estimate the Hessian inverse, these existing stochastic algorithms for problem (2) can be categorized into single loop algorithms and double loop algorithms according to the number of loops.

For single loop algorithms, the upper variable $\boldsymbol{x}$ and the LL variable $\boldsymbol{y}$ are updated simultaneously. As far as we know, the first single loop algorithm proposed for solving problem (2) is the two-timescale stochastic approximation(TTSA) algorithm[7], which is a simple SGD type algorithm. To ensure convergence, a larger step size for the LL problem and a smaller step size for the UL problem are required. Moreover, to reach an $\epsilon$ -stationary point, the computational complexity is $\widetilde{\mathcal{O}}(\epsilon^{-5/2})$ in terms of the target accuracy $\epsilon$ . Recently, based on momentum and variance techniques, some single loop algorithms with improved computational complexity have been proposed, see e.g., [20, 21].

For a double loop algorithm, it contains two loops in which the inner loop is used to update the LL variable $\boldsymbol{y}$ to obtain an approximation to the optimal solution $\boldsymbol{y}^{*}(\boldsymbol{x})$ of the LL problem, and the outer loop is used to update the UL variable $\boldsymbol{x}$ . Following the direction, [17] proposed a bilevel stochastic approximation(BSA) algorithm and the computational complexity is $\widetilde{\mathcal{O}}(\epsilon^{-3})$ in terms of the target accuracy $\epsilon$ . By introducing warm start initializaiton strategy, stocBiO[14] achieves the computational complexity of $\widetilde{\mathcal{O}}(\epsilon^{-2})$ in terms of the target accuracy $\epsilon$ . In addition, in order to improve the computational complexity of the simple SGD type double loop algorithms, variance reduction is used for both UL variable $\boldsymbol{x}$ and LL variable $\boldsymbol{y}$ in [22], and the computational complexity is $\widetilde{\mathcal{O}}(\epsilon^{-1.5})$ in terms of the target accuracy $\epsilon$ .

Although many stochastic algorithms have been proposed to solve problem (2), little attention is paid to the the hypergradient estimation process except a few attempts recently. [23] noticed that to estimate the hypergradient, the existing single loop algorithms and double loop algorithms require an additional loop to estimate Hessian inverse. By studying the well known hypergradient approximation methods such as BP[12], CG[12], and NS[17], they identified a general approximation formulation of hypergradient computation. Based on the formulation and the momentum based variance reduction technique used in [24], they proposed a fully single loop algorithm(FSLA) in which an additional loop to estimate the Hessian inverse is not required.

In the paper, our goal is to design a new simple stochastic gradient descent(SGD) type algorithm for problem (2), which can achieve lower computational complexity compared with the existing simple SGD type algorithms(simple SGD type algorithms refer to stochastic algorithms that do not use momentum and variance techniques). We notice that the existing simple SGD type algorithms(i.e., BSA, TTSA, and stocBiO) for solving problem (2) estimate the hypergradient via stochastic estimation of Neumann series, this kind of hypergradient estimation method termed as Stochastic NS in the paper. Inspired by the work in [23], we study some hypergradient estimation methods: Stochastic BP, Stochastic NS, and SGD-based Estimation, and find that the algorithm has the potential to converge faster by using SGD-based Estimation to estimate the hypergradient. Based on the finding, by using warm start initialization strategy, a new simple SGD type algorithm SSGD based on SGD-based Estimation is proposed. Through our theoretical analysis, SSGD can achieve lower complexity compared with the existing simple SGD type algorithms. The main contributions of this paper are:

1.

We propose to use SGD-based Estimation to estimate the hypergradient.
2.

We propose a new simple SGD type algorithm SSGD.
3.

We prove that SSGD can achieve lower computational complexity compared with the existing simple SGD type algorithms. As shown in Table I, the computational complexity of SSGD to reach an $\epsilon$ -stationary point outperforms that of BSA and stocBiO by an order of $\mathcal{O}(\epsilon^{-1}\log({1}/{\epsilon}))$ and $\mathcal{O}(\log({1}/{\epsilon}))$ , respectively. In addition, in terms of the target accuracy $\epsilon$ , the computational complexity of SSGD outperforms that of TTSA by an order of $\mathcal{O}(\epsilon^{-0.5}\log({1}/{\epsilon}))$ .
4.

We perform experiments to validate our theoretical results and demonstrate the application of our algorithm in hyperparameter optimization.

The rest of this paper is organized as follows. In Section II, we first make some settings for problem (2), and specify the problem to study in this paper. Then, for problem (2), we introduce some hypergradient estimation methods and analyze them. In Section III, we propose a new simple SGD type algorithm for problem (2). The convergence and complexity results of our proposed algorithm are provided in Section IV. Experimental results are provided in Section V. In Section VI, we conclude and summarize the paper.

Notations. We use $\|\cdot\|$ to denote the $l_{2}$ norm for vectors and spectral norm for matrices. For a twice differentiable function $f(\boldsymbol{x},\boldsymbol{y}):\mathbb{R}^{p}\times\mathbb{R}^{q}\rightarrow\mathbb{R}$ , $\nabla f(\boldsymbol{x},\boldsymbol{y})$ denotes the gradient of $f(\boldsymbol{x},\boldsymbol{y})$ taken w.r.t. all the variables $(\boldsymbol{x},\boldsymbol{y})$ , $\nabla_{\boldsymbol{x}}f(\boldsymbol{x},\boldsymbol{y})$ (resp. $\nabla_{\boldsymbol{y}}f(\boldsymbol{x},\boldsymbol{y})$ ) denotes its partial derivate taken w.r.t. $\boldsymbol{x}$ (resp. $\boldsymbol{y}$ ), $\nabla_{\boldsymbol{x}}\nabla_{\boldsymbol{y}}f(\boldsymbol{x},\boldsymbol{y})$ denotes the Jacobian of $\nabla_{\boldsymbol{x}}f(\boldsymbol{x},\boldsymbol{y})$ w.r.t. $\boldsymbol{y}$ , $\nabla_{\boldsymbol{y}}^{2}f(\boldsymbol{x},\boldsymbol{y})$ denotes the Hessian matrix of $f(\boldsymbol{x},\boldsymbol{y})$ w.r.t. $\boldsymbol{y}$ . Furthermore, for a twice differentiable function $f(\boldsymbol{x},\boldsymbol{y}):=\mathbb{E}_{\boldsymbol{\xi}}[F(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\xi})]:\mathbb{R}^{p}\times\mathbb{R}^{q}\rightarrow\mathbb{R}$ , which takes expected value w.r.t. random variable $\boldsymbol{\xi}$ , we define

\nabla_{\boldsymbol{x}}F(\boldsymbol{x},\boldsymbol{y};\mathcal{D})=\frac{1}{|\mathcal{D}|}\sum\limits_{i=1}^{|\mathcal{D}|}\nabla_{\boldsymbol{x}}F(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\xi}_{i})

where $\mathcal{D}=\{\boldsymbol{\xi}_{1},\ldots,\boldsymbol{\xi}_{|\mathcal{D}|}\}$ is a sample set with samples $\boldsymbol{\xi}_{i}$ , $i=1,\ldots,|\mathcal{D}|$ sampled over the distribution of $\boldsymbol{\xi}$ , and $|\mathcal{D}|$ is the sample size of $\mathcal{D}$ . Similarly, the definitions of $\nabla_{\boldsymbol{y}}F(\boldsymbol{x},\boldsymbol{y};\mathcal{D})$ , $\nabla_{\boldsymbol{x}}\nabla_{\boldsymbol{y}}F(\boldsymbol{x},\boldsymbol{y};\mathcal{D})$ , and $\nabla_{\boldsymbol{y}}^{2}F(\boldsymbol{x},\boldsymbol{y};\mathcal{D})$ can be obtained.

TABLE I: Computational complexity of the simple SGD type algorithms to reach an

\epsilon

-stationary point of

\Phi(\boldsymbol{x})

for problem (2)(the computational complexity of BSA, TTSA, and stocBiO can be obtained from Table 1 in [25]). Note that

\kappa

denotes the condition number; the notation

\kappa^{p}

denotes a polynomial function of

\kappa

since [7] does not provide the dependence on the condition number

\kappa

Algorithm	loop	batch size	computational complexity
BSA[17]	double	$\mathcal{O}(1)$	$\mathcal{O}(\kappa^{9}\epsilon^{-3}\log(1/\epsilon))$
TTSA [7]	single	$\mathcal{O}(1)$	$\mathcal{O}(\kappa^{p}\epsilon^{-{5}/{2}}\log(1/\epsilon))$
stocBiO [14]	double	$\mathcal{O}(\epsilon^{-1})$	$\mathcal{O}(\kappa^{9}\epsilon^{-2}\log(1/\epsilon))$
SSGD(Ours)	double	$\mathcal{O}(\epsilon^{-1})$	$\mathcal{O}(\kappa^{9}\epsilon^{-2})$

II Analysis

In this section, we first make some settings for problem (2), and then introduce some hypergradient estimation methods for problem (2). Finally, we analyze these hypergradient estimation methods.

II-A Preliminaries

In the following, we make some assumptions on the objective functions in problem (2).

Assumption 1.

The functions $g(\boldsymbol{x},\boldsymbol{y})$ and $G(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\zeta})$ are $\mu$ -strongly convex w.r.t. $\boldsymbol{y}$ for any $\boldsymbol{x}$ and $\boldsymbol{\zeta}$ , and $\Phi(\boldsymbol{x})$ is nonconvex w.r.t. $\boldsymbol{x}$ .

Assumption 2.

The functions $f(\boldsymbol{x},\boldsymbol{y})$ and $g(\boldsymbol{x},\boldsymbol{y})$ satisfy

•

$f(\boldsymbol{x},\boldsymbol{y})$ is $M$ -Lipschitz, i.e., for $\boldsymbol{z}_{1}:=[\boldsymbol{x}_{1};\boldsymbol{y}_{1}]$ , $\boldsymbol{z}_{2}:=[\boldsymbol{x}_{2};\boldsymbol{y}_{2}]$ , we have

\|f(\boldsymbol{x}_{1},\boldsymbol{y}_{1})-f(\boldsymbol{x}_{2},\boldsymbol{y}_{2})\|\leq M\|\boldsymbol{z}_{1}-\boldsymbol{z}_{2}\|.

•

$\nabla f(\boldsymbol{x},\boldsymbol{y})$ and $\nabla g(\boldsymbol{x},\boldsymbol{y})$ are $L$ -Lipschitz, i.e., for $\boldsymbol{z}_{1}:=[\boldsymbol{x}_{1};\boldsymbol{y}_{1}]$ , $\boldsymbol{z}_{2}:=[\boldsymbol{x}_{2};\boldsymbol{y}_{2}]$ , we have

	$\displaystyle\\|\nabla f(\boldsymbol{x}_{1},\boldsymbol{y}_{1})-\nabla f(\boldsymbol{x}_{2},\boldsymbol{y}_{2})\\|\leq L\\|\boldsymbol{z}_{1}-\boldsymbol{z}_{2}\\|,$
	$\displaystyle\\|\nabla g(\boldsymbol{x}_{1},\boldsymbol{y}_{1})-\nabla g(\boldsymbol{x}_{2},\boldsymbol{y}_{2})\\|\leq L\\|\boldsymbol{z}_{1}-\boldsymbol{z}_{2}\\|.$

•

$\nabla_{\boldsymbol{x}}\nabla_{\boldsymbol{y}}g(\boldsymbol{x},\boldsymbol{y})$ and $\nabla_{\boldsymbol{y}}^{2}g(\boldsymbol{x},\boldsymbol{y})$ are $\tau$ - and $\rho$ -Lipschitz, i.e., for $\boldsymbol{z}_{1}:=[\boldsymbol{x}_{1};\boldsymbol{y}_{1}]$ , $\boldsymbol{z}_{2}:=[\boldsymbol{x}_{2};\boldsymbol{y}_{2}]$ , we have

	$\displaystyle\\|\nabla_{\boldsymbol{x}}\nabla_{\boldsymbol{y}}g(\boldsymbol{x}_{1},\boldsymbol{y}_{1})-\nabla_{\boldsymbol{x}}\nabla_{\boldsymbol{y}}g(\boldsymbol{x}_{2},\boldsymbol{y}_{2})\\|\leq\tau\\|\boldsymbol{z}_{1}-\boldsymbol{z}_{2}\\|,$
	$\displaystyle\\|\nabla_{\boldsymbol{y}}^{2}g(\boldsymbol{x}_{1},\boldsymbol{y}_{1})-\nabla_{\boldsymbol{y}}^{2}g(\boldsymbol{x}_{2},\boldsymbol{y}_{2})\\|\leq\rho\\|\boldsymbol{z}_{1}-\boldsymbol{z}_{2}\\|.$

Furthermore, the same assumption holds for $F(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\xi})$ and $G(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\zeta})$ for any given $\boldsymbol{\xi}$ and $\boldsymbol{\zeta}$ .

In fact, for problem $(\ref{eq2})$ , Assumptions 1 and 2 are widely adopted to establish the convergence of the stochastic algorithms such as BSA, TTSA, and stocBiO. Based on Assumptions 1 and 2, the differentiability of $\Phi(\boldsymbol{x})$ , and the Lipschitz continuity of $\boldsymbol{y}^{*}(\boldsymbol{x})$ and $\nabla\Phi(\boldsymbol{x})$ can be ensured, as shown below.

Proposition 1.

For problem (2), suppose Assumptions 1 and 2 hold. Then $\Phi(\boldsymbol{x})$ is differentiable, and the hypergradient $\nabla\Phi(\boldsymbol{x})$ takes the form of

\nabla\Phi(\boldsymbol{x})=\nabla_{\boldsymbol{x}}f(\boldsymbol{x},\boldsymbol{y}^{*}(\boldsymbol{x}))-\nabla_{\boldsymbol{x}}\nabla_{\boldsymbol{y}}g(\boldsymbol{x},\boldsymbol{y}^{*}(\boldsymbol{x}))\boldsymbol{v}(\boldsymbol{x},\boldsymbol{y}^{*}(\boldsymbol{x}))

(3)

where $\boldsymbol{v}(\boldsymbol{x},\boldsymbol{y}^{*}(\boldsymbol{x}))=[\nabla_{\boldsymbol{y}}^{2}g(\boldsymbol{x},\boldsymbol{y}^{*}(\boldsymbol{x}))]^{-1}\nabla_{\boldsymbol{y}}f(\boldsymbol{x},\boldsymbol{y}^{*}(\boldsymbol{x}))$ is the solution of the following linear system:

\nabla_{\boldsymbol{y}}^{2}g\left(\boldsymbol{x},\boldsymbol{y}^{*}(\boldsymbol{x})\right)\boldsymbol{v}=\nabla_{\boldsymbol{y}}f\left(\boldsymbol{x},\boldsymbol{y}^{*}(\boldsymbol{x})\right).

Furthermore, $\boldsymbol{y}^{*}(\boldsymbol{x})$ and $\nabla\Phi(\boldsymbol{x})$ are Lipschitz continuous, i.e., for any $\boldsymbol{x}_{1}$ , $\boldsymbol{x}_{2}\in\mathbb{R}^{p}$ , we have

	$\displaystyle\\|\boldsymbol{y}^{}(\boldsymbol{x}_{1})-\boldsymbol{y}^{}(\boldsymbol{x}_{2})\\|\leq\frac{L}{\mu}\\|\boldsymbol{x}_{1}-\boldsymbol{x}_{2}\\|,$
	$\displaystyle\\|\nabla\Phi(\boldsymbol{x}_{1})-\nabla\Phi(\boldsymbol{x}_{2})\\|\leq L_{\Phi}\\|\boldsymbol{x}_{1}-\boldsymbol{x}_{2}\\|,$

where $L_{\Phi}=L+(2L^{2}+\tau M^{2})/{\mu}+(\rho LM+L^{3}+\tau ML)/{\mu^{2}}+(\rho L^{2}M)/{\mu^{3}}$ .

The proof of Proposition 1 is in the supplementary material. From Assumption 1 and Proposition 1, we know that $\Phi(\boldsymbol{x})$ is nonconvex and differentiable. Therefore, we want the algorithms to find an $\epsilon$ -stationary point defined as follows.

Definition 1.

For problem (2), we call $\bar{\boldsymbol{x}}$ an $\epsilon$ -stationary point of $\Phi(\boldsymbol{x})$ if $\mathbb{E}\|\nabla\Phi(\bar{\boldsymbol{x}})\|^{2}\leq\epsilon$ .

To facilitate our analysis of the complexity of the algorithm when it reaches the $\epsilon$ stationary point, we borrow the metrics of complexity as defined in Definition 2 of [14].

Definition 2.

For a stochastic function $F(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\xi})$ with $\boldsymbol{x}\in\mathbb{R}^{p}$ , $\boldsymbol{y}\in\mathbb{R}^{q}$ , and $\boldsymbol{\xi}$ being a random variable, and a vector $\boldsymbol{v}\in\mathbb{R}^{q}$ , let Gc $(F,\epsilon)$ be the number of partial derivatives $\nabla_{\boldsymbol{x}}F(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\xi})$ or $\nabla_{\boldsymbol{y}}F(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\xi})$ , and let Jv $(F,\epsilon)$ and Hv $(F,\epsilon)$ be the number of Jacobian vector products $\nabla_{\boldsymbol{x}}\nabla_{\boldsymbol{y}}F(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\xi})\boldsymbol{v}$ and Hessian vector products $\nabla_{\boldsymbol{y}}^{2}F(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\xi})\boldsymbol{v}$ , respectively.

To prove the convergence of the stochastic algorithms, as adopted in [7, 14, 17], we make the following assumptions on the stochastic derivatives.

Assumption 3.

For problem (2), the stochastic derivatives $\nabla F(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\xi})$ , $\nabla G(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\zeta})$ , $\nabla_{\boldsymbol{x}}\nabla_{\boldsymbol{y}}G(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\zeta})$ , and $\nabla_{\boldsymbol{y}}^{2}G(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\zeta})$ are unbiased estimates of $\nabla f(\boldsymbol{x},\boldsymbol{y})$ , $\nabla g(\boldsymbol{x},\boldsymbol{y})$ , $\nabla_{\boldsymbol{x}}\nabla_{\boldsymbol{y}}g(\boldsymbol{x},\boldsymbol{y})$ , and $\nabla_{\boldsymbol{y}}^{2}g(\boldsymbol{x},\boldsymbol{y})$ , respectively. Furthermore, their variances satisfy

	$\displaystyle\mathbb{E}\\|\nabla F(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\xi})-\nabla f(\boldsymbol{x},\boldsymbol{y})\\|^{2}\leq\sigma_{f}^{2},$
	$\displaystyle\mathbb{E}\\|\nabla G(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\zeta})-\nabla g(\boldsymbol{x},\boldsymbol{y})\\|^{2}\leq\sigma_{g}^{2},$
	$\displaystyle\mathbb{E}\\|\nabla_{\boldsymbol{x}}\nabla_{\boldsymbol{y}}G(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\zeta})-\nabla_{\boldsymbol{x}}\nabla_{\boldsymbol{y}}g(\boldsymbol{x},\boldsymbol{y})\\|^{2}\leq\sigma_{g,1}^{2},$
	$\displaystyle\mathbb{E}\\|\nabla_{\boldsymbol{y}}^{2}G(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\zeta})-\nabla_{\boldsymbol{y}}^{2}g(\boldsymbol{x},\boldsymbol{y})\\|^{2}\leq\sigma_{g,2}^{2},$

where $\sigma_{f}^{2}$ , $\sigma_{g}^{2}$ , $\sigma_{g,1}^{2}$ , and $\sigma_{g,2}^{2}$ are positive constants.

II-B Introduction of Some Hypergradient Estimation Methods

In the paper, we study the hypergradient estimation methods: Stochastic BP, Stochastic NS, and SGD-based Estimation, where Stochastic NS is adopted to estimate the hypergradient for all the existing simple SGD type algorithms for solving problem (2), and Stochastic BP and SGD-based Estimation are developed in this paper for problem (2). Notice that Stochastic NS is actually a stochastic version of NS. In the similar way, we obtain Stochastic BP by introducing stochastic estimation into BP, and based on AID, we obtain SGD-based Estimation by solving linear system with SGD.

In the following, we introduce how each of these hypergradient estimation methods estimates the hypergradient in a general simple SGD type algorithm shown in Algorithm 1, with problem (2) satisfying Assumptions 1, 2, and 3.

As shown in Algorithm 1, given $\boldsymbol{x}^{k}$ , first, an approximation $\boldsymbol{y}^{k,T}$ to the optimal solution $\boldsymbol{y}^{*}(\boldsymbol{x}^{k})$ of the LL problem is obtained through the SGD updates in line 5 of Algorithm 1. Based on the output $\boldsymbol{y}^{k,T}$ , these hypergradient estimation methods estimate the hypergradient as follows.

Algorithm 1 A General Simple SGD Type Algorithm for Problem (2)

K

T

J

, stepsizes

\alpha

and

\beta

0<\eta<1/L

, initialization

\boldsymbol{x^{0}}\in\mathbb{R}^{p}

1: for

k=0,\ldots,K-1

2: Initialize

\boldsymbol{y}^{k,0}\in\mathbb{R}^{q}

3: for

t=0,\ldots,T-1

4: Draw a sample batch

\mathcal{S}_{t}

from the distribution of

\boldsymbol{\zeta}

5: Update

\boldsymbol{y}^{k,t+1}=\boldsymbol{y}^{k,t}-\beta\nabla_{\boldsymbol{y}}G(\boldsymbol{x}^{k},\boldsymbol{y}^{k,t};\mathcal{S}_{t})

6: end for

7: Hypergradient estimation via one of the following :

8: Stochastic BP: automatic differentiation through the iterations in line 5 of Algorithm 1, and compute

\boldsymbol{h}_{f}^{k}

via (6)

9: Stochastic NS:

10: 1) solve (8) by approximating

[\nabla_{\boldsymbol{y}}^{2}g(\boldsymbol{x}^{k},\boldsymbol{y}^{k,T})]^{-1}

with the first

J

terms of the Neumann series in (9), and obtain an approximate solution

\boldsymbol{v}^{k,J}

in (10)

11: 2) compute

\boldsymbol{h}_{f}^{k}

via (11)

12: SGD-based Estimation:

13: 1) initialize

\boldsymbol{v}^{k,0}\in\mathbb{R}^{q}

14: 2) solve (12) with

J

steps of SGD, and obtain an approximate solution

\boldsymbol{v}^{k,J}

, as shown in (13)

15: 3) compute

\boldsymbol{h}_{f}^{k}

via (14)

16: Update

\boldsymbol{x}^{k+1}=\boldsymbol{x}^{k}-\alpha\boldsymbol{h}_{f}^{k}

17: end for

Stochastic BP. Notice that $\Phi(\boldsymbol{x})=f(\boldsymbol{x},\boldsymbol{y}^{*}(\boldsymbol{x}))$ , and that

	$\displaystyle\nabla\Phi(\boldsymbol{x}^{k})$	$\displaystyle=\nabla_{\boldsymbol{x}}f(\boldsymbol{x}^{k},\boldsymbol{y}^{*}(\boldsymbol{x}^{k}))$		(4)
		$\displaystyle\qquad+(\nabla\boldsymbol{y}^{}(\boldsymbol{x}^{k}))^{\top}\nabla_{\boldsymbol{y}}f(\boldsymbol{x}^{k},\boldsymbol{y}^{}(\boldsymbol{x}^{k})).$

For $\nabla\Phi(\boldsymbol{x}^{k})$ in the form of (4), Stochastic BP estimates the hypergradient by automatic differentiation. Specifically, since $\boldsymbol{y}^{k,T}$ has a dependence on $\boldsymbol{x}^{k}$ through the SGD updates in line 5 of Algorithm 1, Stochastic BP computes $(\nabla\boldsymbol{y}^{k,T}(\boldsymbol{x}^{k}))^{\top}$ by automatic differentiation along these SGD updates, and $(\nabla\boldsymbol{y}^{k,T}(\boldsymbol{x}^{k}))^{\top}$ can be obtained as follows:

(\nabla\boldsymbol{y}^{k,T}(\boldsymbol{x}^{k}))^{\top}=-\beta\sum\limits_{t=0}^{T-1}\nabla_{\boldsymbol{x}}\nabla_{\boldsymbol{y}}G(\boldsymbol{x}^{k},\boldsymbol{y}^{k,t};\mathcal{S}_{t})\mathbf{A}_{T-1-t}

(5)

where $\mathbf{A}_{T-1-t}=\prod_{i=t+1}^{T-1}(\mathbf{I}-\beta\nabla_{\boldsymbol{y}}^{2}G(\boldsymbol{x}^{k},\boldsymbol{y}^{k,i};\mathcal{S}_{i}))$ with $\mathbf{I}$ being the identity matrix, and $\mathbf{A}_{0}=\mathbf{I}$ (for the derivation, see Appendix D.2 in [14]). Then,

	$\displaystyle\boldsymbol{h}_{f}^{k}$	$\displaystyle=\nabla_{\boldsymbol{x}}F(\boldsymbol{x}^{k},\boldsymbol{y}^{k,T};\mathcal{D}_{F})$		(6)
		$\displaystyle\qquad+(\nabla\boldsymbol{y}^{k,T}(\boldsymbol{x}^{k}))^{\top}\nabla_{\boldsymbol{y}}F(\boldsymbol{x}^{k},\boldsymbol{y}^{k,T};\mathcal{D}_{F})$

is used to estimate the hypergradient $\nabla\Phi(\boldsymbol{x}^{k})$ in (4), where $\mathcal{D}_{F}$ is a sample set sampled from the distribution of $\boldsymbol{\xi}$ and $(\nabla\boldsymbol{y}^{k,T}(\boldsymbol{x}^{k}))^{\top}$ is given in (5).

In fact, Stochastic BP is a variant of BP, and $\boldsymbol{h}_{f}^{k}$ in (6) can be directly obtained by introducing stochastic estimation into the hypergradient approximation formula of BP in Proposition 2 in [14].

Stochastic NS. From Proposition 1, we have

	$\displaystyle\nabla\Phi(\boldsymbol{x}^{k})$	$\displaystyle=\nabla_{\boldsymbol{x}}f(\boldsymbol{x}^{k},\boldsymbol{y}^{*}(\boldsymbol{x}^{k}))$		(7)
		$\displaystyle\qquad-\nabla_{\boldsymbol{x}}\nabla_{\boldsymbol{y}}g(\boldsymbol{x}^{k},\boldsymbol{y}^{}(\boldsymbol{x}^{k}))\boldsymbol{v}(\boldsymbol{x}^{k},\boldsymbol{y}^{}(\boldsymbol{x}^{k}))$

with $\boldsymbol{v}(\boldsymbol{x}^{k},\boldsymbol{y}^{*}(\boldsymbol{x}^{k}))=[\nabla_{\boldsymbol{y}}^{2}g(\boldsymbol{x}^{k},\boldsymbol{y}^{*}(\boldsymbol{x}^{k}))]^{-1}\nabla_{\boldsymbol{y}}f(\boldsymbol{x}^{k},\boldsymbol{y}^{*}(\boldsymbol{x}^{k}))$ . For $\nabla\Phi(\boldsymbol{x}^{k})$ in the form of (7), Stochastic NS estimates the hypergradient by solving the following linear system:

\nabla_{\boldsymbol{y}}^{2}g(\boldsymbol{x}^{k},\boldsymbol{y}^{k,T})\boldsymbol{v}=\nabla_{\boldsymbol{y}}f(\boldsymbol{x}^{k},\boldsymbol{y}^{k,T}).

(8)

By approximating Hessian inverse $[\nabla_{\boldsymbol{y}}^{2}g(\boldsymbol{x}^{k},\boldsymbol{y}^{k,T})]^{-1}$ with the first $J$ terms of the following Neumann series:

\sum_{j=0}^{\infty}\eta(\mathbf{I}-\eta\nabla_{\boldsymbol{y}}^{2}g(\boldsymbol{x}^{k},\boldsymbol{y}^{k,T}))^{j}

(9)

and computing $\nabla_{\boldsymbol{y}}F(\boldsymbol{x}^{k},\boldsymbol{y}^{k,T};\mathcal{D}_{F})$ over a sample set $\mathcal{D}_{F}$ sampled from the distribution of $\boldsymbol{\xi}$ , an approximate solution $\boldsymbol{v}^{k,J}$ of (8) is obtained, and takes the form of

\boldsymbol{v}^{k,J}=\eta\sum\limits_{j=0}^{J-1}\prod_{i=J-j}^{J-1}(\mathbf{I}-\eta\nabla_{\boldsymbol{y}}^{2}G(\boldsymbol{x}^{k},\boldsymbol{y}^{k,T};\mathcal{B}_{i}))\boldsymbol{v}_{0}

(10)

where $\prod_{i=J}^{J-1}(\cdot)=\mathbf{I}$ , $\boldsymbol{v}_{0}=\nabla_{\boldsymbol{y}}F(\boldsymbol{x}^{k},\boldsymbol{y}^{k,T};\mathcal{D}_{F})$ , and sample sets $\mathcal{B}_{i},i=1,\ldots,J-1$ are sampled from the distribution of $\boldsymbol{\zeta}$ . Then,

\boldsymbol{h}_{f}^{k}=\nabla_{\boldsymbol{x}}F(\boldsymbol{x}^{k},\boldsymbol{y}^{k,T};\mathcal{D}_{F})-\nabla_{\boldsymbol{x}}\nabla_{\boldsymbol{y}}G(\boldsymbol{x}^{k},\boldsymbol{y}^{k,T};\mathcal{D}_{G})\boldsymbol{v}^{k,J}

(11)

is used to estimate the hypergradient $\nabla\Phi(\boldsymbol{x}^{k})$ in (7), where $\mathcal{D}_{F}$ is the aforementioned sample set, $\mathcal{D}_{G}$ is a sample set sampled from the distribution of $\boldsymbol{\zeta}$ , and $\boldsymbol{v}^{k,J}$ is given in (10). Note that sample sets $\{\mathcal{D}_{F},\mathcal{D}_{G},\mathcal{B}_{i}(i=1,\ldots,J-1)\}$ are mutually independent.

In fact, $\boldsymbol{h}_{f}^{k}$ in (11) can be directly obtained by introducing stochastic estimation into the hypergradient approximation formula of NS in eq.(3) in [22], which is the reason we call this hypergradient estimation method Stochastic NS. Furthermore, for more details about Stochastic NS, see subsection 2.2 in [14].

SGD-based Estimation. SGD-based Estimation estimates the hypergradient in the way similar to Stochastic NS except that the approximate solution to (8) is obtained by using SGD. Notice that the solution to (8) is equivalent to the optimal solution of the following optimization problem:

\min\limits_{\boldsymbol{v}}\frac{1}{2}\boldsymbol{v}^{\top}\nabla_{\boldsymbol{y}}^{2}g(\boldsymbol{x}^{k},\boldsymbol{y}^{k,T})\boldsymbol{v}-\boldsymbol{v}^{\top}\nabla_{\boldsymbol{y}}f(\boldsymbol{x}^{k},\boldsymbol{y}^{k,T}).

(12)

SGD-based Estimation obtains an approximate solution $\boldsymbol{v}^{k,J}$ to (12) by solving (12) with $J$ steps of SGD as follows:

	$\displaystyle\boldsymbol{v}^{k,j+1}$	$\displaystyle=(\mathbf{I}-\eta\nabla_{\boldsymbol{y}}^{2}G(\boldsymbol{x}^{k},\boldsymbol{y}^{k,T};\mathcal{B}_{j}))\boldsymbol{v}^{k,j}$		(13)
		$\displaystyle\qquad+\eta\nabla_{\boldsymbol{y}}F(\boldsymbol{x}^{k},\boldsymbol{y}^{k,T};\mathcal{D}_{F,j}),\quad j=0,\ldots,J-1$

where $\eta$ is the step size, and for each $j\in\{0,\ldots,J-1\}$ , $\mathcal{D}_{F,j}$ and $\mathcal{B}_{j}$ are the sample sets sampled from the distributions of $\boldsymbol{\xi}$ and $\boldsymbol{\zeta}$ , respectively. Then,

\boldsymbol{h}_{f}^{k}=\nabla_{\boldsymbol{x}}F(\boldsymbol{x}^{k},\boldsymbol{y}^{k,T};\mathcal{D}_{F})-\nabla_{\boldsymbol{x}}\nabla_{\boldsymbol{y}}G(\boldsymbol{x}^{k},\boldsymbol{y}^{k,T};\mathcal{D}_{G})\boldsymbol{v}^{k,J}

(14)

is used to estimate the hypergradient, where $\mathcal{D}_{F}$ and $\mathcal{D}_{G}$ are the sample sets sampled from the distributions of $\boldsymbol{\xi}$ and $\boldsymbol{\zeta}$ , respectively.

II-C Analysis of Hypergradient Estimation Methods

Define $\mathcal{F}_{k,1}^{{}^{\prime}}=\sigma\{\boldsymbol{x}^{0},\ldots,\boldsymbol{y}^{0,T},\ldots,\boldsymbol{x}^{k},\ldots,\boldsymbol{y}^{k,T}\}$ , $\mathcal{F}_{k,2}^{{}^{\prime}}=\sigma\{\boldsymbol{x}^{0},\ldots,\boldsymbol{y}^{0,T},\ldots,\boldsymbol{v}^{0,J},\ldots,\boldsymbol{x}^{k},\ldots,\boldsymbol{y}^{k,T}\}$ , where $\sigma\{\cdot\}$ denotes the $\sigma$ -algebra generated by random variables.

Then, similar to the discussion in the proof of Lemma 3 in the supplementary material, from the iteration in line 16 of Algorithm 1, we can obtain formula (19) in the supplementary material in which $\bar{\boldsymbol{h}}_{f}^{k}=\mathbb{E}[\boldsymbol{h}_{f}^{k}|\mathcal{F}_{k,1}^{{}^{\prime}}]$ for Stochastic BP and Stochastic NS, and $\bar{\boldsymbol{h}}_{f}^{k}=\mathbb{E}[\boldsymbol{h}_{f}^{k}|\mathcal{F}_{k,2}^{{}^{\prime}}]$ for SGD-based Estimation. By shifting the terms, we have

	$\displaystyle\mathbb{E}\\|\nabla\Phi(\boldsymbol{x}^{k})\\|^{2}$
	$\displaystyle\leq\frac{2}{\alpha}\big{(}\mathbb{E}[\Phi(\boldsymbol{x}^{k})]-\mathbb{E}[\Phi(\boldsymbol{x}^{k+1})]\big{)}-(1-\alpha{L_{\Phi}})\mathbb{E}\\|\bar{\boldsymbol{h}}_{f}^{k}\\|^{2}$
	$\displaystyle\qquad+\mathbb{E}\\|\nabla\Phi(\boldsymbol{x}^{k})-\bar{\boldsymbol{h}}_{f}^{k}\\|^{2}+\alpha{L_{\Phi}}\mathbb{E}\\|\boldsymbol{h}_{f}^{k}-\bar{\boldsymbol{h}}_{f}^{k}\\|^{2}$

where $\alpha$ is in Algorithm 1, $L_{\Phi}$ is in Proposition 1.

The upper bound of $\mathbb{E}\|\nabla\Phi(\boldsymbol{x}^{k})\|^{2}$ involves the term $\|\nabla\Phi(\boldsymbol{x}^{k})-\bar{\boldsymbol{h}}_{f}^{k}\|$ . In the following, the upper bound of $\|\nabla\Phi(\boldsymbol{x}^{k})-\bar{\boldsymbol{h}}_{f}^{k}\|$ is provided. For the proof, see Lemma 10 in the supplementary material.

Proposition 2.

Apply Algorithm 1 to solve problem (2). Suppose Assumptions 1, 2, and 3 hold. let $\mathbf{H}_{0}=\mathbf{0}$ , $\boldsymbol{v}_{0}=\boldsymbol{v}^{k,0}$ , and define $\mathbf{H}^{*}=\nabla_{\boldsymbol{x}}\nabla_{\boldsymbol{y}}g(\boldsymbol{x}^{k},\boldsymbol{y}^{*}(\boldsymbol{x}^{k}))(\nabla_{\boldsymbol{y}}^{2}g(\boldsymbol{x}^{k},\boldsymbol{y}^{*}(\boldsymbol{x}^{k})))^{-1}$ , $\boldsymbol{v}^{*}=\nabla_{\boldsymbol{y}}^{2}g(\boldsymbol{x}^{k},\boldsymbol{y}^{*}(\boldsymbol{x}^{k}))^{-1}\nabla_{\boldsymbol{y}}f(\boldsymbol{x}^{k},\boldsymbol{y}^{*}(\boldsymbol{x}^{k}))$ , where $\boldsymbol{v}^{k,0}$ is the initial value to solve problem (12)(see line 13 in Algorithm 1), and $\boldsymbol{y}^{*}(\boldsymbol{x}^{k})$ is the optimal solution to the LL problem of problem (2) with $\boldsymbol{x}=\boldsymbol{x}^{k}$ . For $\mu$ , $L$ in Assumptions 1, 2 and $T$ , $J$ , $\beta$ , $\eta$ in Algorithm 1, let $\mu<1$ , $L<1$ , $T$ , $J$ be any positive integers, and $\beta=\eta$ . Then,

•

For Stochastic BP, for $\bar{\boldsymbol{h}}_{f}^{k}:=\mathbb{E}[\boldsymbol{h}_{f}^{k}|\mathcal{F}_{k,1}^{{}^{\prime}}]$ with $\boldsymbol{h}_{f}^{k}$ in (6), we have

	$\displaystyle\\|\nabla\Phi(\boldsymbol{x}^{k})-\bar{\boldsymbol{h}}_{f}^{k}\\|$
	$\displaystyle\leq L\bigg{(}1+\frac{L}{\mu}\bigg{)}\\|\boldsymbol{y}^{k,T}-\boldsymbol{y}^{}(\boldsymbol{x}^{k})\\|+M(1-\eta\mu)^{T}\\|\mathbf{H}^{}\\|$
	$\displaystyle+M\eta\bigg{(}\frac{L}{\mu}\rho+\tau\bigg{)}\sum\limits_{t=0}^{T-1}(1-\eta\mu)^{t}\\|\boldsymbol{y}^{k,T-1-t}-\boldsymbol{y}^{*}(\boldsymbol{x}^{k})\\|;$

•

For Stochastic NS, for $\bar{\boldsymbol{h}}_{f}^{k}:=\mathbb{E}[\boldsymbol{h}_{f}^{k}|\mathcal{F}_{k,1}^{{}^{\prime}}]$ with $\boldsymbol{h}_{f}^{k}$ in (11), we have

	$\displaystyle\\|\nabla\Phi(\boldsymbol{x}^{k})-\bar{\boldsymbol{h}}_{f}^{k}\\|$
	$\displaystyle\leq L\bigg{(}1+\frac{L}{\mu}\bigg{)}\\|\boldsymbol{y}^{k,T}-\boldsymbol{y}^{}(\boldsymbol{x}^{k})\\|+M(1-\eta\mu)^{J}\\|\mathbf{H}^{}\\|$
	$\displaystyle+M\eta\bigg{(}\frac{L}{\mu}\rho+\tau\bigg{)}\sum\limits_{t=0}^{J-1}(1-\eta\mu)^{t}\\|\boldsymbol{y}^{k,T}-\boldsymbol{y}^{*}(\boldsymbol{x}^{k})\\|;$

•

For SGD-based Estimation, for $\bar{\boldsymbol{h}}_{f}^{k}:=\mathbb{E}[\boldsymbol{h}_{f}^{k}|\mathcal{F}_{k,2}^{{}^{\prime}}]$ with $\boldsymbol{h}_{f}^{k}$ in (14), we have

	$\displaystyle\\|\nabla\Phi(\boldsymbol{x}^{k})-\bar{\boldsymbol{h}}_{f}^{k}\\|$
	$\displaystyle\leq\bigg{(}L+\frac{M}{\mu}\tau\bigg{)}\\|\boldsymbol{y}^{k,T}-\boldsymbol{y}^{*}(\boldsymbol{x}^{k})\\|$
	$\displaystyle+L(1-\eta\mu)^{J}\\|\boldsymbol{v}^{k,0}-\boldsymbol{v}^{*}\\|$
	$\displaystyle+L\eta\bigg{(}\frac{M}{\mu}\rho+L\bigg{)}\sum\limits_{t=0}^{J-1}(1-\eta\mu)^{t}\\|\boldsymbol{y}^{k,T}-\boldsymbol{y}^{*}(\boldsymbol{x}^{k})\\|;$

$\nabla\Phi(\boldsymbol{x}^{k})$ is in (3) with $\boldsymbol{x}=\boldsymbol{x}^{k}$ , $\mu$ , $L$ , $M$ , $\tau$ , $\rho$ are given in Assumptions 1, 2.

From Proposition 2, it is easy to find that there are some differences for Stochastic BP, Stochastic NS, and SGD-based Estimation regarding the upper bound of $\|\nabla\Phi(\boldsymbol{x}^{k})-\bar{\boldsymbol{h}}_{f}^{k}\|$ . For Stochastic BP and Stochastic NS, the upper bound for $\|\nabla\Phi(\boldsymbol{x}^{k})-\bar{\boldsymbol{h}}_{f}^{k}\|$ involves $M(1-\eta\mu)^{T}\|\mathbf{H}^{*}\|$ and $M(1-\eta\mu)^{J}\|\mathbf{H}^{*}\|$ , respectively. While for SGD-based Estimation, the term $L(1-\eta\mu)^{J}\|\boldsymbol{v}^{k,0}-\boldsymbol{v}^{*}\|$ is involved in the upper bound of $\|\nabla\Phi(\boldsymbol{x}^{k})-\bar{\boldsymbol{h}}_{f}^{k}\|$ .

Thus, to ensure that $\|\nabla\Phi(\boldsymbol{x}^{k})-\bar{\boldsymbol{h}}_{f}^{k}\|$ is a sufficiently small number, a sufficiently large $T$ for Stochastic BP and a sufficiently large $J$ for Stochastic NS may be required. In contrast, for SGD-based Estimation, we notice that $\boldsymbol{v}^{k,0}$ in the term $L(1-\eta\mu)^{J}\|\boldsymbol{v}^{k,0}-\boldsymbol{v}^{*}\|$ is the initial value to solve problem (12). Apply warm start initialization strategy to $\boldsymbol{v}^{k,0}$ (i.e., $\boldsymbol{v}^{k,0}=\boldsymbol{v}^{k-1,J}$ , where $\boldsymbol{v}^{k-1,J}$ is the output of the $J$ -th iteration to solve problem (12) given $\boldsymbol{x}^{k-1}$ ) can allow us to track the errors, such as $\|\boldsymbol{v}^{k-1,0}-\boldsymbol{v}^{*}\|$ , …, $\|\boldsymbol{v}^{0,0}-\boldsymbol{v}^{*}\|$ , in the preceeding loops. Thus, a smaller number $J$ may be able to ensure $\|\nabla\Phi(\boldsymbol{x}^{k})-\bar{\boldsymbol{h}}_{f}^{k}\|$ to be a sufficiently small number, and there is no restriction on $T$ . Notice that warm start initialization strategy is used in stocBiO[14] to initialize the initial values for solving the LL problem, and the improved computational complexity is obtained.

The above discussion inspires us to think that SGD-based Estimation may be able to allow $\mathbb{E}\|\nabla\Phi(\boldsymbol{x}^{k})\|^{2}$ to be sufficiently small under weaker restrictions on $T$ than Stochastic BP and $J$ than Stochastic NS. Thus, SGD-based Estimation may be able to make Algorithm 1 converge faster than Stochastic BP and Stochastic NS.

Next, we perform experiments on a synthetic bilevel optimization problem, a special case of problem (2), to evaluate the performance of Algorithm 1 under different hypergradient estimation methods.

Synthetic Bilevel Optimization Problem. We consider the following synthetic bilevel optimization problem:

\min\limits_{\boldsymbol{x}\in\mathbb{R}^{p}}\Phi(\boldsymbol{x}):=f(\boldsymbol{x},\boldsymbol{y}^{*}(\boldsymbol{x}))\quad s.t.~{}\boldsymbol{y}^{*}(\boldsymbol{x})={\underset{\boldsymbol{y}\in\mathbb{R}^{p}}{\operatorname{arg\,min}}\,g(\boldsymbol{x},\boldsymbol{y})}

(15)

with

	$\displaystyle f(\boldsymbol{x},\boldsymbol{y})=\frac{1}{\|\mathcal{D}_{\text{val}}\|}\sum\limits_{(\boldsymbol{u}_{i},v_{i})\in\mathcal{D}_{\text{val}}}\frac{1}{2}({\boldsymbol{y}}^{\top}\boldsymbol{u}_{i}-v_{i})^{2}+\\|\boldsymbol{x}\\|^{3},$
	$\displaystyle g(\boldsymbol{x},\boldsymbol{y})=\frac{1}{\|\mathcal{D}_{\text{tr}}\|}\sum\limits_{(\boldsymbol{u}_{i},v_{i})\in\mathcal{D}_{\text{tr}}}\frac{1}{2}(\boldsymbol{y}^{\top}\boldsymbol{u}_{i}-v_{i})^{2}+\frac{r}{2}\\|\boldsymbol{y}-\boldsymbol{x}\\|^{2}$

where $r=0.5$ , and $\mathcal{D}_{\text{tr}}$ and $\mathcal{D}_{\text{val}}$ are the data sets constructed as follows: given $\boldsymbol{w}_{0}\in\mathbb{R}^{p}$ , first randomly sample $20000$ data points $\boldsymbol{e}_{i}\in\mathbb{R}^{p-1}$ , $i=1,\ldots,20000$ from normal distribution with a mean of 0 and a variance of 0.01; then, for each $i\in\{1,\ldots,20000\}$ , we set $\boldsymbol{u}_{i}=(\boldsymbol{e}_{i},1)$ , and construct $v_{i}\in\mathbb{R}$ by $v_{i}=\boldsymbol{w}_{0}^{\top}\boldsymbol{u}_{i}+\sigma_{i}$ , where $\sigma_{i}$ is the Gaussian noise with mean 0 and variance 1; finally, dividing the dataset $\{(\boldsymbol{u}_{i},v_{i})\}_{i=1}^{20000}$ equally into two parts, we obtain training set $\mathcal{D}_{\text{tr}}$ and validation set $\mathcal{D}_{\text{val}}$ , with the training set size $|\mathcal{D}_{\text{tr}}|$ and validation set size $|\mathcal{D}_{\text{val}}|$ being 10000, respectively.

It is easy to verify that $g(\boldsymbol{x},\boldsymbol{y})$ is strongly convex w.r.t. $\boldsymbol{y}$ . Moreover, for each $\boldsymbol{x}$ , the minimum solution of the LL problem is

\boldsymbol{y}^{*}(\boldsymbol{x})=(\mathbf{A}_{\text{tr}}+r\mathbf{I})^{-1}(\boldsymbol{b_{\text{tr}}}+r\boldsymbol{x})

and the hypergradient is

\nabla\Phi(\boldsymbol{x})=r(\mathbf{A}_{\text{tr}}+r\mathbf{I})^{-1}(\mathbf{A}_{\text{val}}\boldsymbol{y}^{*}(\boldsymbol{x})-\boldsymbol{b_{\text{val}}})+3\|\boldsymbol{x}\|\boldsymbol{x}

where

	$\displaystyle\mathbf{A}_{\text{val}}=\frac{1}{\|\mathcal{D}_{\text{val}}\|}\sum\limits_{(\boldsymbol{u}_{i},v_{i})\in\mathcal{D}_{\text{val}}}\boldsymbol{u}_{i}(\boldsymbol{u}_{i})^{\top},$
	$\displaystyle\boldsymbol{b_{\text{val}}}=\frac{1}{\|\mathcal{D}_{\text{val}}\|}\sum\limits_{(\boldsymbol{u}_{i},v_{i})\in\mathcal{D}_{\text{val}}}v_{i}\boldsymbol{u}_{i},$

and $\mathbf{A}_{\text{tr}}$ , $\boldsymbol{b_{\text{tr}}}$ can be obtained by replacing $\mathcal{D}_{\text{val}}$ and $|\mathcal{D}_{\text{val}}|$ in $\mathbf{A}_{\text{val}}$ and $\boldsymbol{b_{\text{val}}}$ with $\mathcal{D}_{\text{tr}}$ and $|\mathcal{D}_{\text{tr}}|$ .

In the following, we use Algorithm 1 to solve problem (15), and compare the performance of Algorithm 1 under different hypergradient estimation methods. To be specific, the hypergradient estimation methods that we consider are Stochastic BP, Stochastic NS, and SGD-based Estimation, where for SGD-based Estimation, we use two initialization strategies to initialize $\boldsymbol{v}^{k,0}$ in line 13 of Algorithm 1, i.e., non-warm start initialization strategy and warm start initialization strategy. For non-warm start initialization strategy, we set $\boldsymbol{v}^{k,0}$ , $k=0,\ldots,K-1$ to be zero vectors. For warm-start initialization strategy, we set $\boldsymbol{v}^{k,0}=\boldsymbol{v}^{k-1,J}$ , $k=1,\ldots,K-1$ , where for each $k\in\{1,\ldots,K-1\}$ , $\boldsymbol{v}^{k-1,J}$ is obtained by the iteration steps starting from $\boldsymbol{v}^{k-1,0}$ in (13). The experimental details are in the supplementary material.

To simplify the narrative, we refer to Algorithm 1 as Stochastic BP(resp. Stochastic NS) if Algorithm 1 estimates the hypergradient via Stochastic BP(resp. Stochastic NS). Similarly, if Algorithm 1 estimates the hypergradient via SGD-based Estimation and using the non-warm start initialization strategy(resp. warm start initialization strategy) to initialize $\boldsymbol{v}^{k,0}$ , we refer to Algorithm 1 as SGD_W_Start_False(resp. SGD_W_Start_True). Furthermore, we use P_ $T$ _ $i$ (P_ $J$ _ $i$ ) to denote that the iteration step $T$ ( $J$ ) is set to be $i$ for P.

In the experiments, data sets $\mathcal{D}_{\text{tr}}$ and $\mathcal{D}_{\text{val}}$ are constructed by setting $\boldsymbol{w}_{0}=(2,5,7)\in\mathbb{R}^{3}$ . From Fig. 1LABEL:sub@fig_first_case, we observe that to converge to a stationary point with $\|\nabla\Phi(\boldsymbol{x})\|=0$ , $T=1$ is sufficient for SGD_W_Start_True, and a larger $T$ ( $T=20$ ) is required for Stochastic BP. Furthermore, SGD_W_Start_True always converges faster than Stochastic BP for $T\in\{1,5,20\}$ . From Fig. 1LABEL:sub@fig_second_case, we observe that to converge to a stationary point with $\|\nabla\Phi(\boldsymbol{x})\|=0$ , $J=1$ is sufficient for SGD_W_Start_True, and a larger $J$ ( $J=20$ ) is required for Stochastic NS. Furthermore, SGD_W_Start_True with $J=1$ has a faster convergence speed than Stochastic NS with $J=20$ . From Fig. 1LABEL:sub@fig_third_case, we observe that to converge to a stationary point with $\|\nabla\Phi(\boldsymbol{x})\|=0$ , $J=1$ is sufficient for SGD_W_Start_True, and a larger $J$ ( $J=20$ ) is required for SGD_W_Start_False. Furthermore, SGD_W_Start_True with $J=1$ has a faster convergence speed than SGD_W_Start_False with $J=20$ .

In conclusion, the above experimental results show that compared with Stochastic BP and Stochastic NS, Algorithm 1 can obtain a faster convergence speed by using SGD-based Estimation to estimate the hypergradient and using warm start initialization strategy for SGD-based Estimation.

III Algorithm

In Section II, we find that compared with Stochastic BP and Stochastic NS, SGD-based Estimation can make Algorithm 1 achieve a faster convergence speed. Notice that all the existing simple SGD type algorithms for solving problem (2) use Stochastic NS to estimate the hypergradient. It is natural to think that using SGD-based Estimation to estimate the hypergradient may allow us to obtain an algorithm with lower computational complexity for problem (2).

Algorithm 2 A New Simple SGD Type Algorithm SSGD for Problem (2)

K\geq 1

T\geq 1

J\geq 1

, stepsizes

\alpha

\beta

, and

\eta

, initializations

\boldsymbol{x}^{0}\in\mathbb{R}^{p}

\boldsymbol{y}^{0}\in\mathbb{R}^{q}

\boldsymbol{v}^{0}\in\mathbb{R}^{q}

1: for

k=0,\ldots,K-1

\boldsymbol{x}^{k+1}=\boldsymbol{x}^{k}-\alpha\boldsymbol{h}_{f}^{k}

3: Set

\boldsymbol{y}^{k+1,0}=\boldsymbol{y}^{k}

4: for

t=0,\ldots,T-1

5: Draw a sample batch

\mathcal{S}_{t}

\boldsymbol{y}^{k+1,t+1}=\boldsymbol{y}^{k+1,t}-\beta\nabla_{\boldsymbol{y}}G(\boldsymbol{x}^{k+1},\boldsymbol{y}^{k+1,t};\mathcal{S}_{t})

7: end for

8: Set

\boldsymbol{y}^{k+1}=\boldsymbol{y}^{k+1,T}

\boldsymbol{v}^{k+1,0}=\boldsymbol{v}^{k}

9: for

j=0,\ldots,J-1

10: Draw sample batches

\mathcal{B}_{j}

\mathcal{D}_{F,j}

11:

\boldsymbol{v}^{k+1,j+1}=(\mathbf{I}-\eta\nabla_{\boldsymbol{y}}^{2}G(\boldsymbol{x}^{k+1},\boldsymbol{y}^{k+1};\mathcal{B}_{j}))\boldsymbol{v}^{k+1,j}

\qquad+\eta\nabla_{\boldsymbol{y}}F(\boldsymbol{x}^{k+1},\boldsymbol{y}^{k+1};\mathcal{D}_{F,j})

12: end for

13: Set

\boldsymbol{v}^{k+1}=\boldsymbol{v}^{k+1,J}

14: Draw sample batches

\mathcal{D}_{F}

\mathcal{D}_{G}

15: Compute

\boldsymbol{h}_{f}^{k+1}

via

(\ref{eq15})

16: end for

As shown in Algorithm 2, by using warm start initialization strategy, a new simple SGD type algorithm SSGD based on SGD-based Estimation is proposed, in which warm start is not only used for SGD-based Estimation, but also for solving LL problem.

In Algorithm 2, given $\boldsymbol{x}^{k+1}$ , SSGD runs $T$ steps of SGD to obtain an approximation $\boldsymbol{y}^{k+1,T}$ to the optimal solution $\boldsymbol{y}^{*}(\boldsymbol{x}^{k+1})$ of the LL problem(see line 6 in Algorithm 2). Note that the initial value $\boldsymbol{y}^{k+1,0}$ is initialized by using warm start and is set to be $\boldsymbol{y}^{k}$ , where $\boldsymbol{y}^{k}$ is equal to $\boldsymbol{y}^{k,T}$ , which is the output of the $T$ -th iteration starting from $\boldsymbol{y}^{k,0}$ in line 6 of Algorithm 2.

After setting $\boldsymbol{y}^{k+1}=\boldsymbol{y}^{k+1,T}$ , the hypergradient is estimated via SGD-based Estimation. Specifically, SSGD first obtains an approximation $\boldsymbol{v}^{k+1,J}$ to the optimal solution of the following optimization problem

\min\limits_{\boldsymbol{v}}\frac{1}{2}\boldsymbol{v}^{\top}\nabla_{\boldsymbol{y}}^{2}g(\boldsymbol{x}^{k+1},\boldsymbol{y}^{k+1})\boldsymbol{v}-\boldsymbol{v}^{\top}\nabla_{\boldsymbol{y}}f(\boldsymbol{x}^{k+1},\boldsymbol{y}^{k+1})

using $J$ steps of SGD starting from $\boldsymbol{v}^{k+1,0}$ (see line 11 in Algorithm 2), where we also adopt a warm start with $\boldsymbol{v}^{k+1,0}=\boldsymbol{v}^{k}$ . Then, by setting $\boldsymbol{v}^{k+1}=\boldsymbol{v}^{k+1,J}$ ,

	$\displaystyle\boldsymbol{h}_{f}^{k+1}$	$\displaystyle=\nabla_{\boldsymbol{x}}F(\boldsymbol{x}^{k+1},\boldsymbol{y}^{k+1};\mathcal{D}_{F})$		(16)
		$\displaystyle\qquad-\nabla_{\boldsymbol{x}}\nabla_{\boldsymbol{y}}G(\boldsymbol{x}^{k+1},\boldsymbol{y}^{k+1};\mathcal{D}_{G})\boldsymbol{v}^{k+1}$

is constructed as an estimate to $\nabla\Phi(\boldsymbol{x}^{k+1})$ . Note that sample sets $\{\mathcal{D}_{F},\mathcal{D}_{G},\mathcal{B}_{j}(j=0,\ldots,J-1),\mathcal{D}_{F,j}(j=0,\ldots,J-1)\}$ are mutually independent in the hypergradient estimation process.

For the sample sets in Algorithm (2), we suppose the sample sets $\mathcal{S}_{t}$ for all $t=0,\ldots,T-1$ , $\mathcal{B}_{j}$ for all $j=0,\ldots,J-1$ , and $\mathcal{D}_{G}$ have the sizes of $S$ , $D$ , and $D_{g}$ , respectively. Furthermore, we suppose sample set $\mathcal{D}_{F}$ and the sample sets $\mathcal{D}_{F,j}$ for all $j=0,\ldots,J-1$ have the same sample size $D_{f}$ .

Remark 1.

From Algorithm 2, it is easy to observe that there are three loops. However, since in the classification of the existing stochastic algorithms for solving problem (2), the loop to estimate the hypergradient, i.e., approximate the Hessian inverse, is typically not counted. To be consistent with the existing literature, when performing classification, the loop to solve linear system in lines 9-12 of Algorithm 2 is not counted. Then, it is obvious that SSGD can be either a single loop algorithm(i.e., $T=1$ ) or a double loop algorithm(i.e., $T>1$ ). Although in the analysis that follows, it can be concluded that the convergence of SSGD is guaranteed for any $T\geq 1$ . To obtain a lower computational complexity than the existing simple SGD type algorithms, $T=\mathcal{O}(\kappa)$ is required for SSGD, where $\kappa$ is the condition number. Therefore, we call SSGD a double loop algorithm, as shown in Table I.

IV Theoretical Results

In this section, we provide the convergence and complexity results of SSGD. Detailed theoretical analysis can be found in the supplementary material.

We first show that the convergence of SSGD can be guaranteed for any $T\geq 1$ and $J\geq 1$ . Furthermore, the complexity result is provided.

Theorem 1.

Apply SSGD to solve problem (2). Suppose Assumptions 1, 2, and 3 hold. Define $r_{v}={2\mu L}/{(L^{2}+\mu^{2})}=\mathcal{O}(\kappa^{-1})$ , $r_{w}={\eta\mu}/{(7(2-\eta\mu))}$ , $\rho_{2}=\mathcal{O}(\kappa^{-4})$ , $\rho_{1}={(2\rho_{2}\left((1+r_{w})C_{23}+2C_{24}\right))}/{(\rho_{y}(1+r_{v})-r_{v})}$ , and

L_{31}=\frac{L_{\Phi}}{2}+2\left((\rho_{1}+C_{24}\rho_{2})\frac{2}{r_{v}}\frac{L^{2}}{\mu^{2}}+(1+r_{w})C_{22}\rho_{2}\right),

where $\kappa={L}/{\mu}$ is the condition number, $\rho_{y}=2\beta\mu L/(\mu+L)$ , $L_{\Phi}=\mathcal{O}(\kappa^{3})$ is defined in Proposition 1, $C_{22}$ , $C_{23}$ , and $C_{24}$ are defined in Lemma 7 in the supplementary material, and $\eta$ , $\beta$ are the stepsizes. Let $\|\boldsymbol{v}^{0}\|\leq M/\mu$ , and choose stepsizes $\eta={1}/{(2L)}$ , $\beta={3}/{(2(L+\mu))}=\mathcal{O}(1)$ , and

\alpha=\min\left\{\frac{1}{2L_{31}},~{}\frac{1}{4L^{2}}\rho_{2}\eta\mu,~{}\left((1+r_{w})C_{23}+2C_{24}\right)\frac{\rho_{2}}{\bar{C}_{1}^{2}}\right\},

where $\bar{C}_{1}=L+{M\tau}/{\mu}=\mathcal{O}({\kappa})$ is defined in Lemma 4 in the supplementary material. Then, for any $T\geq 1$ and $J\geq 1$ , we have

	$\displaystyle\frac{1}{K}\sum_{k=0}^{K-1}\mathbb{E}\\|\nabla\Phi(\boldsymbol{x}^{k})\\|^{2}$
	$\displaystyle\leq\mathcal{O}\left(\frac{\kappa^{5}}{K}+\kappa^{2}\frac{\sigma_{f}^{2}}{D_{f}}+\kappa^{2}\frac{\sigma_{g,1}^{2}}{D_{g}}+\kappa^{8}\frac{1}{S}\sigma_{g}^{2}+\kappa^{4}\frac{\sigma_{g,2}^{2}}{D}\right).$

Furthermore, to achieve an $\epsilon$ -accurate stationary point, the computational complexity is $\mathcal{O}(\kappa^{13}\epsilon^{-2}\max\{T,J\})$ .

Theorem 1 shows that SSGD converges sublinearly w.r.t. the number $K$ of iterations, the batch sizes $D_{f}$ , $S$ for gradient estimation, $D_{g}$ for Jacobian estimation, and $D$ for Hessian matrix estimation. In addition, from Theorem 1, we know that when $T=1$ and $J=1$ , SSGD can achieve the lowest computational complexity $\mathcal{O}(\kappa^{13}\epsilon^{-2})$ . Compared with the existing simple SGD type algorithms listed in Table I, this computational complexity of SSGD outperforms that of BSA, TTSA, and stocBiO by an order of $\mathcal{O}(\epsilon^{-1}\log({1}/{\epsilon}))$ , $\mathcal{O}(\epsilon^{-0.5}\log({1}/{\epsilon}))$ , and $\mathcal{O}(\log({1}/{\epsilon}))$ in terms of the target accuracy $\epsilon$ , respectively. While, in terms of the condition number $\kappa$ , it has a higher complexity compared with the algorithms listed in Table I.

Notice that the convergence result and the complexity result of Theorem 1 are established under the assumption that $T\geq 1$ and $J\geq 1$ . In the following, we show that by imposing some restrictions on $T$ and $J$ , SSGD can obtain a lower computational complexity.

(a)

(b)

(c)

Figure 2: Comparison of the simple SGD type algorithms for solving problem (15) whose data sets are constructed by setting

\boldsymbol{w}_{0}=(2,5,7)

(a)

(b)

(c)

Figure 3: Comparison of the simple SGD type algorithms for solving problem (15) whose data sets are constructed by setting

\boldsymbol{w}_{0}=(1,4,6,\ldots,6)\in\mathbb{R}^{500}

Theorem 2.

Apply SSGD to solve problem (2). Suppose Assumptions 1, 2, and 3 hold. Define $r_{v}=1$ , $r_{w}=1$ , $\rho_{2}=\mathcal{O}(\kappa^{-3})$ , $\rho_{1}=2\rho_{2}(1+r_{w})C_{23}$ , and

L_{31}=\frac{L_{\Phi}}{2}+2\left(\frac{\rho_{1}}{4\mu^{2}}L^{2}+\rho_{2}(1+r_{w})C_{22}\right),

where $\kappa={L}/{\mu}$ , $L_{\Phi}=\mathcal{O}(\kappa^{3})$ is defined in Proposition 1, $C_{22}=\left(1+{1}/{r_{w}}\right)\left(2+8{L^{2}}/{\mu^{2}}\right)C_{1}^{2}$ , $C_{23}=\left(1+{1}/{r_{w}}\right)(5+r_{w})C_{1}^{2}$ , and $C_{1}={M}/{\mu^{2}}\rho+{L}/{\mu}$ . Let $\|\boldsymbol{v}^{0}\|\leq M/\mu$ , and select $\eta={1}/{(2L)}$ , $\beta={1}/{(2L)}$ ,

\alpha=\min\left\{\frac{1}{2L_{31}},~{}\frac{1}{4L^{2}}\rho_{2},~{}\frac{1}{2}\rho_{2}(1+r_{w})\frac{C_{23}}{\bar{C}_{1}^{2}}\right\},

T\geq\frac{\log(\rho_{1}/(8(\rho_{1}+C_{24}\rho_{2})))}{\log(L/(\mu+L))}=\mathcal{O}(\kappa),

J\geq\frac{\log(1/(4(1+r_{w})^{2}))}{\log(1-\eta\mu)}=\mathcal{O}(\kappa),

where $\bar{C}_{1}=L+{M\tau}/{\mu}=\mathcal{O}({\kappa})$ , and $C_{24}=\left(1+{1}/{r_{w}}\right)(9+8r_{w})C_{1}^{2}$ . Then, we have

	$\displaystyle\frac{1}{K}\sum_{k=0}^{K-1}\mathbb{E}\\|\nabla\Phi(\boldsymbol{x}^{k})\\|^{2}$
	$\displaystyle\leq\mathcal{O}\left(\frac{\kappa^{3}}{K}+\kappa\frac{\sigma_{f}^{2}}{D_{f}}+\kappa^{2}\frac{\sigma_{g,1}^{2}}{D_{g}}+\kappa^{5}\frac{1}{S}\sigma_{g}^{2}+\kappa^{3}\frac{\sigma_{g,2}^{2}}{D}\right).$

Furthermore, a computational complexity of the order of $\mathcal{O}(\kappa^{9}\epsilon^{-2})$ is sufficient for SSGD to reach an $\epsilon$ -accurate stationary point.

In Theorem 2, the convergence of SSGD is established under the assumption that $T\geq\mathcal{O}(\kappa)$ and $J\geq\mathcal{O}(\kappa)$ . Furthermore, a computational complexity by an order of $\mathcal{O}(\kappa^{9}\epsilon^{-2})$ is obtained in Theorem 2, which is lower than the lowest computational complexity $\mathcal{O}(\kappa^{13}\epsilon^{-2})$ obtained in Theorem 1. Compared with the algorithms listed in Table I, SSGD achieves the lowest computational complexity.

V Experimental Results

In this section, we first validate our theoretical results on synthetic bilevel optimization problems. Then, we apply our proposed algorithm on a hyperparameter optimization problem.

V-A Synthetic Bilevel Optimization Problems

In the following, we perform experiments on the bilevel problem in (15). The experimental details are in the supplementary material.

Firstly, we compare our proposed algorithm SSGD with the existing simple SGD type algorithms BSA, TTSA, and stocBiO. In the experiments, we consider two datasets, i.e., datasets constructed by setting $\boldsymbol{w}_{0}=(2,5,7)$ and datasets constructed by setting $\boldsymbol{w}_{0}=(1,4,6,\ldots,6)\in\mathbb{R}^{500}$ .

From Fig. 2LABEL:sub@fig_third_case and Fig. 3LABEL:sub@fig_third_case, we observe that SSGD and stocBiO run almost the same number of iteration steps $K$ to converge to the stationary point with $\|\nabla\Phi(\boldsymbol{x})\|=0$ , less than that of BSA and TTSA. However, compared with SSGD, a larger $J$ to estimate the hypergradient is required for stocBiO to converge to the stationary point with $\|\nabla\Phi(\boldsymbol{x})\|=0$ . As a result, Fig. 2LABEL:sub@fig_second_case and Fig. 3LABEL:sub@fig_second_case show that SSGD converges fastest among all the compared methods. Furthermore, from Fig. 2LABEL:sub@fig_first_case and Fig. 3LABEL:sub@fig_first_case, we observe that SSGD converges to the point with the smallest upper objective function value for both cases.

Then, we study the influence of iteration step $T$ , iteration step $J$ , batch size, and stepsize $\alpha$ on the convergence behavior of SSGD in Algorithm 2. In the experiments, the datasets are constructed by setting $\boldsymbol{w}_{0}=(2,5,7,\ldots,7)\in\mathbb{R}^{1000}$ .

Fig. 4LABEL:sub@fig_first_case shows the convergence behavior of SSGD under different choices of $T\in\{1,5,30,50\}$ , where SSGD_ $T$ _ $i$ indicates that $T$ is set to be $i$ . It shows that beginning from 1, increasing $T$ can speed up the convergence of SSGD, but when $T$ reaches $5$ , increasing $T$ can slow down the convergence speed of SSGD. Fig. 4LABEL:sub@fig_second_case shows the convergence behavior of SSGD under different choices of $J\in\{1,5,20,30\}$ , where SSGD_ $J$ _ $i$ indicates that $J$ is set to be $i$ . We observe that SSGD converges fastest when $J$ is chosen to be 5.

Fig. 5LABEL:sub@fig_first_case shows the convergence behavior of SSGD under different batch sizes with SSGD_ $i$ indicating that the batch size is set to be $i$ . It shows that SSGD converges fastest when batch size is set to be 1, and increasing batch size can slow down the convergence speed of SSGD. Fig. 5LABEL:sub@fig_second_case shows the convergence behavior of SSGD under different stepsizes $\alpha$ . It shows that choosing a larger stepsize $\alpha$ allows SSGD to converge faster.

Fig. 6 shows the convergence behavior of SSGD under different choices of $T$ and $J$ , where SSGD_ $T$ _ $t$ _ $J$ _ $j$ indicates that $T$ is set to be $t$ and $J$ is set to be $j$ . It shows that SSGD converges faster in the double loop case(i.e., $T>1$ ) than in the single loop case(i.e., $T=1$ ).

V-B Hyper-parameter Optimization

In this section, we use a hyperparameter optimization problem, i.e., data hyper-cleaning[5], to evaluate the performance of our algorithm. Data hyper-cleaning is to train a classifier under the case that some data in the training set is corrupted. Following [5], when we train a linear classifier, data hyper-cleaning task can be expressed as a bilvel optimization problem whose UL function is given by

f(\boldsymbol{x},\boldsymbol{y})=\frac{1}{|\mathcal{D}_{\text{val}}|}\sum\limits_{(\boldsymbol{u}_{i},v_{i})\in\mathcal{D}_{\text{val}}}L({\boldsymbol{y}}^{\top}\boldsymbol{u}_{i},v_{i})

and whose LL function is given by

g(\boldsymbol{x},\boldsymbol{y})=\frac{1}{|\mathcal{D}_{\text{tr}}|}\sum\limits_{(\boldsymbol{u}_{i},v_{i})\in\mathcal{D}_{\text{tr}}}\sigma(x_{i})L(\boldsymbol{y}^{\top}\boldsymbol{u}_{i},v_{i})+c\|\boldsymbol{y}\|^{2}

where $L$ denotes the cross-entropy loss, $\sigma(\cdot)$ is the sigmoid function, $c$ is the coefficient of the regular term and is set to be $0.001$ , and $\mathcal{D}_{\text{val}}$ and $\mathcal{D}_{\text{tr}}$ denote the validation set and training set, respectively.

TABLE II: Comparison of the simple SGD type algorithms for data hyper-cleaning tasks on MNIST, FashionMNIST and CIFAR10 datasets. Acc.(%) and F1 score(%) denote the test accuracy and the harmonic mean of the precision and recall, respectively

Method	MNIST		FashionMNIST		CIFAR10
Method	Acc.	F1 score	Acc.	F1 score	Acc.	F1 score
stocBiO	89.90	63.92	82.18	65.17	38.36	64.07
BSA	80.67	14.20	72.57	16.51	22.72	21.09
TTSA	88.88	33.94	80.62	32.54	35.69	40.03
SSGD	90.32	83.81	82.32	80.45	38.35	66.16

In the following experiments, we compare our proposed algorithm SSGD with BSA, TTSA, stocBiO on three datasets: MNIST [26], FashionMNIST [27] and CIFAR10. Each sample in the training set is corrupted with probability $0.3$ . More results and details are in the supplementary material.

Table II shows the experimental results on three datasets over 270s. We observe that SSGD achieves the highest test accuracy on MNIST and FashionMNIST datasets, and achieves the highest F1 score on all the datasets. It shows that SSGD can obtain competetive results on data hyper-cleaning problem.

VI Conclusion

In the paper, we propose to use SGD-based Estimation to estimate the hypergradient for the nonconvex-strongly-convex bilevel optimization problems, and propose a novel simple SGD type algorithm. We provide the convergence guarantee for our algorithm, and prove that our proposed algorithm can achieve lower computational complexity compared with all the existing simple SGD type algorithms for bilevel optimization. Moreover, the experiments validate our theoretical results.

Acknowledgments

The work is partially supported by National Key R&D Program of China(2020YFB1313503 and 2018AAA0100300), National Natural Science Foundation of China(Nos. 61922019, 61733002, 61672125, and 61976041), and LiaoNing Revitalization Talents Program(XLYC1807088).

References

[1] L. Franceschi, P. Frasconi, S. Salzo, R. Grazzi, and M. Pontil, “Bilevel programming for hyperparameter optimization and meta-learning,” in Proc. 35th Int. Conf. Mach. Learn. (ICML), Jul. 2018, pp. 1568–1577.
[2] L. Bertinetto, J. F. Henriques, P. Torr, and A. Vedaldi, “Meta-learning with differentiable closed-form solvers,” in Proc. Int. Conf. Learn. Represent. (ICLR), 2019.
[3] A. Rajeswaran, C. Finn, S. Kakade, and S. Levine, “Meta-learning with implicit gradients,” in Advances in Neural Information Processing Systems, 2019.
[4] K. Ji, J. D. Lee, Y. Liang, and H. V. Poor, “Convergence of meta-learning with task-specific adaptation over partial parameters,” in Proc. 34th Conf. Neural Inf. Process. Syst. (NeurIPS), 2020.
[5] A. Shaban, C. A. Cheng, N. Hatch, and B. Boots, “Truncated back-propagation for bilevel optimization,” in Proc. 22th Int. Conf. Artif. Intell. Statist., vol. 89, 2019, pp. 1723–1732.
[6] V. Konda and J. Tsitsiklis, “Actor-critic algorithms,” in Advances in Neural Information Processing Systems, vol. 12, 1999.
[7] M. Hong, H. T. Wai, Z. Wang, and Z. Yang, “A two-timescale framework for bilevel optimization: Complexity analysis and application to actor-critic,” 2020, arXiv: 2007.05170. [Online]. Available: https://arxiv.org/abs/2007.05170
[8] P. Hansen, B. Jaumard, and G. Savard, “New branch-and-bound rules for linear bilevel programming,” SIAM J. Sci. Stat. Comput., vol. 13, no. 5, pp. 1194–1217, 1992.
[9] C. Shi, J. Lu, and G. Zhang, “An extended kuhn–tucker approach for linear bilevel programming,” Appl. Math. Comput., vol. 162, no. 1, pp. 51–63, 2005.
[10] Y. Lv, T. Hu, G. Wang, and Z. Wan, “A penalty function method based on kuhn–tucker condition for solving linear bilevel programming,” Appl. Math. Comput., vol. 188, no. 1, pp. 808–813, 2007.
[11] S. Qin, X. Le, and J. Wang, “A neurodynamic optimization approach to bilevel quadratic programming,” IEEE Trans. Neural Netw. Learn. Syst., vol. 28, no. 11, pp. 2580–2591, 2017.
[12] J. Domke, “Generic methods for optimization-based modeling,” in Proc. 15th Int. Conf. Artif. Intell. Statist., Apr. 2012, pp. 318–326.
[13] L. Franceschi, M. Donini, P. Frasconi, and M. Pontil, “Forward and reverse gradient-based hyperparameter optimization,” in Proc. 34th Int. Conf. Mach. Learn. (ICML), vol. 70, Aug. 2017, pp. 1165–1173.
[14] Y. L. K. Ji, J. Yang, “Bilevel optimization: Convergence analysis and enhanced design,” in Proc. Int. Conf. Mach. Learn. (ICML), 2021.
[15] F. Pedregosa, “Hyperparameter optimization with approximate gradient,” in Proc. Int. Conf. Mach. Learn. (ICML), 2016, p. 737–746.
[16] J. Lorraine, P. Vicol, and D. Duvenaud, “Optimizing millions of hyperparameters by implicit differentiation,” in Proc. 23th Int. Conf. Artif. Intell. Statist., 2020, pp. 1540–1552.
[17] S. Ghadimi and M. Wang, “Approximation methods for bilevel programming,” 2018, arXiv: 1802.02246. [Online]. Available: https://arxiv.org/abs/1802.02246
[18] D. Maclaurin, D. Duvenaud, and R. P. Adams, “Gradient-based hyperparameter optimization through reversible learning,” in Proc. Int. Conf. Mach. Learn. (ICML), 2015.
[19] R. Liu, P. Mu, X. Yuan, S. Zeng, and J. Zhang, “A generic first-order algorithmic framework for bi-level programming beyond lower-level singleton,” in Proc. Int. Conf. Mach. Learn. (ICML), 2020.
[20] P. Khanduri, S. Zeng, M. Hong, H. T. Wai, Z. Wang, and Z. Yang, “A momentum-assisted single-timescale stochastic approximation algorithm for bilevel optimization,” 2021, arXiv: 2102.07367v1. [Online]. Available: https://arxiv.org/abs/2102.07367v1
[21] Z. Guo, Y. Xu, W. Yin, R. Jin, and T. Yang, “On stochastic moving-average estimators for non-convex optimization,” 2021, arXiv: 2104.14840. [Online]. Available: https://arxiv.org/abs/2104.14840
[22] J. Yang, K. Ji, and Y. Liang, “Provably faster algorithms for bilevel optimization,” in Advances in Neural Information Processing Systems, vol. 34, 2021, pp. 13 670–13 682.
[23] H. H. J. Li, B. Gu, “A fully single loop algorithm for bilevel optimization without hessian inverse,” in Proc. AAAI Conf. Artif. Intell., vol. 36, 2022, pp. 7426–7434.
[24] A. Cutkosky and F. Orabona, “Momentum-based variance reduction in non-convex sgd,” in Advances in Neural Information Processing Systems, vol. 32, 2019.
[25] T. Chen, Y. Sun, and W. Yin, “Closing the gap: Tighter analysis of alternating stochastic gradient methods for bilevel problems,” in Advances in Neural Information Processing Systems, vol. 34, 2021, pp. 25 294–25 307.
[26] Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner, “Gradient-based learning applied to document recognition,” Proc. IEEE, vol. 86, pp. 2278–2324, 1998.
[27] H. Xiao, K. Rasul, and R. Vollgraf, “Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms,” 2017, arXiv: 1708.07747. [Online]. Available: https://arxiv.org/abs/1708.07747

Haimei Huo received the B.S. degree in information and computing science from Changchun University of Science and Technology, China, in 2017. She is currently pursuing the PhD degree in computational mathematics from Dalian University of Technology, Dalian, China. Her research interests include machine learning and optimization.

Risheng Liu received his B.Sc. (2007) and Ph.D. (2012) from Dalian University of Technology, China. From 2010 to 2012, he was doing research as joint Ph.D. in robotics institute at Carnegie Mellon University. From 2016 to 2018, he was doing research as Hong Kong Scholar at the Hong Kong Polytechnic University. He is currently a full professor with the Digital Media Department at International School of Information Science & Engineering, Dalian University of Technology. He was awarded the ”Outstanding Youth Science Foundation” of the National Natural Science Foundation of China. His research interests include optimization, computer vision, and multimedia.

Zhixun Su received the B.S. degree in mathematics from Jilin University, Changchun, China in 1987, the M.S. degree in computer science from Nankai University, Tianjin, China in 1990, and the Ph.D. degree in computational mathematics from Dalian University of Technology, Dalian, China in 1993. He has been a Professor in the School of Mathematical Sciences, Dalian University of Technology, since 1999. He is currently the director of the Key Laboratory for Computational Mathematics and Data Intelligence of Liaoning Province and the vise director of Liaoning Center for Applied Mathematics. His current research interests include computer graphics, computer vision, computational geometry, and machine learning.

	$\displaystyle\mathbb{E}\\|\nabla F(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\xi})-\nabla f(\boldsymbol{x},\boldsymbol{y})\\|^{2}\leq\sigma_{f}^{2},$
	$\displaystyle\mathbb{E}\\|\nabla G(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\zeta})-\nabla g(\boldsymbol{x},\boldsymbol{y})\\|^{2}\leq\sigma_{g}^{2},$
	$\displaystyle\mathbb{E}\\|\nabla_{\boldsymbol{x}}\nabla_{\boldsymbol{y}}G(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\zeta})-\nabla_{\boldsymbol{x}}\nabla_{\boldsymbol{y}}g(\boldsymbol{x},\boldsymbol{y})\\|^{2}\leq\sigma_{g,1}^{2},$
	$\displaystyle\mathbb{E}\\|\nabla_{\boldsymbol{y}}^{2}G(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\zeta})-\nabla_{\boldsymbol{y}}^{2}g(\boldsymbol{x},\boldsymbol{y})\\|^{2}\leq\sigma_{g,2}^{2},$

	$\displaystyle\\|\nabla\Phi(\boldsymbol{x}^{k})-\bar{\boldsymbol{h}}_{f}^{k}\\|$
	$\displaystyle\leq\bigg{(}L+\frac{M}{\mu}\tau\bigg{)}\\|\boldsymbol{y}^{k,T}-\boldsymbol{y}^{*}(\boldsymbol{x}^{k})\\|$
	$\displaystyle+L(1-\eta\mu)^{J}\\|\boldsymbol{v}^{k,0}-\boldsymbol{v}^{*}\\|$
	$\displaystyle+L\eta\bigg{(}\frac{M}{\mu}\rho+L\bigg{)}\sum\limits_{t=0}^{J-1}(1-\eta\mu)^{t}\\|\boldsymbol{y}^{k,T}-\boldsymbol{y}^{*}(\boldsymbol{x}^{k})\\|;$