# A cut-and-project perspective for linearized Bregman iterations

Yu-Hong Dai Kangkang Deng Hui Zhang LSEC, Academy of Mathematics and Systems
Science, Chinese Academy of Sciences, Beijing, China. Email:
<EMAIL_ADDRESS>of Mathematics, National University of Defense
Technology, Changsha, Hunan 410073, China. Email<EMAIL_ADDRESS>Department of Mathematics, National University of Defense Technology,
Changsha, Hunan 410073, China. Email<EMAIL_ADDRESS>
###### Abstract

The linearized Bregman iterations (LBreI) and its variants are powerful tools
for finding sparse or low-rank solutions to underdetermined linear systems. In
this study, we propose a cut-and-project perspective for the linearized
Bregman method via a bilevel optimization formulation, along with a new
unified algorithmic framework. The new perspective not only encompasses
various existing linearized Bregman iteration variants as specific instances,
but also allows us to extend the linearized Bregman method to solve more
general inverse problems. We provide a completed convergence result of the
proposed algorithmic framework, including convergence guarantees to feasible
points and optimal solutions, and the sublinear convergence rate. Moreover, we
introduce the Bregman distance growth condition to ensure linear convergence.
At last, our findings are illustrated via numerical tests.

Keywords: linearized Bregman iteration, bilevel optimization, Bregman
projections, sparse solution, cutting plane method, growth condition

AMS subject classifications. 90C25, 90C30, 65K05, 49M37

##  1 Introduction

The linearized Bregman iterations (LBreI) method, suggested by Darbon and
Osher [12] and introduced in the seminal work [27], has become a popular
method for finding regularized solutions to the underdetermined linear systems
$Ax=b$, where $A\in\mathbb{R}^{m\times n}$ with $m<n$ and $b\in\mathbb{R}^{m}$
are given. For example, to seek sparse solutions, it solves the following
linearly constrained optimization problem

$\operatorname*{minimize\quad}\lambda\|x\|_{1}+\frac{1}{2}\|x\|^{2},~{}~{}\mbox{subject
to}~{}~{}Ax=b,$ (1.1)

where $\lambda>0$ is a regularization parameter. The corresponding algorithmic
scheme is straightforward and consists of two steps

$\displaystyle x^{*}_{k+1}$ $\displaystyle=$ $\displaystyle
x^{*}_{k}-A^{T}(Ax_{k}-b),$ $\displaystyle x_{k+1}$ $\displaystyle=$
$\displaystyle\mathcal{S}_{\lambda}(x_{k+1}^{*})$

initialized with $x_{0}=x^{*}_{0}=0$, where
$\mathcal{S}_{\lambda}(x)=\min\\{|x|-\lambda,0\\}\textrm{sign}(x)$ is the
component-wise soft shrinkage. At the first step, the dual variable
$x_{k}^{*}$ is updated by performing a gradient descent with
$f(x)=\frac{1}{2}\|Ax-b\|^{2}$ to approach the constraint $Ax=b$; at the
second step, the primal variable $x_{k}$ is updated by applying the soft
shrinkage operator to the dual variable $x_{k}^{*}$ so that a sparse solution
can be obtained. The convergence of the two-step iteration has been analyzed
in [9].

In order to speed up the convergence and extend the LBreI method to solve more
general problems than (1.1), several novel insights into this method have
emerged over the past decade. The first insight, elucidated in the work by
[26], identifies LBreI as a gradient descent method for the dual problem of
(1.1). This identification enables the application of gradient-based
optimization techniques, including line search, Barzilai-Borwein, limited
memory BFGS, nonlinear conjugate gradient, and Nesterov’s methods, to
accelerate the LBreI method. A second critical insight, proposed by [11]
considers the following least square problem

$\begin{array}[]{ccc}&\operatorname*{minimize\quad}&~{}~{}\lambda\|x\|_{1}+\frac{1}{2}\|x\|^{2}{}\\\
&\mbox{subject to}&~{}~{}x\in\arg\min_{x}\frac{1}{2}\|Ax-b\|^{2},\end{array}$
(1.2)

where the measurement data $b$ can be out the range of $A$. In particular, if
$b$ is in the range of $A$, then the inner level optimization reduces to
$Az=b$ hence (1.2) reduces to (1.1). The authors [20] introduce an algorithmic
framework using the Bregman projection to address the split feasibility
problem as follows:

$\operatorname*{minimize\quad}~{}\omega(x),~{}\mbox{subject
to}~{}x\in\bigcap_{i=1}^{N}C_{i},$ (1.3)

where $C_{i}$ are closed convex sets that are given by imposing convex
constraints $Q_{i}\subset\mathbb{R}^{N_{i}}$ in the range of a matrix
$A_{i}\in\mathbb{R}^{N_{i}\times n}$, i.e.,
$C_{i}=\\{x\in\mathbb{R}^{n}:A_{i}x\in Q_{i}\\}$. Within this framework, the
system of linear equations $Ax=b$ can be reformulated as a split feasibility
problem, placing the LBreI method as a special case within a broader class of
methods that includes popular techniques like the Kaczmarz method [15] and the
Landweber method [13]. Moreover, (1.3) also covers more generalized problems,
including other objective functions $\omega$ and incorporation of non-Gaussian
noise models, and box constraints.

Despite these advancements, it is important to note the primary theoretical
limitation of this framework[20]: the convergence analysis remains incomplete.
In particular, they only give a monotonic decrease in terms of Bregman
distance. The convergence to optimal solution and convergence rate are
unknown. This naturally gives rise to the following question:

> Can we design a unified algorithmic framework that covers LBreI with a
> completed convergence result?

We affirmatively answer this question by offering a novel perspective on LBreI
through the lens of bilevel optimization. We formulate the above problems as
the bilevel problem and propose a cut-and-project method to encompass LBreI
and the algorithm outlined in [20], providing a comprehensive convergence
result. In particular, we employ the cutting plane and Bregman projection
methods to create a straightforward yet broadly applicable algorithmic
framework, accompanied by a comprehensive convergence analysis. It is
noteworthy that our approach not only provides a clear geometry intuition but
also demonstrates significant potency in surpassing existing theoretical
frameworks.

Besides the previously mentioned several understandings of LBreI, there are
other important extensions. For example, the authors of [8] proposed the
singular value thresholding (SVT) algorithm which is a matrix-form variant of
LBreI and has become a classic method for matrix completion. Additionally, the
authors of [5] extended the LBreI method to a larger class of non-convex
functions for a wider range of imaging applications. In this study, we will
rediscover the main result in [10] as a special case within our proposed
framework, utilizing an essentially different line of thought.

The main contributions of this paper are summarized as follows:

  * •

By introducing a bilevel optimization formulation, we present a cut-and-
project perspective to reexamine Bregman regularized iterations. This method
not only encompasses existing Bregman regularized iterations but also applies
to solving more general inverse problems, including certain sparse noise
models.

  * •

A detailed analysis is conducted on the impact of step-sizes on the algorithm,
particularly identifying the precise range of step-sizes that can guide the
selection of practical step-sizes. Furthermore, we provide a convergence
guarantee to a feasible point with a large range of step-sizes, which are
wider compared to existing literature.

  * •

To establish the convergence of the proposed algorithmic framework to the
optimal solution of the problem, we provide a unified convergence condition.
When applied to specific examples, our condition is less restrictive compared
to conditions proposed in [30].

  * •

A Bregman distance growth condition is introduced to ensure the linear
convergence of the algorithm. Additionally, we derive sublinear convergence
results for the algorithm, which has only been analyzed for specific problems
so far, such as the SVT algorithm [8] or cases where the regularization
function is smooth.

The outline of the paper is as follows. Section 2 presents some preliminaries
and discussions. In Section 3, we give the bilevel optimization formulation,
present the proposed algorithm framework and give the descent lemma. The main
convergence results, including the convergence to a feasible point, the
convergence to the optimal solution, the linear convergence and sublinear
convergence, are provided in Section 4. The numerical experiments of our
algorithm are given in Section 5.

##  2 Preliminaries

###  2.1 Notation

In this paper, we restrict our analysis to real finite dimensional spaces
$\mathbb{R}^{n}$. We use $\langle\cdot,\cdot\rangle$ to denote the inner
product and $\|\cdot\|$ to denote the Euclidean norm. For a multi-variables
function $f(x,y)$, we use $\nabla_{x}f$ (respectively, $\nabla_{y}f$) to
denote the gradient of $f$ with respective to $x$ (respectively, $y$). Let
$Q\subset\mathbb{R}^{n}$ be a nonempty set; the distance function onto the set
$Q$ is defined by ${\mathbf{dist}}(x,Q):=\inf_{y\in Q}\|x-y\|$. The indicator
function of a set $\mathcal{C}$, denoted by $\delta_{\mathcal{C}}$, is set to
be zero on $\mathcal{C}$ and $+\infty$ otherwise.

###  2.2 Convex analysis tools

Some basic notations and facts about convex analysis will be used in our
results.

######  Definition 2.1.

A function $\omega:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is convex if for any
$\alpha\in[0,1]$ and $u,v\in\mathbb{R}^{d}$, we have

$\omega(\alpha u+(1-\alpha)v)\leq\alpha\omega(u)+(1-\alpha)\omega(v);$

and strongly convex with modulus $\mu>0$ if for any $\alpha\in[0,1]$ and
$u,v\in\mathbb{R}^{d}$, we have

$\omega(\alpha
u+(1-\alpha)v)\leq\alpha\omega(u)+(1-\alpha)\omega(v)-\frac{1}{2}\mu\alpha(1-\alpha)\|u-v\|^{2}.$

Furthermore, $\omega$ is strictly convex if for any $u\neq v$ and
$\alpha\in(0,1)$ we have

$\omega(\alpha u+(1-\alpha)v)<\alpha\omega(u)+(1-\alpha)\omega(v).$

######  Definition 2.2.

Let $\omega:\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a convex function. The
subdifferential of $\omega$ at $u\in\mathbb{R}^{d}$ is defined as

$\partial\omega(u):=\\{u^{*}\in\mathbb{R}^{d}:\omega(v)\geq\omega(u)+\langle
u^{*},v-u\rangle,\quad\forall v\in\mathbb{R}^{d}\\}.$

The elements of $\partial\omega(u)$ are called the subgradients of $\omega$ at
$u$.

The subdifferential generalizes the classical concept of differential because
of the well-known fact that $\partial\omega(u)=\\{\nabla\omega(u)\\}$ when the
function $\omega$ is differentiable. In terms of the subdifferential, the
strong convexity in Definition 2.1 can be equivalently stated as [14]: For any
$u,v\in\mathbb{R}^{d}$ and $v^{*}\in\partial\omega(v)$, we have

$\omega(u)\geq\omega(v)+\langle v^{*},u-v\rangle+\frac{\mu}{2}\|u-v\|^{2}.$
(2.1)

######  Definition 2.3.

Let $\omega:\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a convex function. The
conjugate of $\omega$ is defined as

$\omega^{*}(u^{*})=\sup_{v\in\mathbb{R}^{d}}\\{\langle
u^{*},v\rangle-\omega(v)\\}.$

######  Definition 2.4.

A function $\omega:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is gradient-Lipschitz-
continuous with modulus $L>0$ if it is continuously differentiable over
$\mathbb{R}^{n}$ and its gradient is Lipschitz continuous in the sense that
for any $u,v\in\mathbb{R}^{d}$, we have

$\|\nabla\omega(u)-\nabla\omega(v)\|\leq L\|u-v\|.$

The following facts are well-known, and can be found in the classic textbooks
[24] and [14].

######  Lemma 2.1.

Let $\omega:\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a strongly convex function
with modulus $\mu>0$. Then we have that

  * •

its conjugate $\omega^{*}$ is gradient-Lipschitz-continuous with modulus
$\frac{1}{\mu}$;

  * •

the conditions $\omega(u)+\omega^{*}(u^{*})=\langle u,u^{*}\rangle$,
$u^{*}\in\partial\omega(u)$, and $u\in\partial\omega^{*}(u^{*})$ are
equivalent.

###  2.3 Bregman distance tools

The Bregman distance, originally introduced in [6], is a very powerful concept
in many fields where distances are involved. Recently, many variants of
Bregman distances have appeared, as seen in [2, 16, 23]. For simplicity as
well as generality, we choose the Bregman distance defined by a strongly
convex function.

######  Definition 2.5.

