---
abstract: 'We study a class of bilevel convex optimization problems where the goal is to find the minimizer of [[an objective function]{}]{} in the upper level, among the set of all optimal solutions of [[an optimization]{}]{} problem in the lower level. A wide range of problems in convex optimization can be formulated [[using]{}]{} this class. [[An important example is the case where an optimization problem is ill-posed.]{}]{} [[In this paper, our interest lies in addressing the bilevel problems, where the lower level objective is given as a finite sum of separate nondifferentiable convex component functions. This is the case]{}]{} in a variety of applications in distributed optimization, such as large-scale data processing in machine learning and neural networks. [[To the best of our knowledge, this class of bilevel problems, with a finite sum in the lower level, has not been addressed before. Motivated by this gap, we]{}]{} develop an iterative regularized incremental subgradient method[[,]{}]{} [[where the agents update their iterates in a cyclic manner using a regularized subgradient]{}]{}. Under a suitable choice of [[the]{}]{} regularization [[parameter sequence]{}]{}, we establish the convergence of the proposed algorithm and [[derive a]{}]{} rate of $\mathcal{\cal O} \left({1}/k^{0.5-\epsilon}\right)$ [[in terms of]{}]{} the lower level objective function for an arbitrary small $\epsilon>0$. [[We]{}]{} present the performance of the algorithm on a binary text classification problem.'
author:
- 'Mostafa Amini$^{1}$ and Farzad Yousefian$^{2}$[^1]'
title: '**An Iterative Regularized Incremental Projected Subgradient Method for a Class of Bilevel Optimization Problems** '
---

Introduction
============

In this paper, we consider a class of bilevel optimization problems as follows $$\begin{aligned}
\label{def:SL-INC} \tag{$P_f^h$}
\displaystyle \mbox{minimize}& \qquad h(x)\\
\mbox{subject to} & \qquad x \in X^* \triangleq \arg\min_{y\in X}f(y) \notag.\end{aligned}$$ where $f,h :\mathbb{R}^n \rightarrow \mathbb{R}$ [[denote the lower and upper level objective functions, respectively, and $X \subseteq \mathbb{R}^n $ is a constraint set]{}]{}. This is called the *selection problem* ([@Tseng07; @Sabach17]) as we are selecting among optimal solutions of a lower level problem, one that minimizes the objective function $h$. [[In particular]{}]{}, we [[consider the case where the lower level objective function is]{}]{} given as $f(x) \triangleq \sum_{i=1}^{m} f_i(x)$, where $f_i:\mathbb{R}^n \rightarrow \mathbb{R}$ is the $i$th component function for $i=1,\cdots,m$. We make the following basic assumptions.

\[assum:properties\] 

The set $X \subset \mathbb{R}^n$ is nonempty, compact and convex; also $X\subseteq int(dom(f) \cap dom(h))$.

The functions $f_i(x)$ for $i=1,\cdots, m$ are proper, closed, convex, and possibly nondifferentiable.

The function $h$ is strongly convex with parameter $\mu_h>0$ and possibly nondifferentiable.

[[Next]{}]{}, we [[present]{}]{} two instances of the applications of formulation .

Example problems
----------------

[**(i) Constrained nonlinear optimization:**]{} [[Consider a constrained]{}]{} convex optimization problem [[given]{}]{} as $$\begin{aligned}
\displaystyle \mbox{minimize}& \qquad h(x)\\
\mbox{subject to} &\qquad q_i(x) \leq 0, \ \hbox{for } i=1,\cdots,m\\
&\qquad x\in X\end{aligned}$$ where $X \subseteq \mathbb{R}^n$ [[is an easy-to-project constraint set]{}]{}, $h,q_i :\mathbb{R}^n \rightarrow \mathbb{R} $ for all $i=1,\cdots,m$ [[are convex (and possibly nonlinear) funcitons]{}]{}. This problem can be reformulated as by setting (cf. [@Solodov072]) $$\begin{aligned}
f(x) \triangleq \sum_{i=1}^{m} f_i(x)= \sum_{i=1}^{m} \max \{0, q_i(x)\}.\end{aligned}$$ [**(ii) Ill-posed distributed optimization:**]{} An [[optimization]{}]{} problem is called ill-posed when it has multiple optimal solutions or it is very sensitive to data perturbations [@Tikhonov77]. For instance, in applications arising in machine learning, consider the *empirical risk minimization problem* where the goal is to minimize the total loss $\sum_{i=1}^{m} \mathfrak{L} (a_ix,b_i)$, where $a_i$ is the input, $b_i$ is the output of $i$th observed datum and $\mathfrak{L}$ is the loss function. For example, in *logisitic loss regression*, $\mathfrak{L}$ is merely convex. In these cases, another criterion such as sparsity may be taken into account for the optimal solution. So, to induce sparsity, a secondary objective function $h$ is considered in the given problem. For instance, the well-known *elastic net* regularization can be used as function $h$. Hence, to address ill-posedness, the following bilevel optimization model is considered [@Tseng07; @Sabach17]: $$\begin{aligned}
\displaystyle \mbox{minimize}& \qquad \|x\|_1 + \mu \|x\|_2^2\\
\mbox{subject to} &\qquad x \in \arg \min_{y \in X} \sum_{i=1}^{m} \mathfrak{L} (a_i^T y,b_i),\end{aligned}$$\
where $\mu>0$ regulates the trade-off between $\ell_1$ and $\ell_2$ norms.

\[fig:fiveplots\]

Existing methods
----------------

Problem [[,]{}]{} that is also referred to as *hierarchical optimization*[[,]{}]{} is a particular case of *mathematical program with generalized equation (or equilibrium) constraint* [@Kocvara04; @Pang96]. There has been a few approaches to tackle this problem. Note that in all approaches the following minimization problem and its minimizer have been extensively utilized.

\[def:reg minimizer\] [[Given a parameter $\lambda > 0$]{}]{}, the regularized problem corresponding to , [[is defined]{}]{} as $$\begin{aligned}
\label{reg form INC}
    \displaystyle \mbox{minimize}& \qquad f_\lambda (x) \triangleq f(x)+ \lambda h(x)\\
    \mbox{subject to} &\qquad x \in X\notag.
    \end{aligned}$$ Also, let $x_\lambda^*$ denote the unique minimizer of this problem. We may categorize the existing algorithms as follows.

[**(i) Exact regularization:**]{} The regularization technique has been highly used in some applications such as signal processing with $h(x)= \|x\|_2^2$ or $h(x)= \|x\|_1$ [@Tikhonov77; @Beck09]. This technique needs a proper parameter $\lambda$ which is difficult to determine in most of cases. To address this issue, Mangasarian et al. [@Mangasarian79; @Mangasarian85] introduced *exact regularization*. A solution of problem is called *exact* when it is in the set $X^*$. The main drawback of this approach is that the threshold below which for any $\lambda$ the regularization is exact, is very difficult to determine a priori [[(see [@Tseng07])]{}]{}.

[**(ii) Iterative regularization:**]{} In this approach[[,]{}]{} the idea is to develop a single-loop scheme where the regularization parameter is updated iteratively during the algorithm. In [@Solodov07], an *explicit* descent algorithm is proposed, where problem is solved as [[a single]{}]{}-level unconstrained problem iteratively. [[In the smooth case,]{}]{} the convergence is shown when $\sum_{k=1}^{\infty} \lambda_k =\infty$ and $\lim_{k\rightarrow \infty} \lambda_k=0$. For nonsmooth cases, a bundle method was proposed, which has a descent step for the weighted combination of objective functions in the lower and upper levels [@Solodov072]. Another algorithm called *hybrid steepest descent method* (HSDM) and its extensions were developed in [@Yamada11; @Neto11]. The main drawback of all these works is that no complexity analysis was provided for the [[underlying algrotihm]{}]{}. In [@Yousefian17], the rate of $\mathcal{\cal O} (1/k^{1/6-\epsilon})$ is derived for a special case of problem [[,]{}]{} where the objective function is $\|\cdot\|_2$ and the constraint is solutions of a stochastic variational inequality problem. Recently, a primal-dual algorithm in [@Rosasco18] was offered that uses the idea of iterative regularization. Despite [[a sublinear rate in terms of $h$]{}]{}, the algorithm can only be applied to continuously differentiable [[and small-scale (i.e., $m=1$)]{}]{} regimes. [**(iii) Minimal norm gradient:**]{} In [@Beck14], the *minimal norm gradient* (MNG) method was developed for solving [[problem with $m=1$]{}]{}. The rate of $\mathcal{\cal O} (1/\sqrt{k})$ was derived for the convergence with respect to lower level problem. The main disadvantage of the MNG method is that it is a two-loop scheme where at each iteration a minimization problem should be solved.