Let $\omega:\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a strongly convex function
with modulus $\mu>0$. The Bregman distance $D_{\omega}^{v^{*}}(u,v)$ between
$u,v\in\mathbb{R}^{d}$ with respect to $\omega$ and a subgradient
$v^{*}\in\partial\omega(v)$ is defined by

$D_{\omega}^{v^{*}}(u,v):=\omega(u)-\omega(v)-\langle v^{*},u-v\rangle.$ (2.2)

In the following, we state three basic facts about the Bregman distance, which
will be used later in our analysis. It should be pointed out that the results
are well-known (see e.g. [16, 17]). We list them here, along with a brief
proof, for completeness.

######  Lemma 2.2.

Let $\omega:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is strongly convex with
modulus $\mu>0$. For any $u,p,q\in\mathbb{R}^{d}$ and
$p^{*}\in\partial\omega(p),q^{*}\in\partial\omega(q)$, we have that

$D_{\omega}^{p^{*}}(u,p)-D_{\omega}^{q^{*}}(u,q)+D_{\omega}^{q^{*}}(p,q)=\langle
q^{*}-p^{*},u-p\rangle,$ (2.3)
$D_{\omega}^{q^{*}}(p,q)=D_{{\omega}^{*}}^{p}(q^{*},p^{*}),$ (2.4)

and

$D_{\omega}^{q^{*}}(p,q)\geq\frac{\mu}{2}\|p-q\|^{2}.$ (2.5)

###### Proof.

The equality (2.3) follows from (2.2), and the inequality (2.5) from (2.1) and
(2.2). In order to obtain (2.4), we derive that

$\displaystyle\begin{array}[]{lll}D_{\omega}^{q^{*}}(p,q)&=&\omega(p)-\omega(q)-\langle
q^{*},p-q\rangle\\\ &=&\langle p,p^{*}\rangle-\omega^{*}(p^{*})-\langle
q,q^{*}\rangle+\omega^{*}(q^{*})-\langle q^{*},p-q\rangle\\\
&=&\omega^{*}(q^{*})-\omega^{*}(p^{*})-\langle p,q^{*}-p^{*}\rangle\\\
&=&D_{{\omega}^{*}}^{p}(q^{*},p^{*}),\end{array}$ (2.10)

where the second and fourth lines follow by using the second part of Lemma 2.1
and the condition $p^{*}\in\partial\omega(p),q^{*}\in\partial\omega(q)$. ∎

Now, we introduce the concept of Bregman projection.

######  Definition 2.6 (Bregman projection [17]).

Let $\omega:\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a strongly convex function
with modulus $\mu>0$, $C\subset\mathbb{R}^{n}$ be a nonempty closed convex
set. Given $x\in\mathbb{R}^{n}$ and $x^{*}\in\partial\omega(x)$, the Bregman
projection of $x$ onto $C$ with respect to $\omega$ and $x^{*}$ is the point
$\Pi_{C}^{x^{*}}(x)\in C$ such that

$D^{x^{*}}_{\omega}(\Pi_{C}^{x^{*}}(x),x)=\min_{y\in
C}D^{x^{*}}_{\omega}(y,x).$ (2.11)

Note that the point $\Pi_{C}^{x^{*}}(x)$ exists and is unique due to the
strong convexity of the objective function in (2.11); besides, the notation
$\Pi_{C}^{x^{*}}(x)$ is dependent on the function $\omega$ and reduces to the
orthogonal projection, denoted by $\mathcal{P}_{C}(x)$, when
$\omega(x)=\frac{1}{2}\|x\|^{2}.$ Similar to the traditional orthogonal
projection theorem, the Bregman projection can be characterized by a
variational inequality.

######  Lemma 2.3 (Generalized projection theorem[17]).

Let $\omega:\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a strongly convex function
with modulus $\mu>0$ and $C\subset\mathbb{R}^{n}$ be a nonempty closed convex.
Then a point $z\in C$ is the Bregman projection of $x$ onto $C$ with respect
to $\omega$ and $x^{*}\in\partial\omega(x)$ if and only if there is some point
$z^{*}\in\partial\omega(z)$ such that one of the following equivalent
conditions are fulfilled:

$\langle z^{*}-x^{*},y-z\rangle\geq 0,~{}~{}\forall~{}~{}y\in C,$ (2.12)
$D^{z^{*}}_{\omega}(y,z)\leq
D^{x^{*}}_{\omega}(y,x)-D^{x^{*}}_{\omega}(z,x),~{}~{}\forall~{}~{}y\in C.$
(2.13)

Any such $z^{*}$ is called an admissible subgradient for
$z=\Pi_{C}^{x^{*}}(x).$ In particular, when $\omega(x)=\frac{1}{2}\|x\|^{2}$,
the condition $x^{*}\in\partial\omega(x)$ reduces to $x^{*}=x$ and hence the
inequality (2.12) just becomes the orthogonal projection characterization.

##  3 The proposed algorithmic framework

###  3.1 The problem formulation

We are interested in the following bilevel optimization problem which consists
of inner and outer levels. The inner level is the classical optimization
problem

$\operatorname*{minimize\quad}_{x\in\mathbb{R}^{n}}f(x),$ (3.1)

where the optimal set is denoted by $\overline{X}$, assumed to be nonempty.
The out level seeks to find the optimal solution of problem (3.1) which
minimizes a given strongly convex function $\omega$:

$\operatorname*{minimize\quad}_{x\in\overline{X}}\omega(x).$ (3.2)

The inner level objective function $f$ and the outer level objective function
$\omega$ are assumed to satisfy the following:

  * •

$\omega$ is strongly convex over $\mathbb{R}^{n}$ with parameter $\mu>0$;

  * •

$f$ is convex and gradient-Lipschitz-continuous with modulus $L>0$.

By the strong convexity of $\omega$, problem (3.2) has a unique solution,
denoted by $\bar{x}$. Putting the ingredients above together, we obtain the
following bilevel optimization:

$\displaystyle\operatorname*{minimize\quad}~{}~{}\omega(x)~{}$ (3.3)
$\displaystyle\mbox{subject to}~{}~{}x\in\arg\min_{x}f(x).$

The terminology of bilevel optimization, mainly following from papers [4], is
simpler than classical bilevel optimization where the solution to the inner
level optimization usually depends on variables or parameters.

This bilevel optimization encompasses various specific problems. When
$\omega(x)=\lambda\|x\|_{1}+\frac{1}{2}\|x\|^{2}$ and
$f(x)=\frac{1}{2}\|Ax-b\|^{2}$, (3.3) simplifies to the least squares problem
(1.2). Furthermore, if $b\in\mathcal{R}(A)$, the inner-level optimization
reduces to $Ax=b$, and thus (3.3) covers (1.1). Moreover, (3.3) encompasses
the split feasibility problem (1.3) by defining
$f(x):=\sum_{i=1}^{N}{\mathbf{dist}}^{2}(A_{i}x,Q_{i})$, which is convex and
gradient-Lipschitz-continuous. This is due to that minimizing $f$ is
equivalent to finding $x\in\bigcap_{i=1}^{N}C_{i}$.

###  3.2 Sketch of idea

Our proposed iterative scheme will be grounded on a straightforward and
intuitive geometric understanding. Let us start with a given point $x$ which
does not lie in the optimal set $\overline{X}$. We can then construct a
hyperplane to separate this given point $x$ from the optimal set
$\overline{X}$. In order to gradually approach the optimal set $\overline{X}$,
we project the given point $x$ onto the hyperplane, thereby obtaining updated
points that progressively draw nearer to the optimal set. However, our aim is
not solely for the iterates to converge to any solution of the inner-level
optimization but specifically to the unique solution $\bar{x}$ of the bilevel
optimization problem. To achieve this, we select the objective function
$\omega$ of the outer-level optimization as a kernel to generate the Bregman
distance $D_{\omega}$. Subsequently, we utilize the Bregman projection to
replace the orthogonal projection.

From the discussion above, the main idea of our method consists of two steps:
constructing a halfspace (with a hyperplane boundary) and performing the
Bregman projection onto the halfspace. The first step relies on the Baillon-
Haddad theorem [1], which states that the convexity and Lipschitz continuity
of gradient of $f$ is equivalent to the following inequality

$\langle\nabla f(u)-\nabla f(v),u-v\rangle\geq\frac{1}{L}\|\nabla f(u)-\nabla
f(v)\|^{2}.$

Note that for any $x\in\overline{X}$, $\nabla f(x)=0$. Then, for any given
point $z$, we have the following inclusion

$\overline{X}\subset H_{z}:=\\{x\in\mathbb{R}^{n}:\langle\nabla
f(z),z-x\rangle\geq\frac{1}{L}\|\nabla f(z)\|^{2}\\},$

where $H_{z}$ is a halfspace which separates the given point $z$ and the
optimal set $\overline{X}$. In particular, if $\nabla f(z)\neq 0$, then the
separation is strict in the sense that $z$ is not in $H_{z}$ since in that
case $z$ fails to satisfy the inequality $\langle\nabla
f(z),z-x\rangle\geq\frac{1}{L}\|\nabla f(z)\|^{2}$. The next step is to
project $z$ onto the halfspace $H_{z}$. Denote the project operator onto the
halfspace $H_{z}$ as $\mathcal{P}$; then for any point
$\hat{x}\in\overline{X}$, we can derive that

$\|\mathcal{P}(z)-\hat{x}\|=\|\mathcal{P}(z)-\mathcal{P}(\hat{x})\|\leq\|z-\hat{x}\|,$

where the first relationship is due to $\hat{x}\in\overline{X}\subset H_{z}$
and the inequality follows from the nonexpansiveness of the project operator.
Therefore, it is possible to get closer to the optimal set $\overline{X}$
after the project operation.

The idea of our method can be simply summarized as follows:

s1.

Construct cutting halfspaces: $H_{k}=\\{x\in\mathbb{R}^{n}:\langle\nabla
f(x_{k}),x_{k}-x\rangle\geq\frac{1}{L}\|\nabla f(x_{k})\|^{2}\\}.$

s2.

Update via Bregman projections: $x_{k+1}:=\Pi_{H_{k}}^{x^{*}_{k}}(x_{k})$.

Lastly, we would like to mention two closely related works [20] and [3], which
have greatly inspired us. In [20], the authors presented a generalized
algorithmic framework based on Bregman projections, which includes LBreI as a
special case. However, their framework is restricted solely to solving split
feasibility problems (1.3). On the other hand, in [3], the authors introduced
a novel algorithm based on the cutting plane method to address a class of
bilevel convex problems. Nevertheless, their method necessitates the outer-
level function to be smooth, and at each iteration, they are required to solve
a constrained optimization problem, which may pose computational challenges.
Our approach is motivated by amalgamating the strengths of these two works,
aiming to overcome their respective limitations.

###  3.3 The Bregman projection as a primal-dual problem

The key computation in our method involves performing Bregman projections. In
this part, we will formulate it as a primal-dual problem. To this end, we
assume that $\nabla f(x_{k})\neq 0$ and let $a_{k}:=\nabla f(x_{k})$ and

$\beta_{k}:=\langle\nabla f(x_{k}),x_{k}\rangle-\frac{1}{L}\|\nabla
f(x_{k})\|^{2}.$

Then the halfspaces defined in the first step of our method can be simplified
as follows

$H_{k}=\\{x\in\mathbb{R}^{n}:\langle a_{k},x\rangle\leq\beta_{k}\\}.$ (3.4)

From Definition 2.6, the Bregman projections have the following formulation

$\Pi_{H_{k}}^{x^{*}_{k}}(x_{k})=\arg\min_{x\in
H_{k}}D_{\omega}^{x_{k}^{*}}(x,x_{k}),$

which is equivalent to solving the following constrained optimization

$\displaystyle\operatorname*{minimize\quad}$ $\displaystyle\omega(x)-\langle
x_{k}^{*},x\rangle$ subject to $\displaystyle\langle
a_{k},x\rangle\leq\beta_{k}.$

In order to rewrite (3.3) as an unconstrained optimization, we let
$\Omega_{k}=\\{s\in\mathbb{R}:s\leq\beta_{k}\\}$,
$\varphi(s)=\delta_{\Omega_{k}}(s)$, and $\psi(x)=\omega(x)-\langle
x_{k}^{*},x\rangle$; then (3.3) can be equivalently rewritten as

$\operatorname*{minimize\quad}_{x\in\mathbb{R}^{n}}\psi(x)+\varphi(a_{k}^{T}x),$
(3.6)

whose Fenchel-Rockafellar dual problem is

$\operatorname*{minimize\quad}_{t\in\mathbb{R}}\psi^{*}(-ta_{k})+\varphi^{*}(t).$
(3.7)