[**(iv) Sequential averaging:**]{} The *sequential averaging method* (SAM), developed in [@Xu04], was [[employed]{}]{} in [@Sabach17] for solving the problem in a more general setting. The proposed method is proved to have the rate of convergence of $\mathcal{\cal O} (1/k)$ in terms of the function $f$. The method is called *Big-SAM*. [[Despite that it is a single-loop scheme, sequential averaging schemes require smoothness properties of the problem and seem not to lend themselves to distributed implementations.]{}]{}

Main Contributions
------------------

[[For more details on the main distinctions between our work and the existing methods see Table \[fig:fiveplots\]. In none of the existing methods,]{}]{} the finite sum form for the lower level problem is [[considered]{}]{}. The sum structure is very rampant in practice when we have separate objective functions related to different agents in a distributed setting. This is the case for example in machine learning for very large datasets [@Bottou05], where each $f_i$ represents an agent that is cooperating with others. When the complete information of all the agents, i.e. summation of all (sub)gradients is not available, these agents can be treated distinctly. Due to its wide range of applications in distributed optimization, finite sum problem has been extensively studied. Among popular methods are *incremental (sub)gradient* (IG) [@Bertsekas99; @Nedich01] and *incremental aggregated (sub)gradients* (IAG) [@Blatt07; @Tseng14] for deterministic and *stochastic average gradient* (SAG) [@Roux12], SAGA [@Defazio14] and MISO methods [@Mairal13] for stochastic regimes. These algorithms have faster convergence and are computationally efficient in large-scale optimization since a very less amount of memory is required at each step in order to store only one agent’s information and subsequently update the iterate based on that [@Ozdaglar17]. Despite the widespread use of these first-order methods, they do not address [[the]{}]{} bilevel problem .

Motivated by the existing lack in the literature and inspired by the advantages of incremental approaches, in this paper, we let the lower level objective function to be a summation of $m$ components. Then we use the idea of incremental subgradient optimization to address problem . We let functions in both levels to be nondifferentiable. We then prove the convergence of our proposed algorithm as well as the $\mathcal{\cal O} (1/k^{0.5-\epsilon})$ rate of convergence.

An interesting research question that is remained as a future direction to our research is if we can establish the convergence of iterative regularized IAG method in solving problem or similarly SAG, SAGA and MISO in stochastic regimes.

[**Notation**]{} The inner product of two vectors $x,y \in \mathbb{R}^n$, is shown as $x^Ty$. Also, $\|\cdot\|$ denotes Euclidean norm known as $\|\cdot\|_2$. For a convex function $f$ with the domain dom$(f)$, any vector $g_f$ with $f(x)+g_f^T(y-x) \leq f(y)$ for all $x,y \in \hbox{dom}(f)$, is called a subgradient of $f$ at $x$. We let $\partial f(x)$ and $\partial  h(x)$ denote the set of all subgradients of functions $f$ and $h$ at $x$. Let $f^*$ be the optimal value and $X^*$ represent the set of all optimal solutions of [[the]{}]{} lower level problem in and $x^*$ shows any element of this set . Likewise, $x^*_h$ [[denotes]{}]{} the optimal solution of problem [[, and]{}]{} $x^*_{\lambda}$ denotes the optimal solution of problem . [[Also, we let $\mathcal{P}_X(x)$ denote the Euclidean projection of vector $x$ onto the set $X$.]{}]{}

The rest of this paper is organized as follows. In Section \[sec: alg out\], we present the algorithm outline. Then, we discuss the convergence analysis in Section \[sec: conv ana\], and derive the convergence rate in Section \[sec:rate result\]. We present the numerical results in Section \[sec:num\], and conclude in Section \[sec:rem\].