Note that $\psi^{*}(x)=\omega^{*}(x_{k}^{*}+x)$ and

$\displaystyle\varphi^{*}(t)=\delta^{*}_{\Omega_{k}}(t):=\left\\{\begin{array}[]{lll}\beta_{k}t,&\textrm{if}~{}~{}t\geq
0,\\\ +\infty,&\textrm{otherwise}.\end{array}\right.$

The dual problem (3.7) can be further written as

$\operatorname*{minimize\quad}_{t\geq
0}\omega^{*}(x_{k}^{*}-ta_{k})+\beta_{k}t.$ (3.9)

Let $x\in\mathbb{R}^{n},t\in\mathbb{R}$. Then $x$ is a primal solution and $t$
is a dual solution if and only if the KKT condition holds:

$\displaystyle\left\\{\begin{array}[]{rll}-ta_{k}&\in&\partial\psi(x)=\partial\omega(x)-x_{k}^{*},\\\
t&\in&\partial\varphi(a_{k}^{T}x)=N_{\Omega_{k}}(a_{k}^{\top}x),\end{array}\right.$
(3.12)

where the normal cone

$\displaystyle
N_{\Omega_{k}}(a_{k}^{\top}x)=\left\\{\begin{array}[]{ll}\mathbb{R}_{+},&a_{k}^{T}x=\beta_{k},\\\
\\{0\\},&a_{k}^{T}x<\beta_{k}.\end{array}\right.$

We observe that the case of $t\in\\{0\\}$ and $a_{k}^{T}x<\beta_{k}$ can be
excluded. Indeed, let us focus on the dual objective function

$g(t)=\omega^{*}(x_{k}^{*}-ta_{k})+\beta_{k}t.$

Its derivative can be calculated as
$g^{\prime}(t)=-a_{k}^{T}\nabla\omega^{*}(x_{k}^{*}-ta_{k})+\beta_{k}$ and
hence we can derive that

$\displaystyle\begin{array}[]{lll}g^{\prime}(0)&=&-a_{k}^{T}\nabla\omega^{*}(x_{k}^{*})+\beta_{k}\\\
&=&-a_{k}^{T}x_{k}+\beta_{k}\\\ &=&-\frac{1}{L}\|\nabla
f(x_{k})\|^{2}<0,\end{array}$

where we have utilized the fact of $x_{k}^{*}\in\partial\omega(x_{k})$, the
expressions of $a_{k}$ and $\beta_{k}$, and the assumption $\nabla
f(x_{k})\neq 0$. Hence, the function $g(t)$ must be decreasing for small
enough $t>0$. In other words, the dual solution $t$ must be obtained as a
positive number. Therefore, the dual problem can be equivalently rewritten as

$\operatorname*{minimize\quad}_{t>0}\omega^{*}(x_{k}^{*}-ta_{k})+\beta_{k}t,$
(3.15)

and the KKT condition can be reformulated as the following form after some
simple transformations:

$\displaystyle\left\\{\begin{array}[]{rll}x&=&\nabla\omega^{*}(x_{k}^{*}-ta_{k}),\\\
a_{k}^{T}x&=&\beta_{k},\\\ t&>&0.\end{array}\right.$ (3.19)

In particular, the condition $a_{k}^{T}x=\beta_{k}$ implies that the Bregman
projection must lie in the hyperplane $\\{x\in\mathbb{R}^{n}:\langle\nabla
f(x_{k}),x_{k}-x\rangle=\frac{1}{L}\|\nabla f(x_{k})\|^{2}\\}.$ Based on the
KKT condition (3.19), we can deduce the following result.

######  Theorem 3.1.

Let $\hat{t}_{k}$ be a dual solution to (3.15). Then the Bregman projection of
$x_{k}$ onto the halfspace $H_{k}$ can be obtained by the following
formulation

$\Pi_{H_{k}}^{x^{*}_{k}}(x_{k})=\nabla\omega^{*}(x_{k}^{*}-\hat{t}_{k}a_{k})$
(3.20)

and moreover $\hat{t}_{k}$ belongs to $(0,\frac{\mu}{L}]$. On the other hand,
if the point $\nabla\omega^{*}(x_{k}^{*}-a_{k}\tilde{t}_{k})$ with
$\tilde{t}_{k}\in(0,\frac{\mu}{L}]$ lies in the hyperplane
$\\{x\in\mathbb{R}^{n}:\langle\nabla
f(x_{k}),x_{k}-x\rangle=\frac{1}{L}\|\nabla f(x_{k})\|^{2}\\}$, then it must
be the unique point of the Bregman projection, that is

$\nabla\omega^{*}(x_{k}^{*}-a_{k}\tilde{t}_{k})=\nabla\omega^{*}(x_{k}^{*}-\hat{t}_{k}a_{k})=\Pi_{H_{k}}^{x^{*}_{k}}(x_{k}).$
(3.21)

###### Proof.

Since $\hat{t}_{k}$ is an exact solution to (3.15), by the first equation of
the KKT condition (3.19) we deduce that
$\nabla\omega^{*}(x_{k}^{*}-\hat{t}_{k}a_{k})$ must be the unique primal
solution to (3.6) and hence the desired result (3.20) holds. In terms of the
second equation of (3.19), we have

$a_{k}^{T}\nabla\omega^{*}(x_{k}^{*}-\hat{t}_{k}a_{k})=\beta_{k}.$

Using the expressions of $a_{k}$ and $\beta_{k}$ and the fact that
$x_{k}=\nabla\omega^{*}(x_{k}^{*})$, we further get

$\langle\nabla
f(x_{k}),\nabla\omega^{*}(x_{k}^{*})-\nabla\omega^{*}(x_{k}^{*}-\hat{t}_{k}\cdot\nabla
f(x_{k}))\rangle=\frac{1}{L}\|\nabla f(x_{k})\|^{2}.$ (3.22)

By the Lipschitz continuity of $\nabla\omega^{*}$, we can bound the left-hand
side of (3.22) as follows:

$\langle\hat{t}_{k}\cdot\nabla
f(x_{k}),\nabla\omega^{*}(x_{k}^{*})-\nabla\omega^{*}(x_{k}^{*}-\hat{t}_{k}\cdot\nabla
f(x_{k}))\rangle\geq\frac{1}{\mu}\|\hat{t}_{k}\cdot\nabla f(x_{k})\|^{2}.$
(3.23)

Combining (3.22) and (3.23), we get

$\frac{\hat{t}_{k}}{L}\|\nabla
f(x_{k})\|^{2}\geq\frac{1}{\mu}\|\hat{t}_{k}\cdot\nabla f(x_{k})\|^{2},$

which leads to $0<\hat{t}_{k}\leq\frac{\mu}{L}$. Therefore, the proof of the
first part is completed.

In order to show the second part, we introduce two auxiliary points:

$\displaystyle\tilde{x}_{k+1}$ $\displaystyle=$
$\displaystyle\nabla\omega^{*}(x_{k}^{*}-a_{k}\tilde{t}_{k}),$
$\displaystyle\tilde{x}_{k+1}^{*}$ $\displaystyle=$ $\displaystyle
x_{k}^{*}-a_{k}\tilde{t}_{k}.$

Then, $\tilde{x}_{k+1}^{*}\in\partial\omega(\tilde{x}_{k+1})$ and the point
$\tilde{x}_{k+1}$ belongs to $H_{k}$ due to the assumption that
$\tilde{x}_{k+1}$ lies in the hyperplane $\\{x\in\mathbb{R}^{n}:\langle\nabla
f(x_{k}),x_{k}-x\rangle=\frac{1}{L}\|\nabla f(x_{k})\|^{2}\\}$. By the
uniqueness of the Bregman projection, it suffices to show that
$\tilde{x}_{k+1}$ is the Bregman projection of $x_{k}$ onto the halfspace
$H_{k}$. Based on Lemma 2.3, we only need to verify the following
relationship:

$\langle\tilde{x}_{k+1}^{*}-x_{k}^{*},x-\tilde{x}_{k+1}\rangle\geq
0,~{}~{}\forall x\in H_{k}.$

Actually, we can derive that for any $x\in H_{k}$,

$\displaystyle\langle\tilde{x}_{k+1}^{*}-x_{k}^{*},x-\tilde{x}_{k+1}\rangle$
$\displaystyle=$ $\displaystyle\langle-\tilde{t}_{k}\nabla
f(x_{k}),x-\nabla\omega^{*}(x_{k}^{*}-a_{k}\tilde{t}_{k})\rangle$
$\displaystyle=$ $\displaystyle\tilde{t}_{k}(\langle\nabla
f(x_{k}),\nabla\omega^{*}(x_{k}^{*}-a_{k}\tilde{t}_{k})\rangle-\langle\nabla
f(x_{k}),x\rangle)$ $\displaystyle\geq$
$\displaystyle\tilde{t}_{k}(\langle\nabla
f(x_{k}),\nabla\omega^{*}(x_{k}^{*}-a_{k}\tilde{t}_{k})\rangle-\beta_{k})$
$\displaystyle=$ $\displaystyle-\tilde{t}_{k}\left(\langle\nabla
f(x_{k}),\nabla\omega^{*}(x_{k}^{*})-\nabla\omega^{*}(x_{k}^{*}-a_{k}\tilde{t}_{k})\rangle-\frac{1}{L}\|\nabla
f(x_{k})\|^{2}\right)$ $\displaystyle\geq$
$\displaystyle-\tilde{t}_{k}\left(\|\nabla
f(x_{k})\|\cdot\|\nabla\omega^{*}(x_{k}^{*})-\nabla\omega^{*}(x_{k}^{*}-a_{k}\tilde{t}_{k})\|-\frac{1}{L}\|\nabla
f(x_{k})\|^{2}\right)$ $\displaystyle\geq$
$\displaystyle-\tilde{t}_{k}\left(\frac{\tilde{t}_{k}}{\mu}\|\nabla
f(x_{k})\|^{2}-\frac{1}{L}\|\nabla f(x_{k})\|^{2}\right)\geq 0,$

where the first inequality follows from the condition $x\in H_{k}$, the second
from the Cauchy-Schwartz inequality, the third from the Lipschitz continuity
of $\nabla\omega^{*}$, and the last one is due to the condition
$\tilde{t}_{k}\in(0,\frac{\mu}{L}]$. Thus, from the KKT condition (3.19), we
can see that $\tilde{t}_{k}$ must be a solution to the dual problem. Combining
with the deduced result that any exact solutions must belong to
$(0,\frac{\mu}{L}]$, we can further conclude that $(0,\frac{\mu}{L}]$ is just
the set of dual solutions. This completes the proof. ∎

###  3.4 Iterate scheme and descent lemma

Based on the previous study, we are ready to present our method concretely.
The detailed algorithm is provided in Algorithm 1, where the main iterative
scheme consists of two steps: updating the primal variable $x_{k+1}$ and the
dual variable $x_{k+1}^{*}$. Here, we do not specify the step-size $t_{k}$,
which can be obtained through various choices. As shown in Theorem 3.1, when
the step-size $t_{k}$ is chosen as an exact solution of the dual problem
(3.15), the update rule of $x_{k+1}$ in Algorithm 1 is strictly equivalent to
the Bregman projection onto $H_{k}$, i.e.,
$x_{k+1}=\Pi_{H_{k}}^{x^{*}_{k}}(x_{k})$. It should be emphasized that
achieving precise solutions $t_{k}$ necessitates a requirement for computing
the Lipschtz constant $L$ of function $f$, which may be relatively costly in
certain situations. Nevertheless, we will demonstrate that convergence results
can be established across a range of step-size $t_{k}$. Moreover, we also give
the numerical comparison of our algorithm with different step-size $t_{k}$ in
the experiment.

Algorithm 1 Bregman projection method based on cutting halfspaces

1: Initialize:  choose $x_{0}\in\mathbb{R}^{n}$ and
$x_{0}^{*}\in\partial\omega(x_{0})$.

2: for $k=0,\cdots,$ do

3: Update the step-size $t_{k}$ via a certain rule.

4: Update the primal variable $x_{k+1}$:
$x_{k+1}=\nabla\omega^{*}(x_{k}^{*}-t_{k}\cdot\nabla f(x_{k})).$ (3.24)

5: Update the dual variable $x_{k+1}^{*}$:
$x_{k+1}^{*}=x_{k}^{*}-t_{k}\cdot\nabla f(x_{k}).$ (3.25)

6: end for

Algorithm 1 encompasses a broad range of existing algorithms through
appropriate choices of functions $f$ and $\omega$. Specifically, when
$f(x)=\frac{1}{2}\|Ax-b\|^{2}$ and
$\omega(x)=\lambda\|x\|_{1}+\frac{1}{2}\|x\|^{2}$, it reduces to the
linearized Bregman iteration method, as evidenced by Corollary 4.1. Moreover,
setting $f(x)=\frac{1}{2}{\mathbf{dist}}^{2}(Ax,Q)$, where $Q$ is a given
closed convex set in $\mathbb{R}^{m}$, yields the iterative method proposed in
[20] (referenced as (4.13)). This illustrates the versatility of Algorithm 1
in encompassing various optimization methods tailored to specific problem
structures.

The following lemma shows that Algorithm 1 achieves a monotonic decrease in
terms of Bregman distance with different choices of the step-size $t_{k}$.

######  Lemma 3.1 (Basic descent lemma).

Let $(x_{k},x_{k}^{*})$ be the $k$-iterate generated in Algorithm 1. If the
step-size $t_{k}$ is an exact solution to the dual problem (3.15) or equals to
$\frac{\mu}{L}$ exactly, then we have

$D_{\omega}^{x_{k+1}^{*}}(x,x_{k+1})\leq
D_{\omega}^{x_{k}^{*}}(x,x_{k})-\frac{\mu}{2L^{2}}\|\nabla
f(x_{k})\|^{2},~{}\forall x\in H_{k}.$ (3.26)

If the step-size $t_{k}$ lies in the interval $(0,\frac{2\mu}{L}]$, then we
have

$D_{\omega}^{x_{k+1}^{*}}(x,x_{k+1})\leq
D_{\omega}^{x_{k}^{*}}(x,x_{k})-\left(\frac{t_{k}}{L}-\frac{t_{k}^{2}}{2\mu}\right)\|\nabla
f(x_{k})\|^{2},~{}\forall x\in H_{k}.$ (3.27)

###### Proof.

First of all, let us show (3.26), where the iterate sequences $\\{x_{k}\\}$
and $\\{x_{k}^{*}\\}$ are generated by Algorithm 1 with $t_{k}=\hat{t}_{k}$
being the exact solution to the dual problem (3.15). From Lemma 2.1 and the
definition of Bregman distance, we can express the Bregman distance in the
following form

$D_{\omega}^{x^{*}}(y,x)=\omega^{*}(x^{*})-\langle
x^{*},y\rangle+\omega(y),~{}\forall x^{*}\in\partial\omega(x).$

Using this expression and the updated formulation (3.25) of $x_{k+1}^{*}$, we
derive that

$\displaystyle D_{\omega}^{x_{k+1}^{*}}(x,x_{k+1})$ $\displaystyle=$
$\displaystyle\omega^{*}(x_{k+1}^{*})-\langle x_{k+1}^{*},x\rangle+\omega(x)$
$\displaystyle=$ $\displaystyle\omega^{*}(x_{k}^{*}-t_{k}\cdot\nabla
f(x_{k}))-\langle x_{k}^{*},x\rangle+\langle t_{k}\cdot\nabla
f(x_{k}),x\rangle+\omega(x).$

By the definition of $H_{k}$, we know that if $x\in H_{k}$, then
$\langle\nabla f(x_{k}),x\rangle\leq\beta_{k}.$ Thus, combining with the
Lipschitz continuity of $\nabla\omega^{*}$ and the fact that $t_{k}$ is the
minimizer to the dual objective function
$\omega^{*}(x_{k}^{*}-ta_{k})+\beta_{k}t$, we can further derive that for any
$x\in H_{k}$ and $t>0$,

$\displaystyle D_{\omega}^{x_{k+1}^{*}}(x,x_{k+1})$ $\displaystyle\leq$
$\displaystyle\omega^{*}(x_{k}^{*}-t_{k}\cdot\nabla
f(x_{k}))+t_{k}\cdot\beta_{k}-\langle x_{k}^{*},x\rangle+\omega(x)$
$\displaystyle\leq$ $\displaystyle\omega^{*}(x_{k}^{*}-t\cdot\nabla
f(x_{k}))+t\cdot\beta_{k}-\langle x_{k}^{*},x\rangle+\omega(x)$
$\displaystyle\leq$
$\displaystyle\omega^{*}(x_{k}^{*})-\langle\nabla\omega^{*}(x_{k}^{*}),t\nabla
f(x_{k})\rangle+\frac{t^{2}}{2\mu}\|\nabla
f(x_{k})\|^{2}+t\cdot\beta_{k}-\langle x_{k}^{*},x\rangle+\omega(x)$
$\displaystyle=$ $\displaystyle\omega^{*}(x_{k}^{*})-t\langle x_{k},\nabla
f(x_{k})\rangle+\frac{t^{2}}{2\mu}\|\nabla
f(x_{k})\|^{2}+t\cdot\beta_{k}-\langle x_{k}^{*},x\rangle+\omega(x)$
$\displaystyle=$ $\displaystyle\omega^{*}(x_{k}^{*})-\langle
x_{k}^{*},x\rangle+\omega(x)+\left(\frac{t^{2}}{2\mu}-\frac{t}{L}\right)\|\nabla
f(x_{k})\|^{2}$ $\displaystyle=$ $\displaystyle
D_{\omega}^{x_{k}^{*}}(x,x_{k})+\left(\frac{t^{2}}{2\mu}-\frac{t}{L}\right)\|\nabla
f(x_{k})\|^{2},$

where the expression of $\beta_{k}$ and the fact that
$x_{k}=\nabla\omega^{*}(x_{k}^{*})$ have been used. Therefore, for any $x\in
H_{k}$ we have that

$\displaystyle D_{\omega}^{x_{k+1}^{*}}(x,x_{k+1})$ $\displaystyle\leq$
$\displaystyle
D_{\omega}^{x_{k}^{*}}(x,x_{k})+\inf_{t>0}\left(\frac{t^{2}}{2\mu}-\frac{t}{L}\right)\|\nabla
f(x_{k})\|^{2}$ $\displaystyle=$ $\displaystyle
D_{\omega}^{x_{k}^{*}}(x,x_{k})-\frac{\mu}{2L^{2}}\|\nabla f(x_{k})\|^{2}.$

Note that the descent property can also be obtained with
$t_{k}=\frac{\mu}{L}$. Hence, the proof of the first part has been done. The
second part follows by observing that when $t\in(0,\frac{2\mu}{L}]$, the term
$\left(\frac{t}{L}-\frac{t^{2}}{2\mu}\right)$ is nonnegative so that (3.27) is
always a descent property. Thus, we complete the whole proof. ∎

##  4 Convergence Analysis

###  4.1 Convergence to a feasible point

######  Lemma 4.1.

Let $(x_{k},x_{k}^{*})$ be the $k$-iterate generated in Algorithm 1. Assume
that the sequence $\\{x_{k}\\}$ converges to a point $x$. If for a given point
$y$ and a given subgradient sequence $\\{x_{k}^{*}\\}$ satisfying
$x_{k}^{*}\in\partial\omega(x_{k})$, the sequence
$\\{D_{\omega}^{x_{k}^{*}}(y,x_{k})\\}$ converges, then there must exist a
point $x^{*}$ in $\partial\omega(x)$ such that

$\lim_{k\rightarrow\infty}D_{\omega}^{x_{k}^{*}}(y,x_{k})=D_{\omega}^{x^{*}}(y,x).$
(4.1)

###### Proof.

First of all, the convergence of $\\{x_{k}\\}$ implies its boundedness, which
further implies the boundedness of the subgradient sequence $\\{x_{k}^{*}\\}$
by invoking Theorem [3, Theorem 3.16]. Hence, there exists a convergent
subsequence $\\{x_{k_{l}}^{*}\\}$ whose limit point is denoted by $x^{*}$.
Since $\\{x_{k}\\}$ converges to the point $x$, its any subsequence, including
the subsequence $\\{x_{k_{l}}\\}$, must converge to $x$ as well. The function
$\omega$ must be continuous by noting that it is a real-valued convex
function. Based on these facts, we derive that

$\displaystyle\lim_{k\rightarrow\infty}D_{\omega}^{x_{k}^{*}}(y,x_{k})$
$\displaystyle=$
$\displaystyle\lim_{l\rightarrow\infty}D_{\omega}^{x_{k_{l}}^{*}}(y,x_{k_{l}})$
$\displaystyle=$
$\displaystyle\lim_{l\rightarrow\infty}(\omega(y)-\omega(x_{k_{l}})-\langle
x_{k_{l}}^{*},y-x_{k_{l}}\rangle)$ $\displaystyle=$
$\displaystyle\omega(y)-\omega(x)-\langle x^{*},y-x\rangle.$

It remains to show that $x^{*}\in\partial\omega(x)$, or equivalently,
$x=\nabla\omega^{*}(x^{*})$. In fact, noting that
$x^{*}_{k_{l}}\in\partial\omega(x_{k_{l}})$, we derive that

$\displaystyle\|\nabla\omega^{*}(x^{*})-x\|$ $\displaystyle=$
$\displaystyle\lim_{l\rightarrow\infty}\|\nabla\omega^{*}(x^{*})-x_{k_{l}}\|$
$\displaystyle=$
$\displaystyle\lim_{l\rightarrow\infty}\|\nabla\omega^{*}(x^{*})-\nabla\omega^{*}(x_{k_{l}}^{*})\|$
$\displaystyle\leq$
$\displaystyle\lim_{l\rightarrow\infty}\frac{1}{\mu}\|x^{*}-x_{k_{l}}^{*}\|=0,$

where the inequality follows from the gradient-Lipschitz-continuity of
$\omega^{*}$ stated in Lemma 2.1. ∎

######  Theorem 4.1 (Convergence to an inner-level solution).

Let $\\{x_{k}\\}$ be generalized by Algorithm 1 with $t_{k}$ being an exact
solution to the dual problem (3.15) or a constant in $(0,\frac{2\mu}{L})$.
Then, there exists a point $x$ in $\overline{X}$ such that $x_{k}\rightarrow
x$ as $k\rightarrow+\infty$.

###### Proof.

By Lemma 3.1, for any fixed $x\in\overline{X}$ we have that
$\\{D_{\omega}^{x^{*}_{k}}(x,x_{k})\\}$ monotonically decreases and hence must
be a convergence sequence. Combining with the strong convexity of $\omega$
yields that

$\frac{\mu}{2}\|x-x_{k}\|^{2}\leq D_{\omega}^{x^{*}_{k}}(x,x_{k})\leq
D_{\omega}^{x^{*}_{0}}(x,x_{0}),$

which means that $\\{x_{k}\\}$ is bounded. It suffices to show that it has a
unique cluster point. To this end, let $\\{x_{k_{i}}\\}$ and $\\{x_{k_{j}}\\}$
be two subsequences of $\\{x_{k}\\}$ such that $x_{k_{i}}\rightarrow\hat{x}$
and $x_{k_{j}}\rightarrow\tilde{x}$; it reduces to show that
$\hat{x}=\tilde{x}$. In fact, invoking Lemma 4.1, we have that

$\lim_{i\rightarrow\infty}D_{\omega}^{x_{k_{i}}^{*}}(x,x_{k_{i}})=D_{\omega}^{\hat{x}^{*}}(x,\hat{x})$

and

$\lim_{j\rightarrow\infty}D_{\omega}^{x_{k_{j}}^{*}}(x,x_{k_{j}})=D_{\omega}^{\tilde{x}^{*}}(x,\tilde{x}).$

Due to the convergence of $\\{D_{\omega}^{x^{*}_{k}}(x,x_{k})\\}$, we
immediately get

$D_{\omega}^{\hat{x}^{*}}(x,\hat{x})=D_{\omega}^{\tilde{x}^{*}}(x,\tilde{x}).$
(4.2)

On the other hand, from Lemma 3.1 we can get

$\frac{\mu}{2L^{2}}\|\nabla f(x_{k})\|^{2}\leq
D_{\omega}^{x_{k}^{*}}(x,x_{k})-D_{\omega}^{x_{k+1}^{*}}(x,x_{k+1}).$

Using the convergence of $\\{D_{\omega}^{x^{*}_{k}}(x,x_{k})\\}$ and letting
$k$ tend to $\infty$, we get $\lim_{k\rightarrow+\infty}\|\nabla f(x_{k})\|=0$
and hence $\|\nabla f(\hat{x})\|=\lim_{i\rightarrow+\infty}\|\nabla
f(x_{k_{i}})\|=0$ and $\|\nabla
f(\tilde{x})\|=\lim_{j\rightarrow+\infty}\|\nabla f(x_{k_{j}})\|=0$. Thus,
both $\hat{x}$ and $\tilde{x}$ belong to $\overline{X}$. Now, we can take
$x=\hat{x}$ in (4.2) to get $D_{\omega}^{\tilde{x}^{*}}(\hat{x},\tilde{x})=0$
which implies that $\hat{x}=\tilde{x}$ by the strong convexity of $\omega$.
This completes the proof. ∎

###  4.2 Convergence to the optimal solution

######  Theorem 4.2 (Convergence to the optimal solution).

Assume that $f$ is in the form of $f(x)=g(Ax)$, where $A\in\mathbb{R}^{m\times
n}$ is a nonzero matrix and $g:\mathbb{R}^{m}\rightarrow\mathbb{R}$ is a
differentiable convex function satisfying one of the following conditions:

  1. (i)

$A$ is surjective and $g$ has a unique minimizer;

  2. (ii)

$g$ is a strictly convex function.

Let $\\{x_{k}\\}$ be generalized by Algorithm 1 with $t_{k}$ being an exact
solution to the dual problem (3.15) or a constant in $(0,\frac{2\mu}{L})$, and
with initial points $x_{0}$ and $x_{0}^{*}$ satisfying
$x_{0}^{*}\in\partial\omega(x_{0})\bigcap\mathcal{R}(A^{T}).$ Then, the
sequence $\\{x_{k}\\}$ converges to the optimal solution $\bar{x}$ of the
bilevel optimization problem (3.3).

###### Proof.

By Theorem 4.1, we know that there exists $\tilde{x}\in\overline{X}$ such that
$x_{k}\rightarrow\tilde{x}$ as $k\rightarrow\infty$. Hence, $\\{x_{k}\\}$ must
be a bounded sequence, and so is the subgradient sequence $\\{x_{k}^{*}\\}$ by
[3, Theorem 3.16] Moreover, based on the iterative scheme of $x_{k}^{*}$, we
deduce that

$x_{k+1}^{*}=x_{k}^{*}-t_{k}\nabla
f(x_{k})=\cdots=x_{0}^{*}-\sum_{i=0}^{k}t_{i}\nabla f(x_{i}).$

Since $x_{0}^{*}\in\partial\omega(x_{0})\bigcap\mathcal{R}(A^{T})$, there
exists a vector $y_{0}\in\mathbb{R}^{m}$ such that $x_{0}^{*}=A^{T}y_{0}$. Due
to the expression of $f(x)=g(Ax)$, we have $\nabla f(x)=A^{T}\nabla g(Ax)$.
Therefore,

$x_{k+1}^{*}=A^{T}y_{0}-\sum_{i=0}^{k}t_{i}A^{T}\nabla
g(Ax_{i})=A^{T}(y_{0}-\sum_{i=0}^{k}t_{i}\nabla g(Ax_{i})).$

Let $z_{k}=y_{0}-\sum_{i=0}^{k}t_{i}\nabla g(Ax_{i})$ and decompose it into
the direct sum of $z_{k}=z_{k}^{1}+z_{k}^{2}$ with
$z_{k}^{1}\in\mathcal{K}(A^{T})$ and $z_{k}^{2}\in\mathcal{R}(A)$. Then, we
have

$x_{k+1}^{*}=A^{T}z_{k}=A^{T}(z_{k}^{1}+z_{k}^{2})=A^{T}z_{k}^{2}.$

Note that $\\{x_{k}^{*}\\}$ is a bounded sequence and $A^{T}$ is one-to-one
from $\mathcal{A}$ to $\mathbb{R}^{n}$. The sequence $\\{z_{k}^{2}\\}$ must be
bounded, that is, there exists a constant $C>0$ such that

$\|z_{k}^{2}\|\leq C,~{}~{}\forall~{}k\geq 0.$

Now, using the subgradient inequality and Cauchy-Schwartz inequality, we
derive that for $k\geq 1$,

$\displaystyle\omega(x_{k})$ $\displaystyle\leq$
$\displaystyle\omega(\bar{x})-\langle A^{T}z_{k-1}^{2},\bar{x}-x_{k}\rangle$
(4.3) $\displaystyle=$ $\displaystyle\omega(\bar{x})-\langle
z_{k-1}^{2},A\bar{x}-Ax_{k}\rangle$ $\displaystyle\leq$
$\displaystyle\omega(\bar{x})+\|z_{k-1}^{2}\|\cdot\|A\bar{x}-Ax_{k}\|$
$\displaystyle\leq$ $\displaystyle\omega(\bar{x})+C\cdot\|A\bar{x}-Ax_{k}\|.$

In what follows, under the two settings of $g$ and $A$, we will show that

$A\hat{x}=A\bar{x},~{}~{}\forall~{}\hat{x}\in\overline{X},$ (4.4)

which immediately implies that

$\|A\bar{x}-Ax_{k}\|=\|A\tilde{x}-Ax_{k}\|\leq\|A\|\cdot\|\tilde{x}-x_{k}\|\rightarrow
0,$

as $k\rightarrow\infty$ since $x_{k}\rightarrow\tilde{x}$ as
$k\rightarrow\infty$. Together with (4.3), we get
$\omega(\tilde{x})=\lim_{k\rightarrow\infty}\omega(x_{k})\leq\omega(\bar{x}),$
which means that $\tilde{x}$ is also a minimizer of the bilevel optimization
problem (3.3) and hence equals to $\bar{x}$ by the uniqueness of solution. In
other words, the sequence $\\{x_{k}\\}$ converges to the optimal solution
$\bar{x}$ of the bilevel optimization problem (3.3). Therefore, it remains to
show (4.4). In the case (i), we denote the unique minimizer of $g$ by
$\bar{y}$. For any $\hat{x}\in\overline{X}$, we have

$0=\nabla f(\hat{x})=A^{T}\nabla g(A\hat{x}),$

which implies that $\nabla g(A\hat{x})=0$ due to the fact of $A$ being
surjective. Thus, $A\hat{x}$ is a global minimizer of $g$ due to the convexity
of $g$ and hence it must equal to the unique minimizer $\bar{y}$, that is,
$A\hat{x}=\bar{y}$. Therefore, the result (4.4) holds in the first case. In
the case (ii), we will use the following property implied by the strict
convexity of $g$:

$\langle\nabla g(u)-\nabla g(v),u-v\rangle>0,\forall u\neq v.$ (4.5)

Suppose that there exists $\hat{x}\in\overline{X}$ such that $A\hat{x}\neq
A\bar{x}$; then by (4.5), we have

$\langle\nabla g(A\hat{x})-\nabla g(A\bar{x}),A\hat{x}-A\bar{x}\rangle>0.$

On the other hand, we have

$A^{T}\nabla g(A\hat{x})=A^{T}\nabla g(A\bar{x})=0,$

and hence

$\displaystyle\langle\nabla g(A\hat{x})-\nabla
g(A\bar{x}),A\hat{x}-A\bar{x}\rangle$ $\displaystyle=$ $\displaystyle\langle
A^{T}\nabla g(A\hat{x})-A^{T}\nabla g(A\bar{x}),\hat{x}-\bar{x}\rangle=0.$

It is a contradiction. This completes the proof. ∎

Condition (i) in Theorem 4.2 is weaker than the corresponding assumption on
$A$ and $g$ in Theorem 4.1 in [30] since the growth property there must imply
the uniqueness of the minimizer of $g$. Regarding condition (ii), it is
general enough to apply to the first example in Section 3 because the inner-
level objective function $f(x)=\frac{1}{2}\|Ax-b\|^{2}$ can be written as
$f(x)=g(Ax)$ with

$g(y)=\frac{1}{2}\|y-b\|^{2},$

which is a strictly convex function and satisfies the assumption on $f$ in
Theorem 4.2. Therefore, we have the following corollary, which is the main
contribution made in a series of papers [10, 9, 8].

######  Corollary 4.1.

Consider the bilevel optimization problem (3.3) with
$\omega(x)=\lambda\|x\|_{1}+\frac{1}{2}\|x\|_{2}$ and
$f(x)=\frac{1}{2}\|Ax-b\|^{2}$, where $\lambda>0$ is some fixed parameter,
$A\in\mathbb{R}^{m\times n}$ is a given nonzero matrix, and
$b\in\mathbb{R}^{m}$ is a given vector. Let $\\{x_{k}\\}$ be generalized by
the following iterative scheme:

$\displaystyle\left\\{\begin{array}[]{rll}x_{k+1}&=&\mathcal{S}_{\lambda}(x_{k}^{*}-t_{k}\cdot
A^{T}(Ax_{k}-b)),\\\ x_{k+1}^{*}&=&x_{k}^{*}-t_{k}\cdot
A^{T}(Ax_{k}-b),\end{array}\right.$ (4.8)

initialized with $x_{0}=x_{0}^{*}=0$, where
$\mathcal{S}_{\lambda}(x)=\min(|x|-\lambda,0)\textrm{sign}(x)$ is the
componentwise soft shrinkage, and $t_{k}$ is a constant in
$\left(0,\frac{2}{\|AA^{T}\|}\right)$ or an exact solution to the dual problem

$\operatorname*{minimize\quad}_{t>0}g(t)=\frac{1}{2}\|\mathcal{S}_{\lambda}(x_{k}^{*}-t\cdot
A^{T}(Ax_{k}-b))\|^{2}+t\cdot\beta_{k},$ (4.9)

with $\beta_{k}=\langle
A^{T}(Ax_{k}-b),x_{k}\rangle-\frac{1}{\|AA^{T}\|}\|A^{T}(Ax_{k}-b)\|^{2}$.
Then, the sequence $\\{x_{k}\\}$ converges to the optimal solution $\bar{x}$
of the bilevel optimization problem (3.3).

The proof follows directly from Theorem 4.2 and the fact that

$\nabla\omega^{*}(x^{*})=\mathcal{S}_{\lambda}(x^{*}),~{}~{}\omega^{*}(x^{*})=\frac{1}{2}\|\mathcal{S}_{\lambda}(x^{*})\|^{2}.$

The iterative scheme (4.8) is precisely the linearized Bregman method up to
the step-size $t_{k}$; please see (63d)-(63e) in [18]. In original linearized
Bregman methods, the step-size is usually set to be the constant $1/\|A\|^{2}$
to deduce convergence; see [10]. Recently, it has been allowed to be the exact
solution to the dual problem (4.9), but in this setting the convergence to the
optimal solution $\bar{x}$ is absent; see [20]. Our result in the corollary
covers and supplements these existing studies.

Let us turn to the problem (1.3). For simplicity, we restrict our attention to
the case of $N=1$, that is,

$f(x)=\frac{1}{2}{\mathbf{dist}}^{2}(Ax,Q),$

where $Q\subset\mathbb{R}^{m}$ is a given closed convex set. This form is
general enough to cover the different noise models

$\operatorname*{minimize\quad}_{x}\omega(x),~{}\mbox{subject
to}~{}~{}\|Ax-b\|_{\diamond}\leq\sigma,$ (4.10)

where the choice of the norm $\|\cdot\|_{\diamond}$ depends on the
characteristics of noise involved in measurement vector $b$; e.g., we use
$\ell_{2}$-norm for Gaussian noise, $\ell_{1}$-norm for impulsive noise and
$\ell_{\infty}$-norm for uniformly distributed noise. Let
$Q:=\\{y\in\mathbb{R}^{m}:\|y-b\|_{\diamond}\leq\sigma\\}$; then the
constraint $\|Ax-b\|_{\diamond}\leq\sigma$ can be written as
$f(x)={\mathbf{dist}}^{2}(Ax,Q)$ and hence the model (4.10) can be written as
a bilevel convex optimization (3.3) so that our proposed Algorithm 1 can be
applied. In particular, the model (4.10) with
$\omega(x)=\lambda\|x\|_{1}+\frac{1}{2}\|x\|_{2}$ was discussed in [20], with
the following iterative scheme

$\displaystyle\left\\{\begin{array}[]{rll}x_{k+1}&=&\mathcal{S}_{\lambda}(x_{k}^{*}-t_{k}A^{T}(Ax_{k}-\mathcal{P}_{Q}(Ax_{k}))),\\\
x_{k+1}^{*}&=&x_{k}^{*}-t_{k}A^{T}(Ax_{k}-\mathcal{P}_{Q}(Ax_{k})),\end{array}\right.$
(4.13)

where $\mathcal{P}_{Q}$ is the projection operator onto $Q$. This scheme can
be deduced from our proposed Algorithm 1 as well and hence its convergence to
a feasible point and a sublinear rate of convergence for the inner-level
function can be established by Theorem 4.1 and Theorem 4.3 respectively.
However, the answer to the question of whether the iterates $\\{x_{k}\\}$
converges to the optimal solution $\bar{x}$ is negative even when the matrix
$A$ is surjective. On one hand, the function $f(x)=g(Ax)$ with
$g(y)=\frac{1}{2}{\mathbf{dist}}^{2}(y,Q)$ fails to satisfy the conditions (i)
and (ii) in Theorem 4.2; on the other hand, we give numerical tests to
illustrate that the iterates $\\{x_{k}\\}$ may converge to a different point
from the optimal solution $\bar{x}$.

###  4.3 Linear and sublinear convergence

######  Theorem 4.3.

Let $\\{x_{k}\\}$ be generated by Algorithm 1 with $t_{k}$ belonging to
$(0,\frac{2\mu}{L}]$. Then the sequence of function values $\\{f(x_{k})\\}$ is
monotonically nonincreasing. If the step-size $t_{k}$ is restricted to the
interval $(0,\frac{\mu}{L}]$, then the sequence of function values
$\\{f(x_{k})\\}$ monotonically converges to the optimal function value
$\bar{f}:=\min f(x)$ sublinearly in the sense that

$f(x_{T+1})-\bar{f}\leq\frac{\inf_{x\in\overline{X}}D_{\omega}^{x_{0}^{*}}(x,x_{0})}{\sum_{k=0}^{T}t_{k}}.$
(4.14)

If the step-size $t_{k}$ is allowed to belong to a larger interval
$(0,\frac{2\mu}{L})$ and the out-level objective $\omega$ is assumed to be
gradient-Lipschtiz-continuous with constant $\nu>0$, then we have that

$f(x_{T+1})-\bar{f}\leq\frac{f(x_{0})-\bar{f}}{1+(f(x_{0})-\bar{f})\cdot\sum_{k=0}^{T}h(t_{k})},$
(4.15)

where $h(\tau)=\frac{\tau(2\mu-L\tau)\mu}{4\nu\cdot c(x_{0})}$ with
$c(x_{0})=\inf_{x\in\overline{X}}D_{\omega}^{x_{0}^{*}}(x,x_{0}).$

###### Proof.

By the three-point identity (2.3) of Bregman distance in Lemma 2.2, we can
derive that

$\displaystyle
D_{\omega}^{x_{k}^{*}}(x,x_{k})-D_{\omega}^{x_{k+1}^{*}}(x,x_{k+1})-D_{\omega}^{x_{k}^{*}}(x_{k+1},x_{k})$
$\displaystyle=$ $\displaystyle\langle x-x_{k+1},x_{k+1}^{*}-x_{k}^{*}\rangle$
$\displaystyle=$ $\displaystyle-t_{k}\langle x-x_{k+1},\nabla f(x_{k})\rangle$
$\displaystyle=$ $\displaystyle-t_{k}(\langle x-x_{k},\nabla
f(x_{k})\rangle+\langle x_{k}-x_{k+1},\nabla f(x_{k})\rangle).$

Using the subgradient inequality, we have

$\langle x-x_{k},\nabla f(x_{k})\rangle\leq f(x)-f(x_{k}).$

Using the Lipschitz continuity of $\nabla f$, we have

$\langle x_{k}-x_{k+1},\nabla f(x_{k})\rangle\leq
f(x_{k})-f(x_{k+1})+\frac{L}{2}\|x_{k+1}-x_{k}\|^{2}.$

Therefore, the term
$D_{\omega}^{x_{k}^{*}}(x,x_{k})-D_{\omega}^{x_{k+1}^{*}}(x,x_{k+1})-D_{\omega}^{x_{k}^{*}}(x_{k+1},x_{k})$
can be bounded as follows:

$D_{\omega}^{x_{k}^{*}}(x,x_{k})-D_{\omega}^{x_{k+1}^{*}}(x,x_{k+1})-D_{\omega}^{x_{k}^{*}}(x_{k+1},x_{k})\geq
t_{k}(f(x_{k+1})-f(x))-\frac{Lt_{k}}{2}\|x_{k+1}-x_{k}\|^{2}.$ (4.16)

Letting $x=x_{k}$ in (4.16) and using the relationship (2.5) in Lemma 2.2, we
can derive that

$\displaystyle t_{k}(f(x_{k+1})-f(x_{k}))$ $\displaystyle\leq$ $\displaystyle-
D_{\omega}^{x_{k}^{*}}(x_{k+1},x_{k})-D_{\omega}^{x_{k+1}^{*}}(x_{k},x_{k+1})+\frac{Lt_{k}}{2}\|x_{k+1}-x_{k}\|^{2}$
(4.17) $\displaystyle\leq$
$\displaystyle-\mu\|x_{k+1}-x_{k}\|^{2}+\frac{Lt_{k}}{2}\|x_{k+1}-x_{k}\|^{2}$
$\displaystyle=$
$\displaystyle\left(\frac{Lt_{k}}{2}-\mu\right)\|x_{k+1}-x_{k}\|^{2}.$

If the step-size $t_{k}\in(0,\frac{2\mu}{L}]$, then
$\left(\frac{Lt_{k}}{2}-\mu\right)\leq 0$ and hence $f(x_{k+1})\leq f(x_{k})$
which means that the sequence of function values $\\{f(x_{k})\\}$ is
monotonically nonincreasing.

Now, we return to (4.16) and derive that for $t_{k}\in(0,\frac{\mu}{L}]$,

$\displaystyle t_{k}(f(x_{k+1})-f(x))$ $\displaystyle\leq$ $\displaystyle
D_{\omega}^{x_{k}^{*}}(x,x_{k})-D_{\omega}^{x_{k+1}^{*}}(x,x_{k+1})-D_{\omega}^{x_{k}^{*}}(x_{k+1},x_{k})+\frac{Lt_{k}}{2}\|x_{k+1}-x_{k}\|^{2}$
$\displaystyle\leq$ $\displaystyle
D_{\omega}^{x_{k}^{*}}(x,x_{k})-D_{\omega}^{x_{k+1}^{*}}(x,x_{k+1})+\left(\frac{Lt_{k}}{2}-\frac{\mu}{2}\right)\|x_{k+1}-x_{k}\|^{2}$
$\displaystyle\leq$ $\displaystyle
D_{\omega}^{x_{k}^{*}}(x,x_{k})-D_{\omega}^{x_{k+1}^{*}}(x,x_{k+1}).$

Summing the inequality above from $k=0$ to $k=T$, we obtain

$\sum_{k=0}^{T}t_{k}(f(x_{k+1})-f(x))\leq
D_{\omega}^{x_{0}^{*}}(x,x_{0})-D_{\omega}^{x_{T+1}^{*}}(x,x_{T+1})\leq
D_{\omega}^{x_{0}^{*}}(x,x_{0}).$ (4.18)

By the monotonicity of $\\{f(x_{k})\\}$, we can get

$(f(x_{T+1})-f(x))\sum_{k=0}^{T}t_{k}\leq D_{\omega}^{x_{0}^{*}}(x,x_{0}),$
(4.19)

which implies the desired result (4.14) by noting that
$\inf_{x\in\overline{X}}(f(x_{T+1})-f(x))=f(x_{T+1})-\bar{f}.$

Below, let us show (4.15). Since the out-level objective $\omega$ is
additionally assumed to be gradient-Lipschtiz-continuous with constant
$\nu>0$, we can derive that

$\displaystyle\|x_{k+1}-x_{k}\|$ $\displaystyle\geq$
$\displaystyle\nu^{-1}\|\nabla\omega(x_{k+1})-\nabla\omega x_{k}\|$
$\displaystyle=$ $\displaystyle\nu^{-1}\|x_{k+1}^{*}-x_{k}^{*}\|$
$\displaystyle=$ $\displaystyle\nu^{-1}t_{k}\|\nabla f(x_{k})\|.$

Combining this with the inequality (4.17) and noting the step-size condition
$t_{k}\in(0,\frac{2\mu}{L})$, we can get

$f(x_{k+1})-f(x_{k})\leq
t_{k}\nu^{-2}\left(\frac{Lt_{k}}{2}-\mu\right)\|\nabla f(x_{k})\|^{2}.$ (4.20)

On the other hand, using the subgradient inequality, the Cauchy-Schwartz
inequality, and the strong convexity of $\omega$, we derive that for any
$x\in\overline{X}$,

$\displaystyle f(x_{k})-f(x)$ $\displaystyle\leq$ $\displaystyle\langle\nabla
f(x_{k}),x_{k}-x\rangle$ (4.21) $\displaystyle\leq$ $\displaystyle\|\nabla
f(x_{k})\|\cdot\|x_{k}-x\|$ $\displaystyle\leq$ $\displaystyle\|\nabla
f(x_{k})\|\cdot\sqrt{\frac{2}{\mu}D_{\omega}^{x_{k}}(x,x_{k})}$
$\displaystyle\leq$ $\displaystyle\|\nabla
f(x_{k})\|\cdot\sqrt{\frac{2}{\mu}D_{\omega}^{x_{0}}(x,x_{0})},$

where the last inequality follows from the descent property (3.27) in Lemma
3.1. Thus, we can get

$f(x_{k})-\bar{f}\leq\|\nabla
f(x_{k})\|\cdot\sqrt{\frac{2}{\mu}\inf_{x\in\overline{X}}D_{\omega}^{x_{0}}(x,x_{0})}=\sqrt{\frac{2c(x_{0})}{\mu}}\cdot\|\nabla
f(x_{k})\|.$ (4.22)

For simplicity of deduction, we let $f_{k}:=f(x_{k})-\bar{f}$. Then, combining
(4.20) and (4.22) and using the definition of $h(\cdot)$, we have

$f_{k+1}-f_{k}\leq
t_{k}\nu^{-2}\left(\frac{Lt_{k}}{2}-\mu\right)\frac{\mu}{2c(x_{0})}\cdot
f_{k}^{2}=-h(t_{k})\cdot f_{k}^{2}.$ (4.23)

Thus, it holds that

$\frac{1}{f_{k+1}}\geq\frac{1}{f_{k}}+h(t_{k})\cdot\frac{f_{k}}{f_{k+1}}\geq\frac{1}{f_{k}}+h(t_{k}),$
(4.24)

where the last inequality follows from the monotonicity of $\\{f(x_{k})\\}$.
Summing (4.24) from $k=0$ to $T$, we get

$\frac{1}{f_{T+1}}\geq\frac{1}{f_{0}}+\sum_{k=0}^{T}h(t_{k}).$ (4.25)

Therefore, we have

$f_{T+1}\leq\frac{f_{0}}{1+f_{0}\cdot\sum_{k=0}^{T}h(t_{k})}.$ (4.26)

This completes the proof. ∎

Note that when $\omega(x)=\frac{1}{2}\|x\|^{2}$, $\nabla\omega^{*}$ is a
identity mapping and Algorithm 1 reduce to the gradient descent method. In
this case, the allowed step-size interval and the sublinear convergence result
are consistent with the classical result [7] for the gradient descent method.
This implies that (4.15) in Theorem 4.3 provides a tight bound.

In order to derive rates of linear convergence, we introduce the following
Bregman distance growth condition, which was independently proposed in [29]
and [21] recently.

######  Definition 4.1.

Let $\overline{X}$ be the minimizer set of $f$ and $\bar{f}$ be its optimal
function value. We say that $f$ satisfies the Bregman distance growth
condition if there exists a constant $\gamma>0$ such that

$f(x)-\bar{f}\geq\gamma\cdot\inf_{z\in\overline{X}}D_{\omega}^{x^{*}}(z,x).$
(4.27)

######  Definition 4.2.

Let $\overline{X}$ be the minimizer set of $f$ and $\bar{f}$ be its optimal
function value. We say that $f$ satisfies the quadratic growth condition if
there exists $\eta>0$ such that for all $x\in\mathbb{R}^{n}$ we have

$f(x)-\bar{f}\geq\frac{\eta}{2}{\mathbf{dist}}^{2}(x,\overline{X}).$ (4.28)

######  Theorem 4.4 (linear convergence).

Let $\\{x_{k}\\}$ be generated by Algorithm 1 with $t_{k}$ being an exact
solution to the dual problem (3.15) or being the constant $\frac{\mu}{L}$.
Assume that $f$ satisfies the Bregman distance growth condition (4.27). Then,
the sequence $\\{x_{k}\\}$ converges to the feasible set $\overline{X}$
linearly in the following sense

$\inf_{x\in\overline{X}}D_{\omega}^{x_{k}^{*}}(x,x_{k})\leq\left(1-\left(\frac{\mu\cdot\gamma}{2L}\right)^{2}\right)^{k}\inf_{x\in\overline{X}}D_{\omega}^{x_{0}^{*}}(x,x_{0}).$
(4.29)

###### Proof.

Starting with the following basic descent property stated in Lemma 3.1,

$D_{\omega}^{x_{k+1}^{*}}(x,x_{k+1})\leq
D_{\omega}^{x_{k}^{*}}(x,x_{k})-\frac{\mu}{2L^{2}}\|\nabla
f(x_{k})\|^{2},~{}\forall x\in H_{k},$

we can get

$\inf_{x\in\overline{X}}D_{\omega}^{x_{k+1}^{*}}(x,x_{k+1})\leq\inf_{x\in\overline{X}}D_{\omega}^{x_{k}^{*}}(x,x_{k})-\frac{\mu}{2L^{2}}\|\nabla
f(x_{k})\|^{2},~{}\forall x\in H_{k},$ (4.30)

where we have used the fact $\overline{X}\subset H_{k}$. Using the Bregman
distance growth condition (4.27) and the strong convexity of $\omega$, we can
derive that

$\displaystyle f(x)-\bar{f}$ $\displaystyle\geq$
$\displaystyle\gamma\cdot\inf_{z\in\overline{X}}D_{\omega}^{x^{*}}(z,x)$
(4.31) $\displaystyle\geq$
$\displaystyle\gamma\cdot\inf_{z\in\overline{X}}\frac{\mu}{2}\|z-x\|^{2}$
$\displaystyle=$
$\displaystyle\frac{\gamma\cdot\mu}{2}{\mathbf{dist}}^{2}(x,\overline{X}).$

In other words, the inner objective function $f$ satisfies the quadratic
growth condition and hence it also satisfies the PL inequality

$\|\nabla f(x)\|\geq\tau\sqrt{f(x)-\bar{f}},$ (4.32)

where $\tau=\sqrt{\frac{\gamma\cdot\mu}{2}}$; see the interplay theorem in
[28]. Now, combining (4.32), (4.30), and (4.27), we derive that

$\displaystyle\inf_{x\in\overline{X}}D_{\omega}^{x_{k+1}^{*}}(x,x_{k+1})$
$\displaystyle\leq$
$\displaystyle\inf_{x\in\overline{X}}D_{\omega}^{x_{k}^{*}}(x,x_{k})-\frac{\mu}{2L^{2}}\tau^{2}(f(x_{k})-\bar{f})$
(4.33) $\displaystyle\leq$
$\displaystyle\inf_{x\in\overline{X}}D_{\omega}^{x_{k}^{*}}(x,x_{k})-\frac{\mu}{2L^{2}}\tau^{2}\cdot\gamma\inf_{x\in\overline{X}}D_{\omega}^{x_{k}^{*}}(x,x_{k})$
$\displaystyle=$
$\displaystyle\left(1-\left(\frac{\mu\cdot\gamma}{2L}\right)^{2}\right)\inf_{x\in\overline{X}}D_{\omega}^{x_{k}^{*}}(x,x_{k}).$

This completes the proof. ∎

The key to obtaining linear convergence is the Bregman distance growth
condition. Next, we will give the validation of the condition on problem
(1.2). We first give the following result from [25].

######  Lemma 4.2 ([25], Lemma 3.1).

Suppose that the linear system $Ax=b$ is consistent. Let
$\omega(x)=\lambda\|x\|_{1}+\frac{1}{2}\|x\|^{2}$, where $\lambda>0$ is fixed
parameter. Denote

$\bar{x}=\arg\min_{x}\omega(x),~{}~{}\mbox{subject to}~{}Ax=b.$

Then, there exists a constant $\tau>0$ such that for all $x\in\mathbb{R}^{n}$
with $\partial\omega(x)\bigcap\mathcal{R}(A^{T})\neq\emptyset$ and
$x^{*}\in\partial\omega(x)\bigcap\mathcal{R}(A^{T})\neq\emptyset$ we have

$D_{\omega}^{x^{*}}(\bar{x},x)\leq\tau\cdot\|Ax-b\|^{2}.$ (4.34)

######  Corollary 4.2.

Under the same setting of Corollary 4.1 and assuming that $A^{T}b\neq 0$, then
$f$ satisfies the Bregman distance growth condition with a constant
$\frac{\eta}{2\tau L^{2}}$. Moreover, the sequence $\\{x_{k}\\}$ converges to
the feasible set $\overline{X}$ linearly in the following sense that

$\inf_{x\in\overline{X}}D_{\omega}^{x_{k}^{*}}(x,x_{k})\leq\left(1-\left(\frac{\mu\eta}{4\tau
L^{3}}\right)^{2}\right)^{k}\inf_{x\in\overline{X}}D_{\omega}^{x_{0}^{*}}(x,x_{0}).$
(4.35)

###### Proof.

By the first-order optimality condition, the feasible set $\overline{X}$ can
be written as

$\overline{X}=\\{x\in\mathbb{R}^{n}:A^{T}(Ax-b)=0\\}.$

Thus, the corresponding bilevel convex optimization problem (3.3) can be
written as

$\bar{x}=\arg\min_{x}\omega(x),~{}~{}\mbox{subject to}~{}A^{T}Ax=A^{T}b.$

Using the basic fact that $\mathcal{R}(A^{T}A)=\mathcal{R}(A^{T})$, we have
that $A^{T}b\neq 0$ and $A^{T}b$ lies in $\mathcal{R}(A^{T}A)$. Therefore,
invoking Lemma 4.2 and using the fact $\mathcal{R}(A^{T}A)=\mathcal{R}(A^{T})$
again, we conclude that there exists a constant $\tau>0$ such that for all
$x\in\mathbb{R}^{n}$ with
$\partial\omega(x)\bigcap\mathcal{R}(A^{T})\neq\emptyset$ and
$x^{*}\in\partial\omega(x)\bigcap\mathcal{R}(A^{T})\neq\emptyset$ we have

$D_{\omega}^{x^{*}}(\bar{x},x)\leq\tau\cdot\|A^{T}Ax-A^{T}b\|^{2}.$ (4.36)

With this inequality, we further derive that

$\inf_{z\in\overline{X}}D_{\omega}^{x^{*}}(z,x)\leq
D_{\omega}^{x^{*}}(\bar{x},x)\leq\tau\cdot\|A^{T}Ax-A^{T}b\|^{2}=\tau\cdot\|\nabla
f(x)\|^{2}.$ (4.37)

Using the gradient Lipschitz continuity, we have that

$\|\nabla f(x)\|=\|\nabla f(x)-\nabla f(\hat{x})\|\leq
L\cdot\|x-\hat{x}\|,~{}~{}\forall~{}\hat{x}\in\overline{X}.$

Since that $f$ has the form of $g(Ax)$, it follows from [22] that $f$
satisfies the quadratic growth condition, i.e., there exists a constant $\eta$
such that (4.28) holds. This implies that

$\|\nabla
f(x)\|\leq\inf_{\hat{x}\in\overline{X}}L\cdot\|x-\hat{x}\|=L\cdot{\mathbf{dist}}(x,\overline{X})\leq
L\cdot\sqrt{\frac{2}{\eta}(f(x)-\bar{f})}.$

Therefore, combining with (4.37), we obtain that

$\inf_{z\in\overline{X}}D_{\omega}^{x^{*}}(z,x)\leq\frac{2\tau
L^{2}}{\eta}(f(x)-\bar{f}).$

In other words, $f$ satisfies the Bregman distance growth condition and hence
the conclusion follows from Theorem 4.4. This completes the proof. ∎

In [25], the authors provide a similar linear convergence result for the
sparse Kaczmarz method, which is equivalent to the iterative scheme (4.8). It
should be emphasized that they require that the linear system $Ax=b$ is
consistent, as shown in Lemma 4.2. In our case, we drop this condition since
we apply Lemma 4.2 with the linear system $A^{\top}Ax=A^{\top}b$, which is
always consistent. In contrast, we require that $A^{\top}b\neq 0$, this
condition is more checkable compared with the condition in [25].

##  5 Experiment results

In this section, we conduct four numerical experiments to illustrate different
aspects of Algorithm 1. In Section 5.1 we investigate the use of the different
step-sizes with a compressed sensing example. In Section 5.2 we show that it
is simple to adapt Algorithm 1 to feasibility problems by defining the
distance function. In particular, we consider several sparse recovery problems
under different types of noise, including Gaussian noise and uniform noise.
Moreover, we will also answer the question of whether Algorithm 1 computes an
optimal solution when the constraint does not satisfy the condition in Theorem
4.2.

###  5.1 Linear inverse problem

We conduct an experiment on the comparison of step-size rules for solving
problem (1.2). The gradient Lipschitz constant is $L:=\|A^{\top}A\|_{2}$. Our
algorithm has the following iterative scheme:

$\left\\{\begin{aligned}
x_{k+1}&=\mathcal{S}_{\lambda}(x^{*}_{k}-t_{k}A^{\top}(Ax_{k}-b)),\\\
x_{k+1}^{*}&=x^{*}_{k}-t_{k}A^{\top}(Ax_{k}-b),\end{aligned}\right.$ (5.1)

where $t_{k}$ is the step-size, which can be chosen by the following strategy:

  * •

The exact step-size by the exact solution of the following dual problem

$\operatorname*{minimize\quad}_{t>0}\frac{1}{2}\|\mathcal{S}_{\lambda}(x^{*}_{k}-tA^{\top}(Ax_{k}-b))\|^{2}+t\cdot\beta_{k},$
(5.2)

where $\beta_{k}=\langle
A^{\top}(Ax_{k}-b),x_{k}\rangle-\frac{\|A^{\top}(Ax_{k}-b)\|^{2}}{L}$.

  * •

The constant step-size $t_{k}=\frac{\mu}{L}=\frac{1}{L}$.

  * •

The dynamic step-size
$t_{k}=\frac{\|Ax_{k}-b\|_{2}^{2}}{\|A^{\top}(Ax_{k}-b)\|_{2}^{2}}$.

In order to obtain an optimal solution we use the primal-dual method [3]. The
method considers the following general optimization problem:

$\operatorname*{minimize\quad}_{x\in\mathbb{R}^{n}}g(x)+h(Ax),$ (5.3)

where $g$ is strongly convex, closed, and proper. $h$ is properly closed and
proximable. Employs the fast dual proximal gradient method whose general
update step is

$\left\\{\begin{aligned} u_{k}&=\operatorname{argmax}_{u}\left\\{\left\langle
u,Aw_{k}\right\rangle-g(u)\right\\},\\\
y_{k+1}&=w_{k}-\frac{1}{L_{k}}Au_{k}+\frac{1}{L_{k}}\operatorname{prox}_{\lambda
L_{k}h}\left(Au_{k}-L_{k}w_{k}\right),\\\
t_{k+1}&=\frac{1+\sqrt{1+4t_{k}^{2}}}{2},\\\
w_{k+1}&=y_{k+1}+\left(\frac{t_{k}-1}{t_{k+1}}\right)\left(y_{k+1}-y_{k}\right).\end{aligned}\right.$
(5.4)

where $L_{k}$ is chosen either as a constant or by a backtracking procedure
(default). The last computed vector $u_{k}$ is considered as the obtained
solution and the last computed $y_{k}$ is the obtained dual vector. We apply
the algorithm to solve (1.2) by letting
$g(x)=\lambda\|x\|_{1}+\frac{1}{2}\|x\|_{2}^{2}$ and $h(x)=\delta_{Q}(x)$,
where $Q=\\{x~{}:x=b\\}$.

For a given matrix $A$ we construct a recoverable sparse vector $x^{\dagger}$
using L1TestPack [19], compute the corresponding right-hand side
$b=Ax^{\dagger}$, and choose the regularization parameter
$\lambda=\|x^{\dagger}\|_{1}$. The results are summarized in Figure 1. We
denote Algorithm 1 with dynamic, exact, and constant step-sizes as “Breg-dy”,
“Breg-opt” and “Breg-con”, respectively. The primal-dual method (5.4) is
referred to as “primal-dual”. The reconstruction error is defined as
$\|x_{k}-x^{\dagger}\|$, where $x_{k}$ is the iterate in Algorithm. The
objective value denotes $\lambda\|x\|_{1}+\frac{1}{2}\|x\|_{2}^{2}$. The
feasibility is denoted by $\|Ax-b\|_{2}$. Figure 1 (c) shows that All step-
sizes have successfully attained feasible solutions, with the exact step size
demonstrating linear convergence from the iterative point to the feasible
solution, thereby validating Theorem 4.4. Furthermore, all step-sizes have
identified the optimal solution for problem (1.2), as illustrated in Figure 1
(b). This is attributed to the fulfillment of the conditions shown in Theorem
4.2 for problem (1.2). Regarding the reconstruction error, it is observed that
the exact step size exhibits a smaller reconstruction error compared to other
step-sizes.

(a) Reconstruction error ($\log$) (b) Objective value (c) Feasibility ($\log$)

Figure 1: Comparison with different step-size for least square problem

(a) True (b) Constant (c) Optimal (d) Primal-dual

Figure 2: Comparison for least square problem

(a) Reconstruction error ($\log$) (b) Objective value (c) Feasibility ($\log$)

Figure 3: Comparison for the $\ell_{2}$ regularization problem

###  5.2 Feasibility problems

This subsection considers problem (1.3). In particular, consider the following
feasibility problem:

$\operatorname*{minimize\quad}_{x\in\mathbb{R}^{n}}\omega(x),\;\text{s.t.}\;Ax\in
Q,$ (5.5)

where $\omega$ is strongly convex and $Q$ is an convex set. Denote
$f(x):=\frac{1}{2}\mathbf{dist}^{2}(Ax,Q)$. Then one can construct an
equivalent bilevel optimization problem as follows:

$\operatorname*{minimize\quad}_{x\in\mathbb{R}^{n}}\omega(x),\;\text{s.t.}\;x\in\arg\min_{x}f(x)=\frac{1}{2}\mathbf{dist}^{2}(Ax,Q).$
(5.6)

The gradient Lipschitz constant of $f$ is $L:=\|A^{\top}A\|_{2}$. Our
algorithm has the following iterative scheme:

$\left\\{\begin{aligned}
x_{k+1}&=\mathcal{S}_{\lambda}(x^{*}_{k}-t_{k}A^{\top}(Ax_{k}-\mathcal{P}_{Q}(Ax_{k}))),\\\
x_{k+1}^{*}&=x^{*}_{k}-t_{k}A^{\top}(Ax_{k}-\mathcal{P}_{Q}(Ax_{k})).\end{aligned}\right.$
(5.7)

In order to illustrate the performance of Algorithm 1 to solve problem (5.6),
we consider the problem of recovering sparse solutions of linear equations
$Ax=b$, where only noisy data $b^{\sigma}$ is available for different noise
models. In this case, the convex set
$Q:=\\{x~{}:\|x-b^{\sigma}\|\leq\sigma\\}$, where $\sigma$ denote the noise.
The choice of the norm $\|\cdot\|$ is dictated by the noise characteristics.
In the following, we will consider three types of norms, i.e., the
$\ell_{2}$-norm for Gaussian noise and the $\ell_{\infty}$-norm for uniformly
distributed noise. In order to obtain the optimal solution, we also run the
primal-dual algorithm to solve (5.6) by letting
$g(x)=\lambda\|x\|_{1}+\frac{1}{2}\|x\|_{2}^{2}$ and $g(x)=\delta_{Q}(x)$.

####  5.2.1 Sparse recovery with Gaussian noise

This subsection focus on an special case of problem (5.6), where
$Q=\\{y\in\mathbb{R}^{m}~{}:~{}\|y-b^{\sigma}\|_{2}\leq\sigma\\}$, $\sigma$ is
the noise. In this case

$Ax-\mathcal{P}_{Q}(Ax)=\max\left\\{0,1-\frac{\sigma}{\left\|Ax-b^{\sigma}\right\|_{2}}\right\\}\cdot\left(Ax-b^{\sigma}\right).$
(5.8)

For some matrix $A\in\mathbb{R}^{m\times n}$ we produce a sparse vector
$x^{\dagger}$ with only 30 nonzero entries and calculate $b=Ax^{\dagger}$. We
choose $b^{\sigma}$ by adding the Gaussian noise to $b$, and
$\sigma=\|b-b^{\sigma}\|_{2}$.

Figure 3 shows the result. Despite the successful identification of feasible
solutions for all step-sizes, Figure 3 (b) illustrates that they yield larger
objective function values compared to the primal-dual algorithm. This
indicates that our algorithm may not necessarily find the optimal solution.
The reason lies in the failure of problem (5.6) to satisfy the conditions in
Theorem 4.2, thereby indirectly validating the reasonableness of these
conditions. Despite the absence of an optimal solution, Figure 3 (a)
demonstrates that, under the exact step size, our algorithm accurately
recovers the original solution. This further underscores the practical
significance of our algorithm. In order to show the effect of the level of
noise, we let the noise constant $\sigma=c\|b-b^{\sigma}\|_{2}$, where
$c=0.1,0.5,1$. Figure 4 shows that a low level of noise achieves a lower
reconstruction error.

(a) Reconstruction error ($\log$) (b) Objective value (c) Feasibility ($\log$)

Figure 4: Comparison with different noise level for the $\ell_{2}$
regularization problem

####  5.2.2 Sparse recovery with uniform noise

This subsection focus on an special case of problem (5.6), where
$Q=\\{y\in\mathbb{R}^{m}~{}:~{}\|y-b\|_{\infty}\leq\sigma\\}$, $\sigma$ is the
noise. In this case

$Ax-\mathcal{P}_{Q}(Ax)=\mathcal{S}_{\sigma}(Ax-b^{\sigma}).$ (5.9)

For some matrix $A\in\mathbb{R}^{m\times n}$ we produce a sparse vector
$x^{\dagger}$ with only 30 nonzero entries and calculate $b=Ax^{\dagger}$. We
choose $b^{\sigma}$ by adding uniformly distributed noise with range $[-1,1]$,
and $\sigma=\|b-b^{\sigma}\|_{\infty}$. Figure 5 shows a similar result as the
previous subsection. Our algorithm with the exact step-size achieves better
performance in terms of the reconstruction error and feasibility.

(a) Reconstruction error ($\log$) (b) Objective value (c) Feasibility ($\log$)

Figure 5: Comparison for the $\ell_{\infty}$ regularization problem

##  6 Conclusion

This paper introduces a novel bilevel optimization formulation, offering a
cut-and-project perspective to reevaluate Bregman regularized iterations. This
approach not only encompasses existing methods but also extends to tackle
broader inverse problems, including certain sparse noise models. Through
detailed analysis, we explore the impact of step-sizes on algorithmic
performance, delineating a practical range for their selection and providing a
convergence guarantee to a feasible point. Our unified convergence condition
offers a less restrictive criterion compared to existing literature, enhancing
the applicability of our framework. Additionally, we introduce a Bregman
distance growth condition, proving the linear convergence of our algorithm.

## References

  * [1] Jean-Bernard Baillon and Georges Haddad.  Quelques propriétés des opérateurs angle-bornés et n-cycliquement monotones.  Israel Journal of Mathematics, 26:137–150, 1977.
  * [2] Heinz H Bauschke, Jonathan M Borwein, et al.  Legendre functions and the method of random Bregman projections.  Journal of convex analysis, 4(1):27–67, 1997.
  * [3] Amir Beck.  First-order methods in optimization.  SIAM, 2017.
  * [4] Amir Beck and Shoham Sabach.  A first order method for finding minimal norm-like solutions of convex optimization problems.  Mathematical Programming, 147(1-2):25–46, 2014.
  * [5] Martin Benning, Marta M Betcke, Matthias J Ehrhardt, and Carola-Bibiane Schönlieb.  Choose your path wisely: gradient descent in a Bregman distance framework.  SIAM Journal on Imaging Sciences, 14(2):814–843, 2021.
  * [6] Lev M Bregman.  The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming.  USSR computational mathematics and mathematical physics, 7(3):200–217, 1967.
  * [7] Sébastien Bubeck et al.  Convex optimization: Algorithms and complexity.  Foundations and Trends® in Machine Learning, 8(3-4):231–357, 2015.
  * [8] Jian-Feng Cai, Emmanuel J Candès, and Zuowei Shen.  A singular value thresholding algorithm for matrix completion.  SIAM Journal on optimization, 20(4):1956–1982, 2010.
  * [9] Jian-Feng Cai, Stanley Osher, and Zuowei Shen.  Convergence of the linearized Bregman iteration for $\ell_{1}$-norm minimization.  Mathematics of Computation, 78(268):2127–2136, 2009.
  * [10] Jian-Feng Cai, Stanley Osher, and Zuowei Shen.  Linearized Bregman iterations for compressed sensing.  Mathematics of computation, 78(267):1515–1536, 2009.
  * [11] Jian-Feng Cai, Stanley Osher, and Zuowei Shen.  Linearized bregman iterations for frame-based image deblurring.  SIAM Journal on Imaging Sciences, 2(1):226–252, 2009.
  * [12] Jerome Darbon and Stanley Osher.  Fast discrete optimization for sparse approximations and deconvolutions.  Preprint, 2007.
  * [13] Martin Hanke.  Accelerated landweber iterations for the solution of ill-posed equations.  Numerische mathematik, 60:341–373, 1991.
  * [14] Jean-Baptiste Hiriart-Urruty and Claude Lemaréchal.  Fundamentals of convex analysis.  Springer Science & Business Media, 2004.
  * [15] Stefan Karczmarz.  Angenaherte auflosung von systemen linearer glei-chungen.  Bull. Int. Acad. Pol. Sic. Let., Cl. Sci. Math. Nat., pages 355–357, 1937.
  * [16] Krzysztof C Kiwiel.  Free-steering relaxation methods for problems with strictly convex costs and linear constraints.  Mathematics of Operations Research, 22(2):326–349, 1997.
  * [17] Krzysztof C Kiwiel.  Proximal minimization methods with generalized Bregman functions.  SIAM journal on control and optimization, 35(4):1142–1168, 1997\.
  * [18] Ming-Jun Lai and Wotao Yin.  Augmented $\ell_{1}$ and nuclear-norm models with a globally linearly convergent algorithm.  SIAM Journal on Imaging Sciences, 6(2):1059–1091, 2013.
  * [19] Dirk A Lorenz.  Constructing test instances for basis pursuit denoising.  IEEE Transactions on Signal Processing, 61(5):1210–1214, 2012.
  * [20] Dirk A Lorenz, Frank Schopfer, and Stephan Wenger.  The linearized Bregman method via split feasibility problems: analysis and generalizations.  SIAM Journal on Imaging Sciences, 7(2):1237–1262, 2014.
  * [21] Haihao Lu, Robert M Freund, and Yurii Nesterov.  Relatively smooth convex optimization by first-order methods, and applications.  SIAM Journal on Optimization, 28(1):333–354, 2018.
  * [22] Ion Necoara, Yu Nesterov, and Francois Glineur.  Linear convergence of first order methods for non-strongly convex optimization.  Mathematical Programming, 175:69–107, 2019.
  * [23] Daniel Reem, Simeon Reich, and Alvaro De Pierro.  Re-examination of Bregman functions and new properties of their divergences.  Optimization, 68(1):279–348, 2019.
  * [24] R Tyrrell Rockafellar.  Convex analysis, volume 11.  Princeton university press, 1997.
  * [25] Frank Schöpfer and Dirk A Lorenz.  Linear convergence of the randomized sparse Kaczmarz method.  Mathematical Programming, 173:509–536, 2019.
  * [26] Wotao Yin.  Analysis and generalizations of the linearized Bregman method.  SIAM Journal on Imaging Sciences, 3(4):856–877, 2010.
  * [27] Wotao Yin, Stanley Osher, Donald Goldfarb, and Jerome Darbon.  Bregman iterative algorithms for $\ell_{1}$-minimization with applications to compressed sensing.  SIAM Journal on Imaging sciences, 1(1):143–168, 2008.
  * [28] Hui Zhang.  New analysis of linear convergence of gradient-type methods via unifying error bound conditions.  Mathematical Programming, 180(1):371–416, 2020.
  * [29] Hui Zhang, Yu-Hong Dai, Lei Guo, and Wei Peng.  Proximal-like incremental aggregated gradient method with linear convergence under Bregman distance growth conditions.  Mathematics of Operations Research, 46(1):61–81, 2021.
  * [30] Hui Zhang, Lu Zhang, and Hao-Xing Yang.  Revisiting linearized Bregman iterations under Lipschitz-like convexity condition.  Mathematics of Computation, 92(340):779–803, 2023.