Algorithm outline {#sec: alg out}
=================

In this section, we introduce the *iterative regularized incremental projected (sub)gradient* (IR-IG) for generating a sequence that converges to the unique optimal solution of . See Algorithm \[algorithm:IRIG\].

. $$\begin{aligned}
\label{mainstep}
        x_{k,i+1} := \mathcal{P}_X \left(x_{k,i} - \gamma_k\left(g_{f_{i+1}}(x_{k,i}) + \frac{\lambda_k}{m} g_h(x_{k,i})\right) \right).
        \end{aligned}$$ $$\begin{aligned}
        S_{k+1}:=S_k+\gamma_{k+1}^r, \quad          \bar x_{k+1}:=\frac{S_k \bar x_k+\gamma_{k+1}^r x_{k+1}}{S_{k+1}}.\label{def:averagingII} 
        \end{aligned}$$ $\bar x_{N}.$

IR-IG method includes two main steps. First, the agents update their iterates in an incremental fashion similar to the standard IG method. This step takes a circle around the all components of function $f$ to update the iterate. However, The main difference lies in the secondary objective function $h$, which is added by a vanishing multiplier $\lambda_k$. Second, we do averaging in order to accelerate the convergence speed of the algorithm. For this, we consider a weighted average sequence $\{\bar x_k\}$ defined as below: $$\begin{aligned}
\label{weighted alg INC}
    \bar x_{k+1} := \sum_{t=0}^{k} \psi_{t,k} x_t,  \hbox{ where } \psi_{t,k} \triangleq \frac{\gamma_t^r}{\sum_{i=0}^{k} \gamma_i^r},  
    \end{aligned}$$ in which $r<1$ is a constant, controlling the weights. Note that in Algorithm \[algorithm:IRIG\] follows from the relation by applying induction, (see e.g., Proposition 3 in [@Yousefian18]).

Convergence analysis {#sec: conv ana}
====================

In this section, [[our goal is to show that]{}]{} the generated sequence $\{\bar x_k\}$ by Algorithm \[algorithm:IRIG\] converges to the unique optimal solution of problem [[(see Theorem \[thm conv for xbar INC\])]{}]{}.

\[bound on gradient by beck\] (a) Note that from Theorem 3.16, pg. 42 of [@Beck17], Assumption \[assum:properties\] implies that there exist constants $C_f, C_h \in \mathbb{R}$ such that $\|g_{f_i}(x)\| \leq C_f$ and $\|g_{h}(x)\| \leq C_h$ for all $i=1,\cdots, m$ and $x\in X$, where $g_{f_i}(x) \in \partial f_i(x)$ and $g_{h}(x) \in \partial h(x)$.\
(b) From Theorem 3.61, pg. 71 of [@Beck17], functions $f_i$ and $h$ are Lipschitz over $X$ with parameters $C_f$ and $C_h$, respectively, i.e., for all $i=1,\cdots, m$ and $x,y\in X$ $$\begin{aligned}
        |f_i(x)-f_i(y)| \leq C_f \|x-y\|, \quad   |h(x)-h(y)| \leq C_h \|x-y\|.
        \end{aligned}$$ [[(c) Assumption \[assum:properties\](a,b) imply that the optimal solution set, $X^*$, is nonempty.]{}]{}

Here, we start with a lemma which helps bound the error of optimal solutions of the problem for two different values of $\lambda$. We will make use of this lemma in the convergence analysis. The proof for this lemma can be done in a same fashion to that of Proposition 1 in [@Yousefian17].

\[lemma:err\_inequalitygeneral\] Let Assumption \[assum:properties\] hold. Suppose $\{x_{\lambda_k}^* \}$ be the sequence of the optimal solutions of problem with parameter $\lambda:= \lambda_k$. Then,

-   $\|x_{\lambda_k}^*-x_{\lambda_{k-1}}^*\|\leq \frac{C_h}{\mu_h}\left  |1-\frac{\lambda_{k-1}}{\lambda_k} \right|$.

-   If $\lambda_k \rightarrow 0$, then the sequence $\{x_{\lambda_k}^*\}$ converges to the unique optimal solution of problem , i.e., $x_{h}^*$.

To get started, we also need a recursive upper bound on the term $ \|x_{k+1} - x_{\lambda_k}^*  \|$. This is provided by the following lemma and will be used in Proposition \[conv x\_k INC\] to prove the convergence of sequence $\{x_k\}$ generated by the algorithm to $x_h^*$.

\[[**A recursive error bound**]{}\] \[lemma: an error bound INC\] Let Assumption \[assum:properties\] hold and $0<\mu_k\lambda_k \mu_h \leq 2m$. Then, for the sequence $\{x_k\}$ generated by Algorithm \[algorithm:IRIG\] and for all $k>0$ we have $$\begin{aligned}
    \left \|x_{k+1} - x_{\lambda_k}^* \right \|^2 & \leq \left( 1- \frac{\gamma_k \lambda_k \mu_h}{2m} \right) \left \|x_{k} - x_{\lambda_{k-1}}^* \right \|^2 \\ & + \frac{3m C_h^2}{\gamma_k \lambda_k \mu_h^3}\left  |1-\frac{\lambda_{k-1}}{\lambda_k} \right|^2 +  6 m^2 \gamma_k^2 (C_f^2 + \lambda_k^2 C_h^2),
    \end{aligned}$$ where $x_{\lambda_k}^*$ is the unique optimal solution of problem with $\lambda:=\lambda_k$.

Using and the nonexpansiveness property of projection, we have\
$    \left \|x_{k,i+1} - x_{\lambda_k}^* \right \|^2  \\= \left \|\mathcal{P}_X \left(x_{k,i} - \gamma_k\left(g_{f_{i+1}}(x_{k,i}) +  \frac{\lambda_k}{m} g_h(x_{k,i})\right) \right) - \mathcal{P}_X(x_{\lambda_k}^*) \right \|^2 \\
         \leq \left \|x_{k,i} - x_{\lambda_k}^* \right \|^2 + \gamma_k^2 \left \|g_{f_{i+1}}(x_{k,i}) + \frac{\lambda_k}{m} g_h(x_{k,i}) \right \|^2 \\-2\gamma_k \left ( g_{f_{i+1}}(x_{k,i}) +  \frac{\lambda_k}{m} g_h(x_{k,i})\right )^T\left( x_{k,i} - x_{\lambda_k}^*\right ).\\
    $

By boundedness of subgradients from Remark \[bound on gradient by beck\](a), the definition of subgradient for $f_{i+1}$, and the strong convexity of $h$, we obtain $$\begin{aligned}
    & \left \|x_{k,i+1} - x_{\lambda_k}^* \right \|^2\\ & \leq \left \|x_{k,i} - x_{\lambda_k}^* \right \|^2 +2 \gamma_k^2 (C_f^2 + \lambda_k^2 C_h^2)  - 2\gamma_k \left( f_{i+1}(x_{k,i})- f_{i+1}(x_{\lambda_k}^*) \right) \\
    & -  \frac{2\gamma_k \lambda_k}{m} \left( h(x_{k,i}) - h(x_{\lambda_k}^*) \right) - \frac{\gamma_k \lambda_k \mu_h}{m} \left \|x_{k,i} - x_{\lambda_k}^* \right \|^2 \\
    & = \left( 1- \frac{\gamma_k \lambda_k \mu_h}{m} \right) \left \| x_{k,i} - x_{\lambda_k}^* \right \|^2 - 2\gamma_k \left( f_{i+1}(x_{k,i}) + \frac{\lambda_k}{m} h(x_{k,i}) \right) \\& + 2\gamma_k \left(f_{i+1}(x_{\lambda_k}^*)+ \frac{\lambda_k}{m} h(x_{\lambda_k}^*) \right) +2 \gamma_k^2 (C_f^2 + \lambda_k^2 C_h^2).
    \end{aligned}$$ Taking summation from both sides over $i$, using $x_{k,0}=x_k,  x_{k,m}=x_{k+1}$, and that $\gamma_k \lambda_k \mu_h >0$, we obtain [[$$\begin{aligned}
    &\sum_{i=0}^{m-1} \left \|x_{k,i+1} - x_{\lambda_k}^* \right \|^2 \nonumber \\
                &\leq \left( 1- \frac{\gamma_k \lambda_k \mu_h}{m} \right) \left \|x_k - x_{\lambda_k}^* \right \|^2 + \sum_{i=1}^{m-1} \left \| x_{k,i} - x_{\lambda_k}^* \right \|^2 \nonumber \\& + 2 m \gamma_k^2 (C_f^2 + \lambda_k^2 C_h^2) 
    - 2\gamma_k \sum_{i=0}^{m-1} \left( f_{i+1}(x_{k,i}) + \frac{\lambda_k}{m} h(x_{k,i}) \right) \nonumber\\&+ 2\gamma_k \left(f(x_{\lambda_k}^*)+ \lambda_k h(x_{\lambda_k}^*) \right), \label{ineq1 INC}
    \end{aligned}$$]{}]{} where we used the definition of function $f$ in the second inequality. Now by rearranging the terms and adding and subtracting $f(x_k)+\lambda_k h(x_k)$ we obtain [[$$\begin{aligned}
    &\left \|x_{k+1} - x_{\lambda_k}^* \right \|^2 \\ &\leq \left( 1- \frac{\gamma_k \lambda_k \mu_h}{m} \right) \left \|x_k - x_{\lambda_k}^* \right \|^2 + 2 m \gamma_k^2 (C_f^2 + \lambda_k^2 C_h^2)\\
    & - 2\gamma_k \sum_{i=0}^{m-1} \left((f_{i+1}(x_{k,i}) - f_{i+1}(x_{k}) ) + \frac{\lambda_k}{m} (h(x_{k,i})  -h(x_{k}))  \right) \\& + 2\gamma_k \underbrace{\left(f(x_{\lambda_k}^*)+ \lambda_k h(x_{\lambda_k}^*) -f(x_k)-\lambda_k h(x_k) \right)}_{Term1}\\
    &\leq \left( 1- \frac{\gamma_k \lambda_k \mu_h}{m} \right) \left \|x_k - x_{\lambda_k}^* \right \|^2 + 2 m \gamma_k^2 (C_f^2 + \lambda_k^2 C_h^2)\\
    & + 2\gamma_k \sum_{i=0}^{m-1} \left( \underbrace{|f_{i+1}(x_{k,i}) - f_{i+1}(x_{k}) |}_{Term2} + \frac{\lambda_k}{m} \underbrace{|h(x_{k,i})  -h(x_{k})|}_{Term3}  \right),
    \end{aligned}$$]{}]{} where $Term1\leq 0$ is used due to optimality of $x_{\lambda_k}^*$ for $f+\lambda_k h$. Also, from Remark \[bound on gradient by beck\](b) we know that $Term2 \leq C_f\|x_{k,i} - x_k\|$ and $Term3 \leq C_h\|x_{k,i} - x_k\|$. So, We have $$\begin{aligned}
\label{incrementalineq1}
    &\left \|x_{k+1} - x_{\lambda_k}^* \right \|^2  \leq \left( 1- \frac{\gamma_k \lambda_k \mu_h}{m} \right) \left \|x_{k} - x_{\lambda_k}^* \right \|^2  \nonumber \\ &+ 2 m \gamma_k^2 (C_f^2 + \lambda_k^2 C_h^2)
    + 2 (C_f+\lambda_k C_h) \gamma_k \sum_{i=0}^{m-1} \| x_{k,i} - x_k\|.
    \end{aligned}$$ Next, we find an upper bound for $\| x_{k,i} - x_k\|$. We have\
$    \| x_{k,1} - x_k\|\\ =\left \|\mathcal{P}_X \left(x_{k,0} - \gamma_k\left(g_{f_1}(x_{k,0}) + \frac{\lambda_k}{m} g_h(x_{k,0})\right) \right)- \mathcal{P}_X(x_k) \right \|\\
    \leq \gamma_k \left \|g_{f_1}(x_{k,0}) + \frac{\lambda_k}{m} g_h(x_{k,0})\right \| \leq \gamma_k \left( C_f + \frac{\lambda_k}{m} C_h \right).\\
    $ For $i>0$, in a similar way, we have $$\begin{aligned}
    \| x_{k,i+1} - x_k\|         \leq \| x_{k,i} - x_k\|+ \gamma_k \left( C_f + \frac{\lambda_k}{m} C_h \right).
    \end{aligned}$$ So for $i=0,1, \cdots, m-1$, we have $$\begin{aligned}
 \label{ineq2 INC}
    \| x_{k,i+1} - x_k\| &\leq (i+1)\gamma_k \left( C_f + \frac{\lambda_k}{m} C_h \right) \nonumber \\&\leq  (i+1)\gamma_k \left( C_f + \lambda_k C_h \right).
    \end{aligned}$$ Combining this with , we will obtain $$\begin{aligned}
\label{incrementalineq2}
    \left \|x_{k+1} - x_{\lambda_k}^* \right \|^2          &\leq \left( 1- \frac{\gamma_k \lambda_k \mu_h}{m} \right) \left \|x_{k} - x_{\lambda_k}^* \right \|^2 \nonumber \\ &+ 6 m^2 \gamma_k^2 (C_f^2 + \lambda_k^2 C_h^2).
    \end{aligned}$$ Next, we relate $x_k$ to $x_{\lambda_{k-1}}^*$. We have $$\begin{aligned}
    &\left \|x_{k} - x_{\lambda_k}^* \right \|^2=\left \|x_{k} - x_{\lambda_{k-1}}^* \right\|^2 + \left \|x_{\lambda_k}^* - x_{\lambda_{k-1}}^* \right\|^2 \\ & + \underbrace{2 \left(x_{k} - x_{\lambda_{k-1}}^* \right)^T \left(x_{\lambda_{k-1}}^* - x_{\lambda_k}^*  \right)}_{Term4}.
    \end{aligned}$$ Applying the fact that $2a^Tb \leq \|a\|^2/\alpha + \alpha \|b\|^2 $ for all $a,b \in \mathbb{R}^n$ and $\alpha>0$ for $Term4$ when $\alpha = 2m/{\gamma_k \lambda_k \mu_h}$, we obtain $$\begin{aligned}
    &\left \|x_{k} - x_{\lambda_k}^* \right \|^2 =\left \|x_{k} - x_{\lambda_{k-1}}^* \right\|^2 + \left \|x_{\lambda_k}^* - x_{\lambda_{k-1}}^* \right\|^2 \\ & + \frac{\gamma_k \lambda_k \mu_h}{2m} \left \|x_{k} - x_{\lambda_{k-1}}^* \right\|^2 + \frac {2m} {\gamma_k \lambda_k \mu_h} \left \|x_{\lambda_k}^* - x_{\lambda_{k-1}}^* \right\|^2\\
    & = \left( 1+ \frac{\gamma_k \lambda_k \mu_h}{2m} \right) \left \|x_{k} - x_{\lambda_{k-1}}^* \right\|^2 + \left(1+ \frac{2m}{\gamma_k \lambda_k \mu_h}\right)\left \|x_{\lambda_k}^* - x_{\lambda_{k-1}}^* \right\|^2. 
    \end{aligned}$$ Using Lemma \[lemma:err\_inequalitygeneral\](a), we obtain $$\begin{aligned}
    &\left \|x_{k} - x_{\lambda_k}^* \right \|^2  \leq  \left( 1+ \frac{\gamma_k \lambda_k \mu_h}{2m} \right) \left \|x_{k} - x_{\lambda_{k-1}}^* \right\|^2 \\ & + \left(1+ \frac{2m}{\gamma_k \lambda_k \mu_h}\right) \frac{C_h^2}{\mu_h^2}\left  |1-\frac{\lambda_{k-1}}{\lambda_k} \right|^2.
    \end{aligned}$$ Plugging this inequality into we obtain $$\begin{aligned}
    &\left \|x_{k+1} - x_{\lambda_k}^* \right \|^2  \leq \left( 1- \frac{\gamma_k \lambda_k \mu_h}{m} \right) \left( 1+ \frac{\gamma_k \lambda_k \mu_h}{2m} \right) \left \|x_{k} - x_{\lambda_{k-1}}^* \right \|^2 \\& +\left( 1- \frac{\gamma_k \lambda_k \mu_h}{m} \right) \left(1+ \frac{2m}{\gamma_k \lambda_k \mu_h}\right) \frac{C_h^2}{\mu_h^2}\left    |1-\frac{\lambda_{k-1}}{\lambda_k} \right|^2 \\& +  6 m^2 \gamma_k^2 (C_f^2 + \lambda_k^2 C_h^2).     \end{aligned}$$ The desired result is obtained [[from]{}]{} $0<\mu_k\lambda_k \mu_h \leq 2m$.

We will make use of the following result in Proposition \[conv x\_k INC\].

\[lemma conv lemma\] Let $\{\nu_k\}$ be a sequence of nonnegative scalars and let $\{\alpha_k\}$ and $\{\beta_k\}$ be scalar sequences such that: $$\begin{aligned}
    &\nu_{k+1} \leq (1-\alpha_k)\nu_k + \beta_k \qquad \hbox{for all } k \geq 0,\\
    &0 \leq \alpha_k \leq 1, \ \beta_k \geq0, \ \sum_{k=0}^{\infty}\alpha_k=\infty, \ \sum_{k=0}^{\infty}\beta_k <\infty,  \ \lim_{k \to \infty} \frac{\beta_k}{\alpha_k}=0.
    \end{aligned}$$ Then, $\lim_{k \to \infty} \nu_k = 0$.

\[assumptions for conv INC\] Assume that for all $k \geq 0$ we have $(a)\{\gamma_k\}$ and $\{\lambda_k\}$ are non-increasing positive sequences with $\gamma_0\lambda_0 \leq \frac{2m}{\mu_h}.$\
$(b) \sum_{k=0}^{\infty}\gamma_k\lambda_k= \infty.$ $ \qquad (c) \sum_{k=0}^{\infty}\frac{1}{\gamma_k\lambda_k}\left(\frac{\lambda_{k-1}}{\lambda_k}-1\right)^2 <\infty.$\
$(d) \sum_{k=0}^{\infty} \gamma_k^2 < \infty.$\
$(e) \lim_{k\to \infty} \frac{1}{\gamma_k^2\lambda_k^2}\left(\frac{\lambda_{k-1}}{\lambda_k}-1\right)^2=0.$ $\qquad (f) \lim_{k\to \infty} \frac{\gamma_k}{\lambda_k}=0.$

\[conv x\_k INC\] Consider problem . Let Assumption \[assum:properties\] and \[assumptions for conv INC\] hold and $\{x_k\}$ be generated by Algorithm \[algorithm:IRIG\]. Then,

-   $\lim_{k \rightarrow \infty} \|x_k-x_{\lambda_{k-1}}^*\|^2=0$.

-   If $\lim_{k \rightarrow \infty} \lambda_k=0$, $x_k$ converges to the unique optimal solution of problem , i.e., $x_h^*$.

\(a) Consider the result from Lemma \[lemma: an error bound INC\]. We let $$\begin{aligned}
    \nu_k &\triangleq \|x_k - x_{\lambda_{k-1}}^*\|^2, \ \alpha_k \triangleq \frac{\gamma_k\lambda_k \mu_h}{2m}\\ \beta_k &\triangleq \frac{3m C_h^2}{\gamma_k \lambda_k \mu_h^3}\left(1-\frac{\lambda_{k-1}}{\lambda_k} \right)^2 +  6 m^2 \gamma_k^2 (C_f^2 + \lambda_k^2 C_h^2).
    \end{aligned}$$ From Assumption \[assumptions for conv INC\](a,b,c), since $\{\gamma_k\}$ and $\{\lambda_k\} $ are positive and $\gamma_0\lambda_0 \leq 2m/\mu_h$, we have $0 \leq \alpha_k \leq 1$ and $ \beta_k \geq0$ and also $\sum_{k=1}^{\infty}\alpha_k=\infty$ and $\sum_{k=1}^{\infty}\beta_k <\infty$. To show that all the necessary assumptions for Lemma \[lemma conv lemma\] are satisfied, we have $$\begin{aligned}
    \lim_{k \to \infty} \frac{\beta_k}{\alpha_k}    &=\frac{6m^2 C_h^2}{\mu_h^4}\lim_{k \to \infty} \frac{1}{\gamma_k^2\lambda_k^2}\left(\frac{\lambda_{k-1}}{\lambda_k}-1\right)^2 +\frac{12m^3C_f^2}{\mu_h}\lim_{k\to \infty} \frac{\gamma_k}{\lambda_k} \\& +
    \frac{12m^3C_h^2}{\mu_h}\lim_{k\to \infty} \gamma_k\lambda_k.
    \end{aligned}$$ Considering Assumption \[assumptions for conv INC\](e,f), we only need to show that $\lim_{k\to \infty} \gamma_k\lambda_k=0$. Since $\{\lambda_k\}$ is non-increasing for all $k\geq 0$, we have $\lambda_0^2 {\gamma_k}/{\lambda_k} \geq \gamma_k\lambda_k$. So by Assumption \[assumptions for conv INC\](f), $\lim_{k \to \infty} \gamma_k\lambda_k=0$. Consequently $\lim_{k \to \infty} \frac{\beta_k}{\alpha_k}=0$. Now Lemma \[lemma conv lemma\] can be applied. We have $$\begin{aligned}
    \lim_{k \to \infty} \nu_k = \lim_{k \to \infty} \|x_k - x_{\lambda_{k-1}}^*\|^2=0.
    \end{aligned}$$

\(b) Applying the triangular inequality, we obtain $$\begin{aligned}
    \|x_k -x_h^*\|^2 \leq 2\|x_k - x_{\lambda_{k-1}}^* \|^2 + 2 \|x_{\lambda_{k-1}}^*-x_h^*\|^2, \ \hbox{for all } k \geq 0.
    \end{aligned}$$ From part (a), $\|x_k -x_h^*\|^2$ converges to zero. Also, from Lemma \[lemma:err\_inequalitygeneral\](b), we know that when $\lambda_k \rightarrow 0$ the sequence $\{x_{\lambda_k}^*\}$ converges to the unique optimal solution of problem , i.e., $x_h^*$. Therefore the result holds.

To have previous proposition work, we require that sequences $\{\gamma_k\}$ and $\{\lambda_k\}$ satisfy Assumption \[assumptions for conv INC\]. Below, we provide a set of feasible sequences for this assumption. The proof is analogous to that of Lemma 5 in [@Yousefian17].

\[lemma condition for sequences INC\] Assume $\{\gamma_k\}$ and $\{\lambda_k\}$ are sequences such that $\gamma_k=\frac{\gamma_0}{(k+1)^a}$ and $\lambda_k=\frac{\lambda_0}{(k+1)^b}$ where $a,b,\gamma_0$ and $\lambda_0$ are positive scalars and $\gamma_0\lambda_0 \leq \frac{2m}{\mu_h}$. If $a>b$, $a>0.5$ and $a+b<1$, then the sequences $\{\gamma_k\}$ and $\{\lambda_k\}$ satisfy Assumption \[assumptions for conv INC\].

The following is a useful lemma in proving convergence that we will apply in Theorem \[thm conv for xbar INC\].

\[lemma thm for avr\] Let $\{u_t\}\subset \mathbb{R}^n$ be a convergent sequence with the limit point $\hat u\in\mathbb{R}^n$ and let $\{\alpha_k\}$ be a sequence of positive numbers where $\sum_{k=0}^\infty \alpha_k=\infty$. Suppose $\{v_k\}$ is given by $v_k\triangleq \left({\sum_{t=0}^{k-1} \alpha_t u_t}\right)/{\sum_{t=0}^{k-1} \alpha_t}$ for all $k\ge1$. Then, $\lim_{k \rightarrow \infty} v_k=\hat u.$

Now, we can illustrate our ultimate goal in this section which is showing the convergence of the sequence $\{\bar{x_k}\}$ generated by Algorithm \[algorithm:IRIG\] to $x_h^*$

\[thm conv for xbar INC\] Consider problem . Let Assumption \[assum:properties\] hold. Also assume $\{\gamma_k\}$ and $\{\lambda_k\}$ are sequences such that $\gamma_k=\frac{\gamma_0}{(k+1)^a}$ and $\lambda_k=\frac{\lambda_0}{(k+1)^b}$ where $a,b,\gamma_0$ and $\lambda_0$ are positive scalars and $\gamma_0\lambda_0 \leq \frac{2m}{\mu_h}$. Let $\{\bar x_k\}$ be generated by Algorithm \[algorithm:IRIG\]. If $a>b$, $a>0.5$, $a+b<1$ and $ar \leq 1$, then $\{\bar x_k\}$ converges to $x_h^*$.

Considering the given assumptions, by Lemma \[lemma condition for sequences INC\] we can see that Assumption \[assumptions for conv INC\] holds. We have $$\begin{aligned}
    \|\bar x_{k+1}-x_h^*\|=\left\|\sum_{t=0}^{k} \psi_{t,k} x_t-\sum_{t=0}^{k}\psi_{t,k} x_h^*\right\|\leq \sum_{t=0}^{k} \psi_{t,k}\|x_t-x_h^*\|,
        \end{aligned}$$ applying $\sum_{t=0}^{k} \psi_{t,k}=1$ from and the triangular inequality. Now consider the definition of $\psi_k$ and let $\alpha_t \triangleq \gamma_t^r$, $u_t \triangleq \|x_t-x_h^*\|$ and $v_{k+1} \triangleq \sum_{t=0}^{k} \psi_{t,k} \|x_t-x_h^*\|$. Since $ar \leq1$ we have $\sum_{t=0}^\infty \alpha_t=\sum_{t=0}^\infty \gamma_t^{r}=\sum_{t=0}^\infty (t+1)^{-ar}=\infty$. The sequence $\{\lambda_t \}$ is decreasing to zero due to $b>0$, So, from Proposition \[conv x\_k INC\](b), $u_t=\|x_t-x_h^*\|$ converges to zero. Therefore, for $\hat u=0$ we can apply Lemma \[lemma thm for avr\] and thus, $\|\bar x_{k+1}-x_h^*\|$ converges to zero.

Rate analysis {#sec:rate result}
=============

In this section, we first find an error bound with respect to the optimal values of the lower level function $f$ which indeed shows the feasibility of [[the problem]{}]{} . Then, we apply it to derive a convergence rate for the algorithm.

\[lemma averaging\] Consider the sequence $\{\bar x_N\}$ generated by Algorithm \[algorithm:IRIG\]. Let Assumption \[assum:properties\] hold and $\{\gamma_k\}$ and $\{\lambda_k\}$ be positive and non-increasing sequences. Then, for all $N\geq 1$ and $z\in X$ we have $$\begin{aligned}
    &f(\bar x_N)-f^*\leq \left(\sum_{k=0}^{N-1} \gamma_k^r\right)^{-1} \left( m \sum_{k=0}^{N-1} \gamma_k^{r+1} (C_f^2 + \lambda_k^2 C_h^2) \right. \cr\\& \left. 
    +m^2 C_f \sum_{k=0}^{N-1} \gamma_k^{r+1} \left( C_f + \lambda_k C_h \right) + 2 M_h  \sum_{k=0}^{N-1} \gamma_k^r \lambda_k + 2M^2\gamma_{N-1}^{r-1} \right),
    \end{aligned}$$ where $M_h, M$ are scalars such that $\|h(x)\|\leq M_h$, $\|x\| \leq M$ for all $x \in X$.

[[Similar to]{}]{} relation , we [[can have]{}]{} $$\begin{aligned}
    &\left \|x_{k+1} - x^* \right \|^2  \leq \left \|x_k - x^* \right \|^2 + 2 m \gamma_k^2 (C_f^2 + \lambda_k^2 C_h^2) \\
    & - 2\gamma_k \sum_{i=0}^{m-1} \left( f_{i+1}(x_{k,i}) + \frac{\lambda_k}{m} h(x_{k,i}) \right) + 2\gamma_k \left(f^*+ \lambda_k h(x^*) \right) \\
    & \leq \left \|x_k - x^* \right \|^2 + 2 m \gamma_k^2 (C_f^2 + \lambda_k^2 C_h^2) \\
    & - 2\gamma_k \sum_{i=0}^{m-1} f_{i+1}(x_{k,i}) + 2\gamma_k f^*+ 4\gamma_k \lambda_k M_h.
    \end{aligned}$$ Adding and subtracting $f(x_k)$, we obtain $$\begin{aligned}
    &\left \|x_{k+1} - x^* \right \|^2  \leq \left \|x_k - x^* \right \|^2 + 2 m \gamma_k^2 (C_f^2 + \lambda_k^2 C_h^2) \\
    & - 2\gamma_k \sum_{i=0}^{m-1} \left( f_{i+1}(x_{k,i}) - f_i(x_k) \right) + 2\gamma_k \left( f^* - f(x_k)\right)+ 4\gamma_k \lambda_k M_h\\
    & \leq \left \|x_k - x^* \right \|^2 + 2 m \gamma_k^2 (C_f^2 + \lambda_k^2 C_h^2) \\
    & + 2\gamma_k \sum_{i=0}^{m-1} \left| f_{i+1}(x_{k,i}) - f_i(x_k) \right| + 2\gamma_k \left( f^* - f(x_k)\right)+ 4\gamma_k \lambda_k M_h.
    \end{aligned}$$ Applying Remark \[bound on gradient by beck\](b), $| f_{i+1}(x_{k,i}) - f_i(x_k)| \leq C_f\|x_{k,i} - x_k\|$, we have $$\begin{aligned}
    &\left \|x_{k+1} - x^* \right \|^2  \leq \left \|x_k - x^* \right \|^2 + 2 m \gamma_k^2 (C_f^2 + \lambda_k^2 C_h^2) \\
    & + 2 C_f \gamma_k \sum_{i=0}^{m-1} \|x_{k,i} - x_k\| + 2\gamma_k \left( f^* - f(x_k)\right)+ 4\gamma_k \lambda_k M_h.
    \end{aligned}$$ Using the inequality , we obtain $$\begin{aligned}
    &\left \|x_{k+1} - x^* \right \|^2  -  \left \|x_k - x^* \right \|^2  \leq 2\gamma_k \left( f^* - f(x_k)\right) \nonumber\\
    & + 2 m \gamma_k^2 (C_f^2 + \lambda_k^2 C_h^2)   + 2 m^2 C_f \gamma_k^2 \left( C_f + \lambda_k C_h \right) + 4\gamma_k \lambda_k M_h.\label{ineq3 INC} 
    \end{aligned}$$ Multiplying both sides by $\gamma_k^{r-1}$ and adding and subtracting $\gamma_{k-1}^{r-1} \|x_k - x^*\|^2$ to the left hand side, we have $$\begin{aligned}
    &\gamma_k^{r-1} \left \|x_{k+1} - x^* \right \|^2  - \gamma_{k-1}^{r-1} \|x_k - x^*\|^2  + \left( \gamma_{k-1}^{r-1}-\gamma_k^{r-1}\right)  \left \|x_k - x^* \right \|^2 \\& \leq 2 \gamma_k^r \left( f^* - f(x_k)\right) + 2 m \gamma_k^{r+1} (C_f^2 + \lambda_k^2 C_h^2) \\
    & + 2 m^2 C_f \gamma_k^{r+1} \left( C_f + \lambda_k C_h \right) + 4 \gamma_k^r \lambda_k M_h.
    \end{aligned}$$ Since $\{\gamma_k\}$ in non-increasing and $r<1$ we have $ \gamma_{k-1}^{r-1} \leq \gamma_k^{r-1} $. Also, by the triangle inequality $\left \|x_k - x^* \right \|^2 \leq 2\|x_k\|^2+2\|x^*\|^2 \leq 4M^2$. So, we obtain $$\begin{aligned}
    &\gamma_k^{r-1} \left \|x_{k+1} - x^* \right \|^2  - \gamma_{k-1}^{r-1} \|x_k - x^*\|^2  + 4M^2\left( \gamma_{k-1}^{r-1}-\gamma_k^{r-1}\right)\\ & \leq 2 \gamma_k^r \left( f^* - f(x_k)\right) + 2 m \gamma_k^{r+1} (C_f^2 + \lambda_k^2 C_h^2) \\
    & + 2 m^2 C_f \gamma_k^{r+1} \left( C_f + \lambda_k C_h \right) + 4 \gamma_k^r \lambda_k M_h.
    \end{aligned}$$ Taking summation over $k=1,2, \cdots, N-1$, we obtain $$\begin{aligned}
    &\gamma_{N-1}^{r-1} \left \|x_N - x^* \right \|^2  - \gamma_{0}^{r-1} \|x_1 - x^*\|^2 + 4M^2\left( \gamma_{0}^{r-1}-\gamma_{N-1}^{r-1}\right)   \\& \leq 2 \sum_{k=1}^{N-1} \gamma_k^r \left( f^* - f(x_k)\right)  \label{ineq4 INC}+ 2 m \sum_{k=1}^{N-1} \gamma_k^{r+1} (C_f^2 + \lambda_k^2 C_h^2) 
    \\ & + 2 m^2 C_f \sum_{k=1}^{N-1} \gamma_k^{r+1} \left( C_f + \lambda_k C_h \right) + 4 M_h \sum_{k=1}^{N-1} \gamma_k^r \lambda_k .
    \end{aligned}$$ Removing non-negative terms from the left-hand side of the preceding inequality, we have $$\begin{aligned}
    &- \gamma_{0}^{r-1} \|x_1 - x^*\|^2  - 4M^2\gamma_{N-1}^{r-1}  \leq 2 \sum_{k=1}^{N-1} \gamma_k^r \left( f^* - f(x_k)\right) \nonumber \\ & + 2 m \sum_{k=1}^{N-1} \gamma_k^{r+1} (C_f^2 + \lambda_k^2 C_h^2) 
    + 2 m^2 C_f \sum_{k=1}^{N-1} \gamma_k^{r+1} \left( C_f + \lambda_k C_h \right) \nonumber \\& + 4 M_h \sum_{k=1}^{N-1} \gamma_k^r \lambda_k . 
    \end{aligned}$$ From for $k=0$, we obtain $$\begin{aligned}
    \left \|x_1 - x^* \right \|^2  & \leq 2\gamma_0 \left( f^* - f(x_0)\right) + 2 m \gamma_0^2 (C_f^2 + \lambda_0^2 C_h^2)  \\
    & + 2 m^2 C_f \gamma_0^2 \left( C_f + \lambda_0 C_h \right) + 4\gamma_0 \lambda_0 M_h + 4M^2. 
    \end{aligned}$$ By multiplying both sides of the preceding inequality by $\gamma_{0}^{r-1}$ and summing it with the relation , we obtain $$\begin{aligned}
    &- 4M^2\gamma_{N-1}^{r-1}  \leq 2 \sum_{k=0}^{N-1} \gamma_k^r \left( f^* - f(x_k)\right)  + 4 M_h \sum_{k=0}^{N-1} \gamma_k^r \lambda_k \\& + 2 m \sum_{k=0}^{N-1} \gamma_k^{r+1} (C_f^2 + \lambda_k^2 C_h^2) 
     + 2 m^2 C_f \sum_{k=0}^{N-1} \gamma_k^{r+1} \left( C_f + \lambda_k C_h \right)  . 
    \end{aligned}$$ Rearranging the terms we have $$\begin{aligned}
    &2 \sum_{k=0}^{N-1} \gamma_k^r \left( f(x_k) - f^*\right)  \leq  2 m \sum_{k=0}^{N-1} \gamma_k^{r+1} (C_f^2 + \lambda_k^2 C_h^2) \\& 
    + 2 m^2 C_f \sum_{k=0}^{N-1} \gamma_k^{r+1} \left( C_f + \lambda_k C_h \right) + 4 M_h \sum_{k=0}^{N-1} \gamma_k^r \lambda_k + 4M^2\gamma_{N-1}^{r-1} . 
    \end{aligned}$$ Now we divide both sides by $2\sum_{k=0}^{N-1} \gamma_k^r$ and use the definition of $\psi_{k,N-1}$ in , $$\begin{aligned}
    &\sum_{k=0}^{N-1} \psi_{k,N-1} \left( f(x_k) - f^*\right) \\& \leq   \left(\sum_{k=0}^{N-1} \gamma_k^r\right)^{-1} \left( m \sum_{k=0}^{N-1} \gamma_k^{r+1} (C_f^2 + \lambda_k^2 C_h^2) \right. \cr\\& \left. 
    +m^2 C_f \sum_{k=0}^{N-1} \gamma_k^{r+1} \left( C_f + \lambda_k C_h \right) + 2 M_h  \sum_{k=0}^{N-1} \gamma_k^r \lambda_k + 2M^2\gamma_{N-1}^{r-1} \right). 
    \end{aligned}$$ We know $\sum_{k=0}^{N-1} \psi_{k,N-1}=1$ also $f(\bar x_N) \leq \sum_{k=0}^{N-1} \psi_{k,N-1}f(x_k)$ because of convexity of $f$. So, we obtain the desired result.

The following lemma, will be used to find a convergence rate statement in Theorem \[thm rate\].

\[Lemma 9, page 418 in [@Yousefian17]\]\[lemma:ineqHarmonic\] For any scalar $\alpha\neq -1$ and integers $\ell$ and $N$ where $0\leq \ell \leq N-1$, we have $$\begin{aligned}
    \frac{N^{\alpha+1}-(\ell+1)^{\alpha+1}}{\alpha+1}&\leq \sum_{k=\ell}^{N-1}(k+1)^\alpha \\& \leq (\ell+1)^\alpha+\frac{(N+1)^{\alpha+1}-(\ell+1)^{\alpha+1}}{\alpha+1}.
    \end{aligned}$$

In the following theorem[[,]{}]{} we present a rate statement for Algorithm \[algorithm:IRIG\].

\[**[A rate statement for Algorithm \[algorithm:IRIG\]]{}\] \[thm rate\] Assume $\{\bar x_N\}$ is generated by Algorithm \[algorithm:IRIG\] to solve problem . Let Assumption \[assum:properties\] and \[assumptions for conv INC\] hold and also $0<\epsilon<0.5$ and $r<1$ be arbitrary constants. Assume for $0<\epsilon<0.5$, $\{\gamma_k\}$ and $\{\lambda_k\}$ are sequences defined as $$\gamma_k=\frac{\gamma_0}{(k+1)^{0.5+0.5\epsilon}} \hbox{ and } \lambda_k=\frac{\lambda_0}{(k+1)^{0.5-\epsilon}},$$ such that $\gamma_0$ and $\lambda_0$ are positive scalars and $\gamma_0\lambda_0 \mu_h \leq 2m$. Then,\
(a) The sequence $\{\bar x_N\}$ converges to the unique optimal solution of problem , i.e., $x_h^*$.\
(b) $f(\bar x_N)$ converges to $f^*$ with the rate ${\cal O}\left(1/N^{0.5-\epsilon}\right)$.**

Throughout, we set $a:=0.5+0.5\epsilon$, $b:=0.5-\epsilon$.

\(a) From the values of $a$ and $b$, and that $r<1$ and $0<\epsilon<0.5$, we have $$\begin{aligned}
    a>b>0,\ a>0.5, \ a+b=1-0.5\epsilon<1,\\ ar= 0.5(1+\epsilon)r <0.5(1.5)=0.75<1.
    \end{aligned}$$ This implies that all conditions of Theorem \[thm conv for xbar INC\] are satisfied. Therefore, $\{\bar x_N\}$ converges to $x_h^*$ almost surely.

\(b) Since $\{\lambda_k\}$ is a non-increasing sequence from Lemma \[lemma averaging\] we have $$\begin{aligned}
&f(\bar x_N) - f^*  \leq   \left(\sum_{k=0}^{N-1} \gamma_k^r\right)^{-1} \left( m \left(C_f^2 + \lambda_0^2 C_h^2\right) \sum_{k=0}^{N-1} \gamma_k^{r+1}  \right. \cr\\& \left. 
+m^2 C_f\left( C_f + \lambda_0 C_h \right) \sum_{k=0}^{N-1} \gamma_k^{r+1} + 2 M_h   \sum_{k=0}^{N-1} \lambda_k \gamma_k^r + 2M^2\gamma_{N-1}^{r-1} \right). \end{aligned}$$ We have $\gamma_k=\gamma_0/(k+1)^a$ and $\lambda_k=\lambda_0/(k+1)^b$, thus $$\begin{aligned}
& f(\bar x_N) - f^*  \\& \leq   \left(\sum_{k=0}^{N-1} \frac{\gamma_0^r}{(k+1)^{ar}}\right)^{-1} \left( m \left(C_f^2 + \lambda_0^2 C_h^2\right) \sum_{k=0}^{N-1} \frac{\gamma_0^{r+1}}{(k+1)^{a(r+1)}}  \right. \cr\\& \left.  
+m^2 C_f\left( C_f + \lambda_0 C_h \right) \sum_{k=0}^{N-1} \frac{\gamma_0^{r+1}}{(k+1)^{a(r+1)}} + 2 M_h \sum_{k=0}^{N-1} \frac{\lambda_0 \gamma_0^r}{(k+1)^{ar+b}}  \right. \cr\\& \left. + 2M^2\gamma_0^{r-1}N^{a(1-r)} \right). \end{aligned}$$ Rearranging the terms, we have [[$$\begin{aligned}
&f(\bar x_N) - f^* \leq   \left(\sum_{k=0}^{N-1} \frac{\gamma_0^r}{(k+1)^{ar}}\right)^{-1} \\& \times\left(  2 M_h \sum_{k=0}^{N-1} \frac{\lambda_0 \gamma_0^r}{(k+1)^{ar+b}} + 2M^2\gamma_0^{r-1}N^{a(1-r)}  \right. \cr\\& \left. 
+\left( m^2 C_f\left( C_f + \lambda_0 C_h \right)+ m \left(C_f^2 + \lambda_0^2 C_h^2\right) \right) \sum_{k=0}^{N-1} \frac{\gamma_0^{r+1}}{(k+1)^{a(r+1)}}  \right). \end{aligned}$$]{}]{} Let us define $$\begin{aligned}
&\hbox{ Term1}= \left(\sum_{k=0}^{N-1} \frac{1}{(k+1)^{ar}}\right)^{-1} N^{a(1-r)},\\
&\hbox{Term2}=\left(\sum_{k=0}^{N-1} \frac{1}{(k+1)^{ar}}\right)^{-1}\left(\sum_{k=0}^{N-1}\frac{1}{(k+1)^{ar+b}}\right).\end{aligned}$$ We have $$\begin{aligned}
&\hbox{Term1} \leq\frac{N^{a(1-r)}}{\frac{N^{1-ar}-1}{1-ar}}=\frac{(1-ar)N^{a(1-r)}}{N^{1-ar}-1}={\cal O}\left(N^{-(1-a)}\right), \\
&\hbox{Term2} \leq\frac{\frac{(N+1)^{1-ar-b}-1}{1-ar-b}+1}{\frac{N^{1-ar}-1}{1-ar}} =\frac{(1-ar)\left((N+1)^{1-ar-b}-1\right)}{(1-ar-b)\left(N^{1-ar}-1\right)} \\& + \frac{1-ar}{N^{1-ar}-1}  ={\cal O}\left(N^{-b}\right) + {\cal O}\left(N^{-(1-ar)}\right).\end{aligned}$$ So, we have $$\begin{aligned}
f(\bar x_N)-f^*\leq {\cal O}\left(N^{-\min \{1-ar,1-a,b\}}\right)= {\cal O}\left(N^{-\min \{1-a,b\}}\right),\end{aligned}$$ where we used $1-a\leq 1-ar$. Replacing $a$ and $b$ by their values, we have $$\begin{aligned}
f(\bar x_N)-f^*\leq {\cal O}\left(N^{-\min \{0.5-0.5\epsilon,0.5-\epsilon\}}\right)={\cal O}\left(N^{-(0.5-\epsilon)}\right).\end{aligned}$$

[c| c c c]{} & $x_0=-10\times\mathbf{1}_n$ & $x_0=\mathbf{0}_n$ & $x_0=10\times\mathbf{1}_n$\
\

[90]{} $(10,1)$

&

![image](Accfig1.pdf)

&

![image](Accfig2.pdf)

&

![image](Accfig3.pdf)

\

[90]{} $(1,10)$

&

![image](Accfig4.pdf)

&

![image](Accfig5.pdf)

&

![image](Accfig6.pdf)

\

[90]{} $(0.1,0.1)$

&

![image](Accfig7.pdf)

&

![image](Accfig8.pdf)

&

![image](Accfig9.pdf)

\[fig:fiveplots\]

Numerical results {#sec:num}
=================

In this section, we apply the IR-IG method on a text classification problem. In this problem we assume to have a summation of hinge loss function, i.e., $\mathcal{L}(\langle x,a\rangle, b) \triangleq \max \{0, 1-b \langle x,a\rangle\}$ for any sample $a,b$ in the lower level of . The set of observations $(a_i,b_i)$ are derived from the Reuters Corpus Volume I (RCV1) dataset (see [@Lewis04]). This dataset has categorized Reuters articles, from 1996 to 1997, into four groups: *Corporate/Industrial, Economics, Government/Social* and *Markets*. In this application, we use a subset of this dataset with $N=50,000$ articles and 138,921 tokens to perform a binary classification of articles only with respect to the *Markets* class. Each vector $a_i$ represents the existence of all tokens in article $i$, and $b_i$ shows whether the article belongs to the *Markets* class. To decide that a new article can be placed in the *Markets* class, the problem in the lower level should be solved such that the optimal solution is a weight vector of the tokens regarding the *Markets* class which minimizes the total loss. However, to make such a decision, an optimal solution with a large number of nonzero components are undesirable. [[To induce sparsity,]{}]{} we consider the [[following]{}]{} bilevel problem:

$$\begin{aligned}
 \label{numeric upper2}
\displaystyle \mbox{minimize}& \ h(x)\triangleq  \frac{\mu_h }{2}\|x\|_2^2+ \|x\|_1\\
\mbox{subject to} & \ x \in \arg\min_{y \in \mathbb{R}^n } \sum_{i=1}^{m}\sum_{j=1}^{N/m}{ \mathcal{L}(a_{(i-1)N/m+j}^Ty, b_{(i-1)N/m+j})}, \nonumber\end{aligned}$$

where, we let $\mu_h=0.1$ and we consider each batch of $N/m=1000$ articles to be one component function with the total number of component functions $m=50$. The function $h$ is strongly convex with parameter $\mu_h$. [[We]{}]{} let $\gamma_k$ and $\lambda_k$ be given by the rules in Theorem \[thm rate\]. We [[study the sensitivity of the method by changing]{}]{} $x_0$, $\gamma_0$, $\lambda_0$, and the averaging parameter $r<1$. We finally report the logarithm of average of the loss function $\mathcal{L}$. The plots in Fig. 1 show the convergence of the IR-IG method for the problem . The results show the convergece of Algorithm \[algorithm:IRIG\] with different initial values such as the starting point or parameters $\gamma_0, \lambda_0$ while when we pick a smaller $r$ the algorithm is faster in all the cases.

Concluding remarks {#sec:rem}
==================

Motivated by the applications of incremental [[gradient schemes]{}]{} in distributed optimization, especially in machine learning and large data training, we develop an iterative regularized incremental first-order method, called IR-IG, for solving a class of bilevel convex optimization. We prove the convergence of IR-IG and establish the corresponding rate in terms of the [[lower]{}]{} level objective function. We finally [[apply IR-IG to]{}]{} a binary text classification problem [[and]{}]{} demonstrate the [[performance]{}]{} of [[the proposed algorithm]{}]{}.

[99]{}

A. Beck, First-order methods in optimization, MOS-SIAM Series on Optimization, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2017.

A. Beck and S. Sabach, A first-order method for finding minimal norm-like solutions of convex optimization problems, Mathematical Programming, 147(2) (2014), 25-46.

A. Beck and M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM journal on imaging sciences, 2(1) (2009), 183-202.

D. Bertsekas, Nonlinear programming, Athena Scientific, Belmont, MA, 1999.

D. Blatt, A. Hero, and H. Gauchman, A convergent incremental gradient method with a constant step-size, SIAM Journal Of Optimization, 18 (2007), 29-51.

L. Bottou and Y. Le Cun, On-line learning for very large data sets, Applied Stochastic Models in Business and Industry, 21 (2005), 137-151.

A. Defazio, F. Bach, and S. Lacoste-Julien, SAGA: A fast incremental gradient method with support for non-strongly convex composite objectives, in Advances in Neural Information Processing Systems, MIT Press, Cambridge, MA, 2014, 1646-1654.

M. P. Friedlander and P. Tseng, Exact regularization of convex programs, SIAM Journal Of Optimization, 18(4) (2007), 1326-1350.

G. Garrigos, L. Rosasco, and S. Villa, Iterative regularization via dual diagonal descent, Journal of Mathematical Imaging and Vision 60(2) (2018), 189-215.

M. G[ü]{}rb[ü]{}zbalaban, A. Ozdaglar, and P. A. Parillo, On the convergence rate of incremental aggregated gradient algorithms, SIAM Journal Of Optimization, 27(2) (2018), 640-660.

S. -P. Han and O. L. Mangasarian, Exact penalty functions in nonlinear programming, 17(1) (1979), 251-269.

K. Knopp, Theory and applications of infinite series, Blackie & Son Ltd., Glasgow, Great Britain, 1951.

M. Kočvara and J. V. Outrata, Optimization problems with equilibrium constraints and their numerical solution, Mathematical Programming, 101(1) (2004), 119-149.

D. D. Lewis, Y. Yang, T. G. Rose, and F. Li, RCV1: A new benchmark collection for text categorization research, Journal of Machine Learning Research, 5 (2004), 361-397.

Z.-Q. Luo, J.-S. Pang and D. Ralph, Mathematical programs with equilibrium constraints, Cambridge University Press, Cambridge (1996).

J. Mairal, Optimization with first-order surrogate functions, in ICML, JMLR Proceedings 28, 2013, 783-791.

O. L. Mangasarian, Sufficiency Of exact penalty minimization, SIAM Journal on Control and Optimization, 23(1) (1985), 30-37.

A. Nedić and D. Bertsekas, Incremental subgradient methods for nondifferentiable optimization, SIAM Journal of Optimization, 12(1) (2001), 109-138.

E. S. H. Neto and Á. A. R. De Pierro, On perturbed steepest descent methods with inexact line search for bilevel convex optimization, Optimization, 60 (2011), 991-1008.

B. T. Polyak, Introduction to optimization, Optimization Software, Inc., New York, 1987.

N. L. Roux, M. Schmidt, and F. R. Bach, A stochastic gradient method with an exponential convergence rate for finite training sets, in Advances in Neural Information Processing Systems 25, F. Pereira, C. J. C. Burges, L. Bottou, and K. Q. Weinberger, eds., MIT Press, Cambridge, MA, 2012, 2663-2671.

S. Sabach and S. Shtern, A first-order method for solving convex bilevel optimization Problems, SIAM Journal on Optimization, 27(2) (2017), 640-660.

M. Solodov, A bundle method for a class of bilevel nonsmooth convex minimization problems, SIAM Journal Of Optimization, 18(1) (2007), 242-259.

M. Solodov, An explicit descent method for bilevel convex optimization, Journal of Convex Analysis, 14(2) (2007).

A. N. Tikhonov and V. Y. Arsenin, Solutions of ill-posed problems, V. H. Winston and Sons, Washington, D. C., 1977. Translated from Russian.

P. Tseng and S. Yun, Incrementally updated gradient methods for constrained and regularized optimization, Journal of Optimization Theory and Applications, 160 (2014), 832-853.

H.-K. Xu, Viscosity approximation methods for nonexpansive mappings, Journal of Mathematical Analysis and Applications, 298(1) (2004), 279-291.

I. Yamada, M. Yukawa, and M. Yamagishi, Minimizing the [[M]{}]{}oreau envelope of nonsmooth convex functions over the fixed point set of certain quasi-nonexpansive mappings, in Fixed-Point Algorithms for Inverse Problems in Science and Engineering, Springer, New York, 2011, 345-390.

F. Yousefian, A. Nedić, and U. V. Shanbhag., On Smoothing, regularization and averaging in stochastic approximation methods for stochastic variational inequality problems, Mathematical Programming, 165 (1) (2017), 391-431.

F. Yousefian, A. Nedić, and U. V. Shanbhag., On stochastic mirror-prox algorithms for stochastic [[C]{}]{}artesian variational inequalities: randomized block coordinate, and optimal averaging schemes, Set-Valued and Variational Analysis, (2018), https://doi.org/10.1007/s11228-018-0472-9.

[^1]: $^{1}$Mostafa Amini and $^{2}$Farzad Yousefian are with the Department of Industrial Engineering and Management, Oklahoma State University, Stillwater, OK 74078, USA, [moamini@okstate.edu, farzad.yousefian@okstate.edu]{}
