# Two-step inertial Bregman proximal alternating linearized minimization
algorithm for nonconvex and nonsmooth problems††thanks: Supported by
Scientific Research Project of Tianjin Municipal Education Commission
(2022ZD007).

Chenzheng Guo, Jing Zhao
College of Science, Civil Aviation University of China, Tianjin 300300, China
Email<EMAIL_ADDRESS>author. Email<EMAIL_ADDRESS>
Abstract. In this paper, we study an algorithm for solving a class of
nonconvex and nonsmooth nonseparable optimization problems. Based on proximal
alternating linearized minimization (PALM), we propose a new iterative
algorithm which combines two-step inertial extrapolation and Bregman distance.
By constructing appropriate benefit function, with the help of
Kurdyka–Łojasiewicz property we establish the convergence of the whole
sequence generated by proposed algorithm. We apply the algorithm to sparse
nonnegative matrix factorization, signal recovery, quadratic fractional
programming problem and show the effectiveness of proposed algorithm.

Key words: Nonconvex and nonsmooth nonseparable optimization; Proximal
alternating linearized minimization; Inertial extrapolation; Bregman distance;
Kurdyka–Łojasiewicz property.

##  1 Introduction

Since nonconvex functions usually approximate the original problem better than
convex functions, a large number of problems require to solve nonconvex
minimization problems. In recent years, many scholars have paid attention to
study nonconvex optimization problems, and some effective and stable
algorithms have been proposed. In this paper, we will consider solving the
following nonconvex and nonsmooth nonseparable optimization problem:

$\min_{x\in\mathbb{R}^{l},y\in\mathbb{R}^{m}}L(x,y)=f(x)+Q(x,y)+g(y),$ (1.1)

where $f:\mathbb{R}^{l}\rightarrow{(-\infty,+\infty]}$,
$g:\mathbb{R}^{m}\rightarrow{(-\infty,+\infty]}$ are proper lower
semicontinuous,
$Q(x,y):\mathbb{R}^{l}\times\mathbb{R}^{m}\rightarrow\mathbb{R}$ is
continuously differentiable and $\nabla Q$ is Lipschitz continuous on bounded
subsets. Assume that $\inf_{\mathbb{R}^{l}\times\mathbb{R}^{m}}L>-\infty$,
$\inf_{\mathbb{R}^{l}}f>-\infty$, $\inf_{\mathbb{R}^{m}}g>-\infty$. Note that
here and throughout the paper, no convexity is assumed on the objective
function.

Many nonconvex optimization problems have coupling terms in the objective
function, so many application problems can be modeled as (1.1), e.g., signal
recovery and image processing [1, 2, 3], matrix decomposition [4], quadratic
fractional programming [5], compressed sensing [6, 7], applied statistics and
machine learning [8], etc.

Utilizing the two-block structure, a natural method to solve problem (1.1) is
the alternating minimization method, which, from a given initial point
$\left(x_{0},y_{0}\right)\in\mathbb{R}^{l}\times\mathbb{R}^{m}$, generates the
iterative sequence $\left\\{\left(x_{k},y_{k}\right)\right\\}$ via the scheme:

$\begin{cases}x_{k+1}\in\arg\min_{x\in\mathbb{R}^{l}}\\{L(x,y_{k})\\},\\\
y_{k+1}\in\arg\min_{y\in\mathbb{R}^{m}}\\{L(x_{k+1},y)\\}.\\\ \end{cases}$
(1.2)

In the literature, the alternating minimization method are also called the
Gauss–Seidel method or the block coordinate descent method.

The convergence of the alternating minimization method was first studied for
the convex case. If the function $L$ is convex and continuously
differentiable, and it is strict convex of one argument while the other is
fixed, then every limit point of the sequence
$\left\\{\left(x_{k},y_{k}\right)\right\\}$ generated by (1.2) minimizes $L$
[9, 10, 11]. Removing the strict convexity assumption, one can modify the
alternating minimization algorithm by adding a proximal term [12], resulting
in the scheme:

$\begin{cases}x_{k+1}\in\arg\min_{x\in\mathbb{R}^{l}}\\{L(x,y_{k})+\frac{1}{2\lambda_{k}}\|x-x_{k}\|^{2}_{2}\\},\\\
y_{k+1}\in\arg\min_{y\in\mathbb{R}^{m}}\\{L(x_{k+1},y)+\frac{1}{2\mu_{k}}\|y-y_{k}\|^{2}_{2}\\},\\\
\end{cases}$ (1.3)

where $\\{\lambda_{k}\\}_{k\in\mathbb{N}}$ and
$\\{\mu_{k}\\}_{k\in\mathbb{N}}$ are positive sequences. In [13], Attouch et
al. applied (1.3) to solve nonconvex problem (1.1) and proved the sequence
generated by the proximal alternating minimization (PAM) algorithm (1.3)
converges to a critical point. Because the proximal alternating minimization
algorithm requires an exact solution at each iteration step, the subproblems
in (1.3) may be very hard to solve. The linearization technique is one of the
effective methods to overcome the absence of an analytic solution to the
subproblem. To overcome these drawbacks in the above algorithm, by the
proximal linearization of each subproblem in (1.3), Bolte et al. [4] were
inspired by the proximal forward-backward algorithm and proposed the proximal
alternating linearized minimization (PALM) algorithm under the condition that
the coupling term $Q(x,y)$ is continuously differentiable. That is to say, for
the $x$-subproblem, the function $Q(x,y)$ is linearized at the point $x_{k}$,
and for the $y$-subproblem, the function $Q(x,y)$ is linearized at the point
$y_{k}$ that yields the following algorithm:

$\begin{cases}x_{k+1}\in\arg\min_{x\in\mathbb{R}^{l}}\\{f(x)+\langle
x,\nabla_{x}Q(x_{k},y_{k})\rangle+\frac{1}{2\lambda_{k}}\|x-x_{k}\|^{2}_{2}\\},\\\
y_{k+1}\in\arg\min_{y\in\mathbb{R}^{m}}\\{g(y)+\langle
y,\nabla_{y}Q(x_{k+1},y_{k})\rangle+\frac{1}{2\mu_{k}}\|y-y_{k}\|^{2}_{2}\\}.\\\
\end{cases}$ (1.4)

In this way, if proximal operator of $f$ or $g$ have a closed-form or can be
easily calculated, then the problem (1.4) can be easily solved. They proved
that each bounded sequence generated by PALM globally converges to a critical
point. When $f$ and $g$ are continuously differentiable, a natural idea is to
linearize $f$ and $g$. Corresponding algorithm was proposed by Nikolova et al.
[14], called alternating structure-adapted proximal gradient descent (ASAP)
algorithm. With the help of Kurdyka–Łojasiewicz property they establish the
convergence of the whole sequence generated by proposed algorithm.

The inertial extrapolation technique has been widely used to accelerate the
iterative algorithms for convex and nonconvex optimizations, since the cost of
each iteration stays basically unchanged [15, 16]. The inertial scheme,
starting from the so-called heavy ball method of Polyak [17], was recently
proved to be very efficient in accelerating numerical methods, especially the
first-order methods. The main feature of the idea is that the new iteration
uses the previous two or more iterations.

Recently, there are increasing interests in studying inertial type algorithms,
such as Nicolas et al. [18] proposed inertial versions of block coordinate
descent methods for solving nonconvex nonsmooth optimization problems. Feng et
al. [19] focused on a minimization optimization model that is nonconvex and
nonsmooth and established an inertial Douglas–Rachford splitting (IDRS)
algorithm, which incorporate the inertial technique into the framework of the
Douglas–Rachford splitting algorithm. Specially, for solving problem (1.1),
Zhang and He [20] introduced an inertial version of the proximal alternating
minimization method, Pock and Sabach [21] proposed the following inertial
proximal alternating linearized minimization (iPALM) algorithm:

$\begin{cases}u_{1k}=x_{k}+\alpha_{1k}(x_{k}-x_{k-1}),v_{1k}=x_{k}+\beta_{1k}(x_{k}-x_{k-1}),\\\
x_{k+1}\in\arg\min_{x\in\mathbb{R}^{l}}\\{f(x)+\langle
x,\nabla_{x}Q(v_{1k},y_{k})\rangle+\frac{1}{2\lambda_{k}}\|x-u_{1k}\|^{2}_{2}\\},\\\
u_{2k}=y_{k}+\alpha_{2k}(y_{k}-y_{k-1}),v_{2k}=y_{k}+\beta_{2k}(y_{k}-y_{k-1}),\\\
y_{k+1}\in\arg\min_{y\in\mathbb{R}^{m}}\\{g(y)+\langle
y,\nabla_{y}Q(x_{k+1},v_{2k})\rangle+\frac{1}{2\mu_{k}}\|y-u_{2k}\|^{2}_{2}\\},\\\
\end{cases}$ (1.5)

where $\alpha_{1k},\alpha_{2k},\beta_{1k},\beta_{2k}\in\left[0,1\right]$. They
proved that the generated sequence globally converges to critical point of the
objective function under the condition of the Kurdyka–Łojasiewicz property.
When $\alpha_{1k}=\alpha_{2k}=\beta_{1k}=\beta_{2k}=0$, iPALM reduces to PALM.
Then Cai et al. [22] presented a Gauss–Seidel type inertial proximal
alternating linearized minimization (GiPALM) algorithm for solving problem
(1.1):

$\begin{cases}x_{k+1}\in\arg\min_{x\in\mathbb{R}^{n}}\\{f(x)+\langle
x,\nabla_{x}Q(\tilde{x}_{k},\tilde{y}_{k})\rangle+\frac{1}{2\lambda_{k}}\|x-\tilde{x}_{k}\|^{2}_{2}\\},\\\
\tilde{x}_{k+1}=x_{k+1}+\alpha(x_{k+1}-\tilde{x}_{k}),\alpha\in[0,1),\\\
y_{k+1}\in\arg\min_{y\in\mathbb{R}^{m}}\\{g(y)+\langle
y,\nabla_{y}Q(\tilde{x}_{k+1},\tilde{y}_{k})\rangle+\frac{1}{2\mu_{k}}\|y-\tilde{y}_{k}\|^{2}_{2}\\},\\\
\tilde{y}_{k+1}=y_{k+1}+\beta(y_{k+1}-\tilde{y}_{k}),\beta\in[0,1).\end{cases}$
(1.6)

In order to use the existing information as much as possible to further
improve the numerical performance, Han et al. [23] proposed a new inertial
version of proximal alternating linearized minimization (NiPALM) algorithm,
which inherits both advantages of iPALM and GiPALM. And Xu et al. [24]
proposed accelerated alternating structure-adapted proximal gradient descent
algorithm by using inertial extrapolation technique.

Bregman distance is a useful substitute for a distance, obtained from the
various choices of functions. The applications of Bregman distance instead of
the norm give us alternative ways for more flexibility in the selection of
regularization. Choosing appropriate Bregeman distances can obtain colsed form
of solution for solving subproblem. Bregman distance regularization is also an
effective way to improve the numerical results of the algorithm. In [25], the
authors constructed the following two-step inertial Bregman alternating
minimization (TiBAM) algorithm using the information of the previous three
iterates:

$\begin{cases}x_{k+1}\in\arg\min_{x\in\mathbb{R}^{l}}\\{L(x,y_{k})+D_{\phi_{1}}(x,x_{k})+\alpha_{1k}\langle
x,x_{k-1}-x_{k}\rangle+\alpha_{2k}\langle x,x_{k-2}-x_{k-1}\rangle\\},\\\
y_{k+1}\in\arg\min_{y\in\mathbb{R}^{m}}\\{L(x_{k+1},y)+D_{\phi_{2}}(y,y_{k})+\beta_{1k}\langle
y,y_{k-1}-y_{k}\rangle+\beta_{2k}\langle y,y_{k-2}-y_{k-1}\\},\end{cases}$
(1.7)

where $D_{\phi_{i}}(i=1,2)$ denote the Bregman distance with respect to
$\phi_{i}(i=1,2)$, respectively. The convergence is obtained provided an
appropriate regularization of the objective function satisfies the
Kurdyka–Łojasiewicz inequality. Based on alternating minimization algorithm,
Chao et al. [26] proposed the following inertial alternating minimization with
Bregman distance (BIAM) algorithm:

$\begin{cases}x_{k+1}\in\arg\min_{x\in\mathbb{R}^{l}}\\{f(x)+Q(x,\hat{y}_{k})+\lambda_{k}D_{\phi_{1}}(x,\hat{x}_{k})\\},\\\
\hat{x}_{k+1}=x_{k+1}+\alpha(x_{k+1}-\hat{x}_{k}),\alpha\in[0,1),\\\
y_{k+1}\in\arg\min_{y\in\mathbb{R}^{m}}\\{g(y)+Q(\hat{x}_{k+1},y)+\mu_{k}D_{\phi_{2}}(y,\hat{y}_{k})\\},\\\
\hat{y}_{k+1}=y_{k+1}+\beta(y_{k+1}-\hat{y}_{k}),\beta\in[0,1).\end{cases}$
(1.8)

Suppose that the benefit function satisfies the Kurdyka–Łojasiewicz property
and the parameters are selected appropriately, they proved the convergence of
BIAM algorithm.

In this paper, based on the proximal alternating minimization algorithm, we
combine inertial technique, linearization technique and Bregman distance to
construct two-step inertial Bregman proximal alternating linearized
minimization algorithm using the information from the previous three
iterations, and give a global convergence analysis of the algorithm under the
assumption that the benefit function satisfies the Kurdyka–Łojasiewicz
property.

The article is organized as follows: In Section 2, we recall some concepts and
important lemmas which will be used in the proof of main results. In Section
3, we present the two-step inertial Bregman proximal alternating linearized
minimization algorithm and show its convergence. Finally, in Section 4, the
preliminary numerical examples on sparse nonnegative matrix factorization,
signal recovery and quadratic fractional programming problem are provided to
illustrate the behavior of the proposed algorithm.

##  2 Preliminaries

###  2.1 Subdifferentials of nonconvex and nonsmooth functions

Consider the Euclidean vector space $\mathbb{R}^{d}$ of dimension $d\geq 1$
and denote the standard inner product by $\langle\cdot,\cdot\rangle$ and the
induced norm by $\|\cdot\|_{2}$. We use $\omega(x_{k})=\\{x:\exists
x_{k_{j}}\rightarrow x\\}$ to stand for the limit set of
$\\{x_{k}\\}_{k\in\mathbb{N}}$.

Given a function $f:\mathbb{R}^{d}\rightarrow(-\infty,+\infty]$, the domain
and the epigraph of $f$ are defined by
dom$f:=\\{x\in\mathbb{R}^{d}:f(x)<+\infty\\}$ and
epi$f:=\\{(x,\alpha)\in\mathbb{R}^{d}\times\mathbb{R}:f(x)\leq\alpha\\}$,
respectively. We say that $f$ is proper if dom$f\neq\emptyset$, $f$ is lower
semicontinuous if epi$f$ is closed. Further we recall some generalized
subdifferential notions and the basic properties which are needed in the
paper, see [27, 28].

######  Definition 2.1.

(Subdifferentials) Let $f:\mathbb{R}^{d}\rightarrow(-\infty,+\infty]$ be a
proper and lower semicontinuous function. For $x\in$ dom$f$, the Fréchet
subdifferential of $f$ at $x$, written $\hat{\partial}f(x)$, is the set of
vectors $v\in\mathbb{R}^{d}$ which satisfy

$\liminf_{y\rightarrow x}\frac{1}{\|x-y\|_{2}}[f(y)-f(x)-\langle
v,y-x\rangle]\geq 0.$

If $x\not\in$ dom$f$, let $\hat{\partial}f(x)=\emptyset$. The limiting-
subdifferential [27], or simply the subdifferential for short, of $f$ at
$x\in$ dom$f$, written $\partial f(x)$, is defined as follows:

$\partial f(x):=\\{v\in\mathbb{R}^{d}:\exists x_{k}\rightarrow
x,f(x_{k})\rightarrow f(x),v_{k}\in\hat{\partial}f(x_{k}),v_{k}\rightarrow
v\\}.$

######  Remark 2.1.

(a) The above definition implies that $\hat{\partial}f(x)\subseteq\partial
f(x)$ for each $x\in\mathbb{R}^{d}$, where the first set is convex and closed
while the second one is closed (see[29]).

(b) (Closedness of $\partial f$) Let $\\{x_{k}\\}_{k\in\mathbb{N}}$ and
$\\{v_{k}\\}_{k\in\mathbb{N}}$ be sequences in $\mathbb{R}^{d}$ such that
$v_{k}\in\partial f(x_{k})$ for all $k\in\mathbb{N}$. If
$(x_{k},v_{k})\rightarrow(x,v)$ and $f(x_{k})\rightarrow f(x)$ as
$k\rightarrow\infty$, then $v\in\partial f(x)$.

(c) A necessary (but not sufficient) condition for $x\in\mathbb{R}^{d}$ to be
a minimizer of $f$ is

$0\in\partial f(x).$ $None$

A point satisfying (2.1) is called limiting-critical or simply critical. The
set of critical points of $f$ is denoted by crit$f$.

(d) If $f:\mathbb{R}^{d}\rightarrow(-\infty,+\infty]$ be a proper and lower
semicontinuous and $h:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is a continuously
differentiable function, then $\partial(f+h)(x)=\partial f(x)+\nabla h(x)$ for
all $x\in\mathbb{R}^{d}$.

In what follows, we will consider the problem of finding a critical point
$(x^{\ast},y^{\ast})\in$dom$L$.

###  2.2 The Kurdyka–Łojasiewicz property

Now we state the definition of the Kurdyka–Łojasiewicz property. This class of
functions will play a crucial role in the convergence results of the proposed
algorithm.

Let $f:\mathbb{R}^{d}\rightarrow(-\infty,+\infty]$ be a proper and lower
semicontinuous function. For $\eta_{1}$, $\eta_{2}$ such that
$-\infty<\eta_{1}<\eta_{2}\leq+\infty$, we set

$[\eta_{1}<f<\eta_{2}]=\\{x\in\mathbb{R}^{d}:\eta_{1}<f(x)<\eta_{2}\\}.$

######  Definition 2.2.

(Kurdyka–Łojasiewicz property) The function
$f:\mathbb{R}^{d}\rightarrow(-\infty,+\infty]$ is said to have the
Kurdyka–Łojasiewicz (KL) property at $x^{\ast}\in$dom$f$ if there exist
$\eta\in(0,+\infty]$, a neighborhood $U$ of $x^{\ast}$ and a continuous
concave function $\varphi:[0,\eta)\rightarrow\mathbb{R}_{+}$ such that
$\varphi(0)=0$, $\varphi$ is $C^{1}$ on $(0,\eta)$, for all $s\in(0,\eta)$ it
is $\varphi^{\prime}(s)>0$ and for all $x$ in
$U\cap[f(x^{\ast})<f<f(x^{\ast})+\eta]$ the Kurdyka–Łojasiewicz inequality
holds,

$\varphi^{\prime}(f(x)-f(x^{\ast})){\rm dist}(0,\partial f(x))\geq 1.$

Proper lower semicontinuous functions which satisfy the Kurdyka–Łojasiewicz
inequality at each point of its domain are called KL functions.

A remarkable aspect of KL functions is that they are ubiquitous in
applications, typical functions that have the KL property include strongly
convex functions, real analytic functions, semi-algebraic functions [13],
subanalytic functions [30] and log-exp functions. We note that the $L_{q}$
quasi-norm $\|x\|q:=(\Sigma_{i}|x_{i}|^{q})^{1/q}$ with $0<q<1$, the sup-norm
$\|x\|_{\infty}:=\max_{i}|x_{i}|$, $\|Ax-b\|_{q}^{q}$ with $0<q<1$ and
$\|Ax-b\|_{\infty}$ are all semi-algebraic functions for any matrix $A$, and
so they have the KL property.

It is known that both real analytic and semi-algebraic functions are
subanalytic. Generally speaking, the sum of two subanalytic functions is not
necessarily subanalytic. However, for two subanalytic functions, if at least
one function maps bounded sets to bounded sets, then their sum is also
subanalytic, as shown in [31]. In particular, the sum of a subanalytic
function and an analytic function is subanalytic. Typical subanalytic
functions include: $\|Ax-b\|^{2}_{2}+\lambda\|y\|_{q}^{q}$,
$\|Ax-b\|_{2}^{2}+\lambda\Sigma_{i}(y_{i}^{2}+\varepsilon)^{q/2}$ and so on.

###  2.3 Bregman distance

######  Definition 2.3.

A function $f$ is said convex if dom$f$ is a convex set and if, for all $x$,
$y\in$dom$f$, $\alpha\in[0,1]$,

$f(\alpha x+(1-\alpha)y)\leq\alpha f(x)+(1-\alpha)f(y).$

$f$ is said $\theta$-strongly convex with $\theta>0$ if
$f-\frac{\theta}{2}\|\cdot\|^{2}$ is convex, i.e.,

$f(\alpha x+(1-\alpha)y)\leq\alpha
f(x)+(1-\alpha)f(y)-\frac{1}{2}\theta\alpha(1-\alpha)\|x-y\|^{2}$

for all $x$, $y\in$dom$f$ and $\alpha\in[0,1]$.

Suppose that the function $f$ is differentiable. Then $f$ is convex if and
only if dom$f$ is a convex set and

$f(x)\geq f(y)+\langle\nabla f(y),x-y\rangle$

holds for all $x$, $y\in$dom$f$. Moreover, $f$ is $\theta$-strongly convex
with $\theta>0$ if and only if

$f(x)\geq f(y)+\langle\nabla f(y),x-y\rangle+\frac{\theta}{2}\|x-y\|^{2}$

for all $x$, $y\in$dom$f$.

The Bregman distance is briefly described below, as a generalization of the
squared Euclidean distance, the Bregman distance shares many similar nice
properties of the Euclidean distance. However, the Bregman distance is not a
real metric, since it does not satisfy the triangle inequality nor symmetry.
But the use of the Bregman distance makes the algorithm more widely
applicable, so the Bregman distance plays an important role in various
iterative algorithms.

######  Definition 2.4.

Let $\phi:\mathbb{R}^{d}\rightarrow(-\infty,+\infty]$ be a convex and Gâteaux
differentiable function. The function $D_{\phi}:$ dom$\phi\,\,\times$
intdom$\phi\rightarrow[0,+\infty)$, defined by

$D_{\phi}(x,y)=\phi(x)-\phi(y)-\langle\nabla\phi(y),x-y\rangle,$

is called the Bregman distance with respect to $\phi$.

From the above definition, it follows that

$D_{\phi}(x,y)\geq\frac{\theta}{2}\|x-y\|^{2},$ $None$

if $\phi$ is $\theta$-strongly convex. Furthermore, it is easy to verify that
$D_{\phi}(x,y)\geq 0$ and $D_{\phi}(x,y)\geq D_{\phi}(z,y)$ if $z\in[y,x]$. We
take the non-Euclidean distance [32] as an example. For a positive definite
matrix $M$, let $\|x\|_{M}:=\sqrt{\langle x,Mx\rangle}$ and take
$\phi(x)=\|x\|_{M}^{2}$. Then $D_{\phi}(x,y)=\|x-y\|_{M}^{2}$. In particular,
if $M$ is an identity, $D_{\phi}(x,y)$ reduces to the Euclidean squared
distance.

###  2.4 Proximal operator for nonconvex functions

######  Definition 2.5.

(Well-definedness of proximal operator) Let
$\sigma:\mathbb{R}^{d}\rightarrow(-\infty,+\infty]$ be a proper and lower
semicontinuous function with $\inf_{\mathbb{R}^{d}}\sigma>-\infty$. Given
$x\in\mathbb{R}^{d}$ and $t>0$, the proximal operator associated to $\sigma$
is defined respectively by

$prox_{t\sigma}\left(x\right):=\arg\min_{u\in\mathbb{R}^{d}}\left\\{\sigma\left(u\right)+\frac{1}{2t}\left\|u-x\right\|^{2}_{2}\right\\}.$
$None$

Then, for every $t\in\left(0,\infty\right)$ the set
$prox_{\frac{1}{t}\sigma}\left(x\right)$ is nonempty and compact.

It is well-known that the proximal forward-backward scheme for minimizing the
sum of a smooth function $h$ with a nonsmooth one $\sigma$ can simply be
viewed as the proximal regularization of $h$ linearized at a given point
$x^{k}$:

$x_{k+1}\in\arg\min_{x\in\mathbb{R}^{l}}\\{\sigma(x)+\langle x-x_{k},\nabla
h(x_{k})\rangle+\frac{t}{2}\|x-x_{k}\|^{2}_{2}\\},\ (t>0),$ $None$

that is, using the proximal operator notation defined in (2.3), we get

$x_{k+1}\in prox_{\frac{1}{t}\sigma}\left(x_{k}-\frac{1}{t}\nabla
h(x_{k})\right).$

###  2.5 Some basic convergence properties

In the following we give some basic properties which are needed in the
convergence analysis. We begin by recalling the well-known and important
descent lemma, see [9, 27].

######  Lemma 2.1.

(Descent lemma[9]) Let $h:\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a
continuously differentiable function with gradient $\nabla h$ assumed
$L_{h}$-Lipschitz continuous. Then

$\left|h(u)-h(v)-\left\langle u-v,\nabla
h(v)\right\rangle\right|\leq\frac{L_{h}}{2}\left\|u-v\right\|^{2}_{2},\
\forall u,v\in\mathbb{R}^{d}.$ $None$

######  Assumption 2.1.

Let
$\\{u_{k}\\}_{k\in\mathbb{N}}:=\\{(z_{k},z_{k-1},z_{k-2})\\}_{k\in\mathbb{N}}$
be a sequence in $\mathbb{R}^{3p}$, where $z_{k},z_{k-1},\\\
z_{k-2}\in\mathbb{R}^{p}$. For convenience we use the abbreviation
$\triangle_{k}:=\|z_{k}-z_{k-1}\|_{2}$ for $k\in\mathbb{N}$. Let
$F:\mathbb{R}^{3p}\rightarrow{(-\infty,+\infty]}$ be a proper lower
semicontinuous function. Then, letting $a$ and $b$ are fixed positive
constants, we consider a sequence $\\{u_{k}\\}_{k\in\mathbb{N}}$ satisfying
the following conditions.

(H1) (Sufficient decrease condition). For each $k\in\mathbb{N}$, it holds that

$F(u_{k+1})+a\triangle_{k+1}^{2}\leq F(u_{k}),$

i.e., $F(z_{k+1},z_{k},z_{k-1})+a\|z_{k+1}-z_{k}\|_{2}^{2}\leq
F(z_{k},z_{k-1},z_{k-2})$ for $u_{k}=(z_{k},z_{k-1},z_{k-2})$;

(H2) (Relative error condition). For each $k\in\mathbb{N}$, there exists
$w_{k+1}\in\partial F(u_{k+1})$ such that

$\|w_{k+1}\|_{2}\leq b(\triangle_{k+1}+\triangle_{k}+\triangle_{k-1}),$

i.e., ${\|w_{k+1}\|_{2}\leq
b(\|z_{k+1}-z_{k}\|_{2}+\|z_{k}-z_{k-1}\|_{2}+\|z_{k-1}-z_{k-2}\|_{2})}$;

(H3) (Continuity condition). There exists a subsequence
$\\{u_{k_{j}}\\}_{j\in\mathbb{N}}$ and $\tilde{u}$ such that

$u_{k_{j}}\rightarrow\tilde{u}\ \ {\rm and}\ \ F(u_{k_{j}})\rightarrow
F(\tilde{u})$

as $j\rightarrow\infty$.

An abstract convergence result for inexact descent methods (see[25]) is given
below, which can be applied to the two-step inertial algorithm. So it will be
used to analyze the convergence of the proposed algorithm.

######  Lemma 2.2.

(convergence to a critical point[25]) Let
$F:\mathbb{R}^{3p}\rightarrow{(-\infty,+\infty]}$ be a proper lower
semicontinuous function and
$\\{u_{k}\\}_{k\in\mathbb{N}}=\\{(z_{k},z_{k-1},z_{k-2})\\}_{k\in\mathbb{N}}$
be a sequence which satisfies conditions (H1), (H2) and (H3). Suppose $F$ have
the KL property at the cluster point $\tilde{u}$ specified in (H3). Then the
sequence $\\{z_{k}\\}_{k\in\mathbb{N}}$ has a finite length, i.e.,
$\sum_{k=0}^{\infty}\|z_{k+1}-z_{k}\|_{2}<+\infty$. Moreover the sequence
$\\{z_{k}\\}_{k\in\mathbb{N}}$ converges to $\bar{z}=\tilde{z}$ as
$k\rightarrow\infty$, $\bar{u}=(\bar{z},\bar{z},\bar{z})$ is a critical point
of $F$ where $\tilde{u}=(\tilde{z},\tilde{z},\tilde{z})$.

##  3 Two-step inertial Bregman proximal alternating linearized minimization
algorithm

In this section we investigate the convergence property of the two-step
inertial Bregman proximal alternating linearized minimization algorithm for
solving nonconvex and nonsmooth nonseparable optimization problems. We assume
that $\phi_{i}$ is $\theta_{i}$-strongly convex and $\nabla\phi_{i}$ is
$L_{\nabla\phi_{i}}$-Lipschitz continuous for $i=1,2$. Set
$\theta=\min\\{\theta_{1},\theta_{2}\\}$.

######  Algorithm 3.1.

Choose $(x_{0},y_{0})\in{\text{d}om}L$ and set
$(x_{-i},y_{-i})=(x_{0},y_{0})$, $i=1,2$. Take the sequences
$\\{\alpha_{1k}\\}$, $\\{\beta_{1k}\\}\subseteq[0,\alpha_{1}]$ and
$\\{\alpha_{2k}\\}$, $\\{\beta_{2k}\\}\subseteq[0,\alpha_{2}]$, where
$\alpha_{1}\geq 0$ and $\alpha_{2}\geq 0$. For $k\geq 0$, let

$\begin{cases}\begin{aligned}
x_{k+1}\in\arg\min_{x\in\mathbb{R}^{l}}\\{&f(x)+\langle
x,\nabla_{x}Q(x_{k},y_{k})\rangle+D_{\phi_{1}}(x,x_{k})+\alpha_{1k}\langle
x,x_{k-1}-x_{k}\rangle\\\ &+\alpha_{2k}\langle x,x_{k-2}-x_{k-1}\rangle\\},\\\
y_{k+1}\in\arg\min_{y\in\mathbb{R}^{m}}\\{&g(y)+\langle
y,\nabla_{y}Q(x_{k+1},y_{k})\rangle+D_{\phi_{2}}(y,y_{k})+\beta_{1k}\langle
y,y_{k-1}-y_{k}\rangle\\\ &+\beta_{2k}\langle
y,y_{k-2}-y_{k-1}\rangle\\},\end{aligned}\end{cases}$ $None$

where $D_{\phi_{1}}$ and $D_{\phi_{2}}$ denote the Bregman distance with
respect to $\phi_{1}$ and $\phi_{2}$, respectively.

Since $f$ and $g$ are proper, lower semicontinuous, $f$ and $g$ are bounded
from below, and $\phi_{i}$ ($i=1,2$) is strongly convex, we get that it makes
sense to define the iterative scheme, meaning that the existence of
$(x_{k+1},y_{k+1})$ is guaranteed for each $k\geq 0$.

######  Remark 3.1.

We give special cases of Algorithm 3.1.

  * (i)

Let $\lambda$ and $\mu$ are positive constants. If we take
$\phi_{1}(x)=\frac{1}{2\lambda}\|x\|^{2}_{2}$ and
$\phi_{2}(y)=\frac{1}{2\mu}\|y\|^{2}_{2}$ for all $x\in\mathbb{R}^{l}$ and
$y\in\mathbb{R}^{m}$, then the iterative Algorithm 3.1 becomes

$\begin{cases}\begin{aligned}
x_{k+1}\in\arg\min_{x\in\mathbb{R}^{l}}\\{&f(x)+\langle
x,\nabla_{x}Q(x_{k},y_{k})\rangle+\frac{1}{2\lambda}\|x-x_{k}\|^{2}_{2}+\alpha_{1k}\langle
x,x_{k-1}-x_{k}\rangle\\\ &+\alpha_{2k}\langle x,x_{k-2}-x_{k-1}\rangle\\},\\\
y_{k+1}\in\arg\min_{y\in\mathbb{R}^{m}}\\{&g(y)+\langle
y,\nabla_{y}Q(x_{k+1},y_{k})\rangle+\frac{1}{2\mu}\|y-y_{k}\|^{2}_{2}+\beta_{1k}\langle
y,y_{k-1}-y_{k}\rangle\\\ &+\beta_{2k}\langle
y,y_{k-2}-y_{k-1}\rangle\\}.\end{aligned}\end{cases}$ $None$

In this case, letting $\alpha_{2k}\equiv\beta_{2k}\equiv 0$ for all $k\geq 0$,
(3.2) becomes the one-step inertial algorithm

$\begin{cases}\begin{aligned}
x_{k+1}\in\arg\min_{x\in\mathbb{R}^{l}}\\{&f(x)+\langle
x,\nabla_{x}Q(x_{k},y_{k})\rangle+\frac{1}{2\lambda}\|x-x_{k}\|^{2}_{2}+\alpha_{1k}\langle
x,x_{k-1}-x_{k}\rangle\\},\\\
y_{k+1}\in\arg\min_{y\in\mathbb{R}^{m}}\\{&g(y)+\langle
y,\nabla_{y}Q(x_{k+1},y_{k})\rangle+\frac{1}{2\mu}\|y-y_{k}\|^{2}_{2}+\beta_{1k}\langle
y,y_{k-1}-y_{k}\rangle\\}.\end{aligned}\end{cases}$

Letting $\alpha_{2k}\equiv\beta_{2k}=\alpha_{1k}\equiv\beta_{1k}\equiv 0$ for
all $k\geq 0$, (3.2) becomes the proximal alternating linearized minimization
algorithm (1.4) with $\lambda_{k}\equiv\lambda$ and $\mu_{k}\equiv\mu$ for
$k\geq 0$.

  * (ii)

Letting $\alpha_{2k}\equiv\beta_{2k}\equiv 0$ for all $k\geq 0$, (3.1) becomes
the one-step inertial Bregman proximal alternating linearized minimization
algorithm (iBPALM):

$\begin{cases}\begin{aligned}
x_{k+1}\in\arg\min_{x\in\mathbb{R}^{l}}\\{&f(x)+\langle
x,\nabla_{x}Q(x_{k},y_{k})\rangle+D_{\phi_{1}}(x,x_{k})+\alpha_{1k}\langle
x,x_{k-1}-x_{k}\rangle\\},\\\
y_{k+1}\in\arg\min_{y\in\mathbb{R}^{m}}\\{&g(y)+\langle
y,\nabla_{y}Q(x_{k+1},y_{k})\rangle+D_{\phi_{2}}(y,y_{k})+\beta_{1k}\langle
y,y_{k-1}-y_{k}\rangle\\}.\end{aligned}\end{cases}$

  * (iii)

Letting $\alpha_{2k}\equiv\beta_{2k}=\alpha_{1k}\equiv\beta_{1k}\equiv 0$ for
all $k\geq 0$, (3.1) becomes the Bregman proximal alternating minimization
algorithm (BPALM):

$\begin{cases}\begin{aligned}
x_{k+1}\in\arg\min_{x\in\mathbb{R}^{l}}\\{&f(x)+\langle
x,\nabla_{x}Q(x_{k},y_{k})\rangle+D_{\phi_{1}}(x,x_{k})\\},\\\
y_{k+1}\in\arg\min_{y\in\mathbb{R}^{m}}\\{&g(y)+\langle
y,\nabla_{y}Q(x_{k+1},y_{k})\rangle+D_{\phi_{2}}(y,y_{k})\\}.\end{aligned}\end{cases}$

In order to get convergence analysis of our algorithm, the following
assumptions required are given.

######  Assumption 3.1.

(i) For any fixed $y$, the partial gradient $\nabla_{x}Q$ is globally
Lipschitz with module $L_{1}{\left(y\right)}$, that is,

$\left\|\nabla_{x}Q\left(x_{1},y\right)-\nabla_{x}Q\left(x_{2},y\right)\right\|\leq
L_{1}{\left(y\right)}\left\|x_{1}-x_{2}\right\|,\ \forall
x_{1},x_{2}\in\mathbb{R}^{l}.$

Likewise, for any fixed $x$, the partial gradient $\nabla_{y}Q$ is globally
Lipschitz with module $L_{2}{\left(x\right)}$,

$\left\|\nabla_{y}Q\left(x,y_{1}\right)-\nabla_{y}Q\left(x,y_{2}\right)\right\|\leq
L_{2}{\left(x\right)}\left\|y_{1}-y_{2}\right\|,\ \forall
y_{1},y_{2}\in\mathbb{R}^{m}.$

(ii) For $i=1,2$, there exists $L_{i}^{+}>0$ such that

$\inf\left\\{L_{1}\left(y_{k}\right):k\in N\right\\}\geq L_{1}^{-}\ \ {\rm
and}\ \inf\left\\{L_{2}\left(x_{k}\right):k\in N\right\\}\geq L_{2}^{-},$
$None$

And $\theta_{1}>L_{1}^{+}$, $\theta_{2}>L_{2}^{+}$ hold, let
$\rho=\min\\{\theta_{1}-L_{1}^{+}\ ,\theta_{2}-L_{2}^{+}\\}$.

In all of what follows, let $\\{(x_{k},y_{k})\\}_{k\in\mathbb{N}}$ be the
sequence generated by Algorithm 3.1. Furthermore, for more convenient notation
we abbreviate $z_{k}=(x_{k},y_{k})$, $\triangle_{k}=\|z_{k}-z_{k-1}\|_{2}$ for
$k\geq 0$, and introduce an benefit function

$H(u,v,w):=L(u)+\frac{\alpha_{1}+\alpha_{2}}{2}\|u-v\|_{2}^{2}+\frac{\alpha_{2}}{2}\|v-w\|_{2}^{2}.$
$None$

Note that $H(u,v,w)=L(u)$ for $u=v=w$.

######  Lemma 3.1.

Let $2(\alpha_{1}+\alpha_{2})<\rho$ and the sequence
$\\{z_{k}\\}_{k\in\mathbb{N}}$ be generated by Algorithm 3.1. Then

(i) the sequence $\\{H(z_{k},z_{k-1},z_{k-2})\\}_{k\in\mathbb{N}}$ is
monotonically nonincreasing and convergent, and it holds that

$H(z_{k+1},z_{k},z_{k-1})+a\triangle_{k+1}^{2}\leq H(z_{k},z_{k-1},z_{k-2}),\
k\geq 0,$ $None$

where $a=\frac{\rho-2(\alpha_{1}+\alpha_{2})}{2}>0$;

(ii) it holds that $\sum_{k=0}^{\infty}\triangle_{k}^{2}<+\infty$, and thus
$\lim_{k\rightarrow\infty}\triangle_{k}=0$;

(iii) the sequence $\\{L(x_{k},y_{k})\\}_{k\in\mathbb{N}}$ is convergent.

###### Proof.

(i) From Algorithm 3.1, we have

$\displaystyle f(x_{k+1})+\langle
x_{k+1},\nabla_{x}Q(x_{k},y_{k})\rangle+D_{\phi_{1}}(x_{k+1},x_{k})+\alpha_{1k}\langle
x_{k+1},x_{k-1}-x_{k}\rangle$ $None$ $\displaystyle+\alpha_{2k}\langle
x_{k+1},x_{k-2}-x_{k-1}\rangle$ $\displaystyle\leq$ $\displaystyle
f(x_{k})+\langle
x_{k},\nabla_{x}Q(x_{k},y_{k})\rangle+D_{\phi_{1}}(x_{k},x_{k})+\alpha_{1k}\langle
x_{k},x_{k-1}-x_{k}\rangle$ $\displaystyle+\alpha_{2k}\langle
x_{k},x_{k-2}-x_{k-1}\rangle$ $\displaystyle=$ $\displaystyle f(x_{k})+\langle
x_{k},\nabla_{x}Q(x_{k},y_{k})\rangle+\alpha_{1k}\langle
x_{k},x_{k-1}-x_{k}\rangle+\alpha_{2k}\langle x_{k},x_{k-2}-x_{k-1}\rangle.$

From Lemma 2.1, we have

$Q\left(x_{k+1},y_{k}\right)\leq Q\left(x_{k},y_{k}\right)+\left\langle
x_{k+1}-x_{k},\nabla_{x}Q\left(x_{k},y_{k}\right)\right\rangle+\frac{L_{1}(y_{k})}{2}{\left\|x_{k+1}-x_{k}\right\|^{2}_{2}}.$
$None$

By (3.5) and (3.6), we have

$\displaystyle
f(x_{k+1})+Q(x_{k+1},y_{k})+D_{\phi_{1}}(x_{k+1},x_{k})+\alpha_{1k}\langle
x_{k+1},x_{k-1}-x_{k}\rangle+\alpha_{2k}\langle
x_{k+1},x_{k-2}-x_{k-1}\rangle$ $None$ $\displaystyle\leq$ $\displaystyle
f(x_{k+1})+Q\left(x_{k},y_{k}\right)+\left\langle
x_{k+1}-x_{k},\nabla_{x}Q\left(x_{k},y_{k}\right)\right\rangle+\frac{L_{1}(y_{k})}{2}{\left\|x_{k+1}-x_{k}\right\|^{2}_{2}}+D_{\phi_{1}}(x_{k+1},x_{k})$
$\displaystyle+\alpha_{1k}\langle
x_{k+1},x_{k-1}-x_{k}\rangle+\alpha_{2k}\langle
x_{k+1},x_{k-2}-x_{k-1}\rangle$ $\displaystyle\leq$ $\displaystyle
f(x_{k})+Q\left(x_{k},y_{k}\right)+\frac{L_{1}(y_{k})}{2}{\left\|x_{k+1}-x_{k}\right\|^{2}_{2}}+\alpha_{1k}\langle
x_{k},x_{k-1}-x_{k}\rangle+\alpha_{2k}\langle x_{k},x_{k-2}-x_{k-1}\rangle,$

which implies that

$\displaystyle
L(x_{k+1},y_{k})+D_{\phi_{1}}(x_{k+1},x_{k})-\frac{L_{1}(y_{k})}{2}{\left\|x_{k+1}-x_{k}\right\|^{2}_{2}}$
$None$ $\displaystyle\leq$ $\displaystyle L(x_{k},y_{k})+\alpha_{1k}\langle
x_{k+1}-x_{k},x_{k}-x_{k-1}\rangle+\alpha_{2k}\langle
x_{k+1}-x_{k},x_{k-1}-x_{k-2}\rangle.$

By (2.2) and Assumption 3.1 (ii), it follows from (3.8) that

$\displaystyle
L(x_{k+1},y_{k})+\frac{1}{2}\left(\theta_{1}-L_{1}^{+}\right)\|x_{k+1}-x_{k}\|_{2}^{2}$
$None$ $\displaystyle\leq$ $\displaystyle L(x_{k},y_{k})+\alpha_{1k}\langle
x_{k+1}-x_{k},x_{k}-x_{k-1}\rangle+\alpha_{2k}\langle
x_{k+1}-x_{k},x_{k-1}-x_{k-2}\rangle.$

Similarly, we have

$\displaystyle
L(x_{k+1},y_{k+1})+\frac{1}{2}\left(\theta_{2}-L_{2}^{+}\right)\|y_{k+1}-y_{k}\|_{2}^{2}$
$None$ $\displaystyle\leq$ $\displaystyle L(x_{k+1},y_{k})+\beta_{1k}\langle
y_{k+1}-y_{k},y_{k}-y_{k-1}\rangle+\beta_{2k}\langle
y_{k+1}-y_{k},y_{k-1}-y_{k-2}\rangle.$

Adding up the inequalities (3.9) and (3.10), we get

$\displaystyle
L(x_{k+1},y_{k+1})+\frac{1}{2}\left(\theta_{1}-L_{1}^{+}\right)\|x_{k+1}-x_{k}\|_{2}^{2}+\frac{1}{2}\left(\theta_{2}-L_{2}^{+}\right)\|y_{k+1}-y_{k}\|_{2}^{2}$
$None$ $\displaystyle\leq$ $\displaystyle L(x_{k},y_{k})+\alpha_{1k}\langle
x_{k+1}-x_{k},x_{k}-x_{k-1}\rangle+\alpha_{2k}\langle
x_{k+1}-x_{k},x_{k-1}-x_{k-2}\rangle$ $\displaystyle+\beta_{1k}\langle
y_{k+1}-y_{k},y_{k}-y_{k-1}\rangle+\beta_{2k}\langle
y_{k+1}-y_{k},y_{k-1}-y_{k-2}\rangle$ $\displaystyle\leq$ $\displaystyle
L(x_{k},y_{k})+\frac{\alpha_{1k}}{2}(\|x_{k+1}-x_{k}\|_{2}^{2}+\|x_{k}-x_{k-1}\|_{2}^{2})+\frac{\alpha_{2k}}{2}(\|x_{k+1}-x_{k}\|_{2}^{2}+\|x_{k-1}-x_{k-2}\|_{2}^{2})$
$\displaystyle+\frac{\beta_{1k}}{2}(\|y_{k+1}-y_{k}\|_{2}^{2}+\|y_{k}-y_{k-1}\|_{2}^{2})+\frac{\beta_{2k}}{2}(\|y_{k+1}-y_{k}\|_{2}^{2}+\|y_{k-1}-y_{k-2}\|_{2}^{2})$
$\displaystyle\leq$ $\displaystyle
L(x_{k},y_{k})+(\frac{\alpha_{1k}+\alpha_{2k}}{2}\|x_{k+1}-x_{k}\|_{2}^{2}+\frac{\beta_{1k}+\beta_{2k}}{2}\|y_{k+1}-y_{k}\|_{2}^{2})+(\frac{\alpha_{1k}}{2}\|x_{k}-x_{k-1}\|_{2}^{2}$
$\displaystyle+\frac{\beta_{1k}}{2}\|y_{k}-y_{k-1}\|_{2}^{2})+(\frac{\alpha_{2k}}{2}\|x_{k-1}-x_{k-2}\|_{2}^{2}+\frac{\beta_{2k}}{2}\|y_{k-1}-y_{k-2}\|_{2}^{2}).$

Since $\rho=\min\\{\theta_{1}-L_{1}^{+}\ ,\theta_{2}-L_{2}^{+}\\}$,
$\\{\alpha_{1k}\\}$, $\\{\beta_{1k}\\}\subseteq[0,\alpha_{1}]$ and
$\\{\alpha_{2k}\\}$, $\\{\beta_{2k}\\}\subseteq[0,\alpha_{2}]$, it follows
that

$\displaystyle L(z_{k+1})+\frac{\rho}{2}\|z_{k+1}-z_{k}\|_{2}^{2}$
$\displaystyle\leq$ $\displaystyle
L(z_{k})+\frac{\alpha_{1}+\alpha_{2}}{2}\|z_{k+1}-z_{k}\|_{2}^{2}+\frac{\alpha_{1}}{2}\|z_{k}-z_{k-1}\|_{2}^{2}+\frac{\alpha_{2}}{2}\|z_{k-1}-z_{k-2}\|_{2}^{2},$

which implies that

$L(z_{k+1})+\frac{\alpha_{1}+\alpha_{2}}{2}\triangle_{k+1}^{2}+\frac{\alpha_{2}}{2}\triangle_{k}^{2}+\frac{\rho-2(\alpha_{1}+\alpha_{2})}{2}\triangle_{k+1}^{2}\leq
L(z_{k})+\frac{\alpha_{1}+\alpha_{2}}{2}\triangle_{k}^{2}+\frac{\alpha_{2}}{2}\triangle_{k-1}^{2}.$

So, (3.4) holds and the sequence
$\\{H(z_{k},z_{k-1},z_{k-2})\\}_{k\in\mathbb{N}}$ is monotonically
nonincreasing. By assumption, $L$ is bounded from below, and hence
$\\{H(z_{k},z_{k-1},z_{k-2})\\}_{k\in\mathbb{N}}$ converges.

(ii) It follows from (3.4) that

$a\triangle_{k+1}^{2}\leq H(z_{k},z_{k-1},z_{k-2})-H(z_{k+1},z_{k},z_{k-1}).$
$None$

Summing up (3.12) from $k=0,\cdots,K$ yields (note that
$H(z_{0},z_{-1},z_{-2})=L(z_{0})$)

$\sum_{k=0}^{K}a\triangle_{k+1}^{2}\leq L(z_{0})-H(z_{K+1},z_{K},z_{K-1})\leq
L(z_{0})-\inf_{\mathbb{R}^{l}\times\mathbb{R}^{m}}L<\infty.$

Let $K\rightarrow+\infty$, we have that (ii) holds.

(iii) By

$L(x_{k},y_{k})=L(z_{k})=H(z_{k},z_{k-1},z_{k-2})-\frac{\alpha_{1}+\alpha_{2}}{2}\triangle_{k}^{2}-\frac{\alpha_{2}}{2}\triangle_{k-1}^{2},$

it follows from (i) and (ii) that the sequence
$\\{L(x_{k},y_{k})\\}_{k\in\mathbb{N}}$ converges and

$\lim_{k\rightarrow\infty}L(x_{k},y_{k})=\lim_{k\rightarrow\infty}H(z_{k},z_{k-1},z_{k-2}).$

∎

######  Lemma 3.2.

Let the sequence $\\{z_{k}\\}_{k\in\mathbb{N}}$ be generated by Algorithm 3.1.
For $k\geq 0$, define

$\displaystyle w_{k+1}$ $None$ $\displaystyle=$
$\displaystyle(\nabla_{x}Q(x_{k+1},y_{k+1})-\nabla_{x}Q(x_{k},y_{k})+\nabla\phi_{1}(x_{k})-\nabla\phi_{1}(x_{k+1}),\nabla_{y}Q(x_{k+1},y_{k+1})-\nabla_{y}Q(x_{k+1},y_{k})$
$\displaystyle+\nabla\phi_{2}(y_{k})-\nabla\phi_{2}(y_{k+1}),0,0,0,0)$
$\displaystyle+(\alpha_{1k}(x_{k}-x_{k-1})+\alpha_{2k}(x_{k-1}-x_{k-2}),\beta_{1k}(y_{k}-y_{k-1})+\beta_{2k}(y_{k-1}-y_{k-2}),0,0,0,0)$
$\displaystyle+((\alpha_{1}+\alpha_{2})(x_{k+1}-x_{k}),(\alpha_{1}+\alpha_{2})(y_{k+1}-y_{k}),0,0,0,0)$
$\displaystyle+(0,0,(\alpha_{1}+\alpha_{2})(x_{k}-x_{k+1})+\alpha_{2}(x_{k}-x_{k-1}),(\alpha_{1}+\alpha_{2})(y_{k}-y_{k+1})+\alpha_{2}(y_{k}-y_{k-1}),0,0)$
$\displaystyle+(0,0,0,0,\alpha_{2}(x_{k-1}-x_{k}),\alpha_{2}(y_{k-1}-y_{k})).$

Then

$w_{k+1}\in\partial
H(z_{k+1},z_{k},z_{k-1})={\partial}H(x_{k+1},y_{k+1},x_{k},y_{k},x_{k-1},y_{k-1}).$
$None$

Moreover, if the sequence $\\{z_{k}\\}_{k\in\mathbb{N}}$ is bounded, then
there exists $M>0$ such that

$\displaystyle\|w_{k+1}\|_{2}$ $None$ $\displaystyle\leq$
$\displaystyle(M+L_{\nabla\phi_{1}}+2\alpha_{1}+2\alpha_{2})\|x_{k+1}-x_{k}\|_{2}+(M+L_{2}{\left(x_{k+1}\right)}+L_{\nabla\phi_{2}}+2\alpha_{1}+2\alpha_{2})\|y_{k+1}-y_{k}\|_{2}$
$\displaystyle+(\alpha_{1}+2\alpha_{2})(\|x_{k}-x_{k-1}\|_{2}+\|y_{k}-y_{k-1}\|_{2})+\alpha_{2}(\|x_{k-1}-x_{k-2}\|_{2}+\|y_{k-1}-y_{k-2}\|_{2}),$

###### Proof.

By the definition of $x_{k+1}$, and by Remark 2.1 (c), $0$ must lie in the
subdifferential at point $x_{k+1}$ of the function

$x\longmapsto f(x)+\langle
x,\nabla_{x}Q(x_{k},y_{k})\rangle+D_{\phi_{1}}(x,x_{k})+\alpha_{1k}\langle
x,x_{k-1}-x_{k}\rangle+\alpha_{2k}\langle x,x_{k-2}-x_{k-1}\rangle.$

Since $\phi$ are differential, we have

$0\in\partial
f(x_{k+1})+\nabla_{x}Q(x_{k},y_{k})+\nabla\phi_{1}(x_{k+1})-\nabla\phi_{1}(x_{k})+\alpha_{1k}(x_{k-1}-x_{k})+\alpha_{2k}(x_{k-2}-x_{k-1}),$

which implies that

$-\nabla_{x}Q(x_{k},y_{k})+\nabla\phi_{1}(x_{k})-\nabla\phi_{1}(x_{k+1})+\alpha_{1k}(x_{k}-x_{k-1})+\alpha_{2k}(x_{k-1}-x_{k-2})\in\partial
f(x_{k+1}).$ $None$

Similarly, we have

$0\in\partial
g(y_{k+1})+\nabla_{y}Q(x_{k+1},y_{k})+\nabla\phi_{2}(y_{k+1})-\nabla\phi_{2}(y_{k})+\beta_{1k}(y_{k-1}-y_{k})+\beta_{2k}(y_{k-2}-y_{k-1}),$

and so

$-\nabla_{y}Q(x_{k+1},y_{k})+\nabla\phi_{2}(y_{k})-\nabla\phi_{2}(y_{k+1})+\beta_{1k}(y_{k}-y_{k-1})+\beta_{2k}(y_{k-1}-y_{k-2})\in\partial
g(y_{k+1}).$ $None$

Because of the structure of $L$, it follows that
$\partial_{x}L(x_{k+1},y_{k+1})=\partial
f(x_{k+1})+\nabla_{x}Q(x_{k+1},y_{k+1})$ and
$\partial_{y}L(x_{k+1},y_{k+1})=\nabla_{y}Q(x_{k+1},y_{k+1})+\partial
g(y_{k+1})$. Hence, from (3.16) and (3.17), we have

$\displaystyle\nabla_{x}Q(x_{k+1},y_{k+1})-\nabla_{x}Q(x_{k},y_{k})+\nabla\phi_{1}(x_{k})-\nabla\phi_{1}(x_{k+1})+\alpha_{1k}(x_{k}-x_{k-1})+\alpha_{2k}(x_{k-1}-x_{k-2})$
$None$ $\displaystyle\in\partial_{x}L(x_{k+1},y_{k+1}),$

and

$\displaystyle\nabla_{y}Q(x_{k+1},y_{k+1})-\nabla_{y}Q(x_{k+1},y_{k})+\nabla\phi_{2}(y_{k})-\nabla\phi_{2}(y_{k+1})+\beta_{1k}(y_{k}-y_{k-1})+\beta_{2k}(y_{k-1}-y_{k-2})$
$None$ $\displaystyle\in\partial_{y}L(x_{k+1},y_{k+1}).$

Since

$\displaystyle H(z_{k+1},z_{k},z_{k-1})$ $\displaystyle=$ $\displaystyle
L(z_{k+1})+\frac{\alpha_{1}+\alpha_{2}}{2}\|z_{k+1}-z_{k}\|_{2}^{2}+\frac{\alpha_{2}}{2}\|z_{k}-z_{k-1}\|_{2}^{2}$
$\displaystyle=$ $\displaystyle
L(x_{k+1},y_{k+1})+\frac{\alpha_{1}+\alpha_{2}}{2}(\|x_{k+1}-x_{k}\|_{2}^{2}+\|y_{k+1}-y_{k}\|_{2}^{2})+\frac{\alpha_{2}}{2}(\|x_{k}-x_{k-1}\|_{2}^{2}+\|y_{k}-y_{k-1}\|_{2}^{2}),$

we have (3.14) holds. If the sequence $\\{z_{k}\\}_{k\in\mathbb{N}}$ is
bounded, noting that $\nabla Q$ is Lipschitz continuous on bounded subsets,
then there exists $M>0$ such that

$\displaystyle\left\|\nabla_{x}Q\left(x_{k+1},y_{k+1}\right)-\nabla_{x}Q\left(x_{k},y_{k}\right))\right\|\leq\left\|\nabla
Q\left(x_{k+1},y_{k+1}\right)-\nabla Q\left(x_{k},y_{k}\right))\right\|$
$\displaystyle\leq$ $\displaystyle
M(\left\|x_{k+1}-x_{k}\right\|^{2}+\left\|y_{k+1}-y_{k}\right\|^{2})^{\frac{1}{2}}\leq
M(\left\|x_{k+1}-x_{k}\right\|+\left\|y_{k+1}-y_{k}\right\|).$

Hence, the inequality (3.15) follows from the definition of the sequence
$\\{w_{k+1}\\}_{k\in\mathbb{N}}$. ∎

######  Lemma 3.3.

Suppose that the function $L$ is coercive. Let $2(\alpha_{1}+\alpha_{2})<\rho$
and the sequence $\\{z_{k}\\}_{k\in\mathbb{N}}$ be generated by Algorithm 3.1.
Then,

(i) the sequence $\\{z_{k}\\}_{k\in\mathbb{N}}$ is bounded and the set
$\omega(z_{k})$ is nonempty;

(ii) for any $\tilde{z}\in\omega(z_{k})$ and $z_{k_{j}}\rightarrow\tilde{z}$
as $j\rightarrow\infty$, the sequence
$H(z_{k_{j}},z_{k_{j}-1},z_{k_{j}-2})\rightarrow
H(\tilde{z},\tilde{z},\tilde{z})$ as $j\rightarrow\infty$ and $0\in\partial
L(\tilde{z})$.

###### Proof.

(i) By Lemma 3.1 and $H(z_{0},z_{-1},z_{-2})=L(z_{0})$, we have

$L(z_{k})\leq H(z_{k},z_{k-1},z_{k-2})\leq H(z_{0},z_{-1},z_{-2})=L(z_{0})$

for $k\geq 0$. It is clear that the whole sequence
$\\{z_{k}\\}_{k\in\mathbb{N}}$ is contained in the level set
$\\{z\in\mathbb{R}^{l}\times\mathbb{R}^{m}:\inf_{z\in\mathbb{R}^{l}\times\mathbb{R}^{m}}L(z)\leq
L(z)\leq L(z_{0})\\}$. It follows from the coercivity of $L$ and
$\inf_{z\in\mathbb{R}^{l}\times\mathbb{R}^{m}}L(z)>-\infty$ that the sequence
$\\{z_{k}\\}_{k\in\mathbb{N}}$ is bounded, and so the set $\omega(z_{k})$ is
nonempty.

(ii) For any $\tilde{z}\in\omega(z_{k})$ and $z_{k_{j}}\rightarrow\tilde{z}$
as $j\rightarrow\infty$, we show that $L(z_{k_{j}})\rightarrow L(\tilde{z})$
as $j\rightarrow\infty$. Let $\tilde{z}=(\tilde{x},\tilde{y})$. It following
from Lemma 3.1 (ii) that

$\lim_{k\rightarrow\infty}\triangle_{k}=\lim_{k\rightarrow\infty}\|z_{k}-z_{k-1}\|_{2}=0,$
$None$

and so
$\lim_{k\rightarrow\infty}\|x_{k}-x_{k-1}\|_{2}=\lim_{k\rightarrow\infty}\|y_{k}-y_{k-1}\|_{2}=0$.
From the definition of $x_{k}$, we have

$\displaystyle f(x_{k})+\langle
x_{k},\nabla_{x}Q(x_{k},y_{k-1})\rangle+D_{\phi}(x_{k},x_{k-1})+\alpha_{1,k-1}\langle
x_{k},x_{k-2}-x_{k-1}\rangle$ $\displaystyle+\alpha_{2,k-1}\langle
x_{k},x_{k-3}-x_{k-2}\rangle$ $\displaystyle\leq$ $\displaystyle
f(\tilde{x})+\langle\tilde{x},\nabla_{x}Q(\tilde{x},y_{k-1})\rangle+D_{\phi}(\tilde{x},x_{k-1})+\alpha_{1,k-1}\langle\tilde{x},x_{k-2}-x_{k-1}\rangle$
$\displaystyle+\alpha_{2,k-1}\langle\tilde{x},x_{k-3}-x_{k-2}\rangle,$

which implies that

$\displaystyle f(x_{k_{j}})+\langle
x_{k_{j}},\nabla_{x}Q(x_{k_{j}},y_{k_{j}-1})\rangle+D_{\phi}(x_{k_{j}},x_{k_{j}-1})+\alpha_{1,k_{j}-1}\langle
x_{k_{j}},x_{k_{j}-2}-x_{k_{j}-1}\rangle$
$\displaystyle+\alpha_{2,k_{j}-1}\langle
x_{k_{j}},x_{k_{j}-3}-x_{k_{j}-2}\rangle$ $\displaystyle\leq$ $\displaystyle
f(\tilde{x})+\langle\tilde{x},\nabla_{x}Q(\tilde{x},y_{k_{j}-1})\rangle+D_{\phi}(\tilde{x},x_{k_{j}-1})+\alpha_{1,k_{j}-1}\langle\tilde{x},x_{k_{j}-2}-x_{k_{j}-1}\rangle$
$\displaystyle+\alpha_{2,k_{j}-1}\langle\tilde{x},x_{k_{j}-3}-x_{k_{j}-2}\rangle.$

Letting $j\overset{}{\rightarrow}\infty$, we have

$\limsup_{j\rightarrow\infty}f(x_{k_{j}})\leq f(\tilde{x}).$

By the lower semicontinuity of $f$, it follows that

$f(\tilde{x})\leq\liminf_{j\rightarrow\infty}f(x_{k_{j}}),$

and so $\lim_{j\rightarrow\infty}f(x_{k_{j}})=f(\tilde{x})$. Similarly, we
have $\lim_{j\rightarrow\infty}g(y_{k_{j}})=g(\tilde{y})$, and hence
$\lim_{j\rightarrow\infty}L(z_{k_{j}})=L(\tilde{z})$. By (3.20) and

$H(z_{k_{j}},z_{k_{j}-1},z_{k_{j}-2})=L(z_{k_{j}})+\frac{\alpha_{1}+\alpha_{2}}{2}\|z_{k_{j}}-z_{k_{j}-1}\|_{2}^{2}+\frac{\alpha_{2}}{2}\|z_{k_{j}-1}-z_{k_{j}-2}\|_{2}^{2},$

we have
$\lim_{j\rightarrow\infty}H(z_{k_{j}},z_{k_{j}-1},z_{k_{j}-2})=L(\tilde{z})=H(\tilde{z},\tilde{z},\tilde{z})$.
It follows from (3.18) and (3.19) that

$s_{k}\in\partial L(x_{k},y_{k}),$

where

$\displaystyle s_{k}=$
$\displaystyle(\nabla_{x}Q(x_{k},y_{k})-\nabla_{x}Q(x_{k-1},y_{k-1})+\nabla\phi_{1}(x_{k-1})-\nabla\phi_{1}(x_{k}),\nabla_{y}Q(x_{k},y_{k})-\nabla_{y}Q(x_{k},y_{k-1})$
$\displaystyle+\nabla\phi_{2}(y_{k-1})-\nabla\phi_{2}(y_{k}))+(\alpha_{1,k-1}(x_{k-1}-x_{k-2})+\alpha_{2,k-1}(x_{k-2}-x_{k-3}),$
$\displaystyle\beta_{1,k-1}(y_{k-1}-y_{k-2})+\beta_{2,k-1}(y_{k-2}-y_{k-3})).$

Since $\nabla Q$ is Lipschitz continuous on bounded subsets and
$\nabla\phi_{i}$ is Lipschitz continuous ($i=1,2$), it follows that
$s_{k}\rightarrow 0$ as $k\rightarrow\infty$. By the closure property of the
subdifferential $\partial L$, we have $0\in\partial L(\tilde{z})$. ∎

Now, using Lemma 2.2, we can verify the convergence of the sequence
$\\{(x_{k},y_{k})\\}_{k\in\mathbb{N}}$ generated by Algorithm 3.1.

######  Theorem 3.1.

Let $2(\alpha_{1}+\alpha_{2})<\rho$ and $\\{(x_{k},y_{k})\\}_{k\in\mathbb{N}}$
be generated by Algorithm 3.1. Suppose that the function $L$ is coercive and
$z_{k}=(x_{k},y_{k})$ for all $k\in\mathbb{N}$. Then, the sequence
$\\{(z_{k+1},z_{k},z_{k-1})\\}_{k\in\mathbb{N}}$ satisfies (H1), (H2), and
(H3) for the function

$H(u,v,w):\mathbb{R}^{3(l+m)}\rightarrow{(-\infty,+\infty]},$

which is defined by (3.3). Moreover, if $H(u,v,w)$ has the KL property at a
cluster point $(z^{\ast},z^{\ast},z^{\ast})$ of the sequence
$\\{(z_{k+1},z_{k},z_{k-1})\\}_{k\in\mathbb{N}}$. Then the sequence
$\\{z_{k}\\}_{k\in\mathbb{N}}$ has finite length, $z_{k}\rightarrow z^{\ast}$
as $k\rightarrow\infty$, and $(z^{\ast},z^{\ast},z^{\ast})$ is a critical
point of $H$. Hence, $z^{\ast}$ is a critical point of $L$.

###### Proof.

First, we verify that assumptions (H1), (H2) and (H3) are satisfied. We
consider the sequence $u_{k}=(z_{k},z_{k-1},z_{k-2})$ for all $k\in\mathbb{N}$
and the proper lower semicontinuous function $F=H$. Condition (H1) is proved
in Lemma 3.1 (i). To prove condition (H2), from Lemma 3.2 we have

$w_{k+1}\in\partial H(z_{k+1},z_{k},z_{k-1})=\partial
H(x_{k+1},y_{k+1},x_{k},y_{k},x_{k-1},y_{k-1}),$

where $w_{k+1}$ is defined by (3.13) for $k\in\mathbb{N}$. Lemma 3.3 (i)
implies that the sequence $\\{z_{k}\\}_{k\in\mathbb{N}}$ is bounded. Let
$c=\max\\{L_{\nabla\phi_{1}},L_{\nabla\phi_{2}}\\},d=M+L_{2}(x_{k+1})$. It
follows from

$(\|x_{k+1}-x_{k}\|_{2}+\|y_{k+1}-y_{k}\|_{2})^{2}\leq
2(\|x_{k+1}-x_{k}\|_{2}^{2}+\|y_{k+1}-y_{k}\|_{2}^{2})=2\|z_{k+1}-z_{k}\|_{2}^{2}$

and (3.15) that

$\displaystyle\|w_{k+1}\|_{2}$ $\displaystyle\leq$
$\displaystyle({c+d+2\alpha_{1}+2\alpha_{2}})\sqrt{2}\|z_{k+1}-z_{k}\|_{2}+(\alpha_{1}+2\alpha_{2})\sqrt{2}\|z_{k}-z_{k-1}\|_{2}+\alpha_{2}\sqrt{2}\|z_{k-1}-z_{k-2}\|_{2},$

which implies that condition (H2) holds for
$b=\sqrt{2}(c+d+\alpha_{1}+2\alpha_{2})$. From Lemma 3.3 (i), we have
$\omega(z_{k})\neq\emptyset$, and so there exists a subsequence
$\\{z_{k_{j}}\\}_{j\in\mathbb{N}}\subseteq\\{z_{k}\\}_{k\in\mathbb{N}}$ such
that $z_{k_{j}}\rightarrow\tilde{z}$ as $j\rightarrow\infty$. By Lemma 3.1
(ii), we know $\triangle_{k_{j}}=\|z_{k_{j}}-z_{k_{j}-1}\|_{2}\rightarrow 0$
as $j\rightarrow\infty$. Taking
$u_{k_{j}}=(z_{k_{j}},z_{k_{j}-1},z_{k_{j}-2})$ and
$\tilde{u}=(\tilde{z},\tilde{z},\tilde{z})$, we have
$u_{k_{j}}\rightarrow\tilde{u}$ as $j\rightarrow\infty$. Lemma 3.3 (ii)
implies that

$\lim_{j\rightarrow\infty}F(u_{k_{j}})=\lim_{j\rightarrow\infty}H(z_{k_{j}},z_{k_{j}-1},z_{k_{j}-2})=H(\tilde{z},\tilde{z},\tilde{z})=F(\tilde{u})$

and $0\in\partial L(\tilde{z})$. So we get that (H3) holds. Now, Lemma 2.2
concludes the proof. ∎

##  4 Numerical experiments

In all of the following experiments, the codes are written in MATLAB and are
performed on a personal Dell computer with Pentium(R) Dual-Core CPU @ 2.4GHz
and RAM 2.00GB.

###  4.1 Sparse nonnegative matrix factorization

Given a matrix $A$, the target of nonnegative matrix factorization (NMF) [33,
34, 35] is to decompose database $A$ into a number of sparse basis, such that
each database can be approximated by using a small number of those parts,
i.e., seeking a factorization $A\approx XY$, where $X\in\mathbb{R}^{n\times
r}$, $Y\in\mathbb{R}^{r\times d}$ are nonnegative with $r\leq d$ and $X$
sparse. In order to enforce sparsity on the basis, we additionally consider
$l_{0}$ sparsity constraint [36] in this problem. Thus sparse-NMF can be
formulated as the following model

$\displaystyle\underset{X,Y}{\min}\left\\{\left\|A-XY\right\|_{F}^{2}:\
X,Y\geq 0,\ \left\|X_{i}\right\|_{0}\leq s,\ i=1,2,\dots,r\right\\}.$ (4.1)

Here, $X_{i}$ denotes the ith column of $X$. In dictionary learning and sparse
coding, $X$ is called the learned dictionary with coefficients $Y$. In this
formulation, the sparsity on $X$ is strictly enforced using the nonconvex
$l_{0}$ constraint, restricting $75\%$ of the entries to be 0. This model
(4.1) can be converted to (1.1), where
$Q(X,Y)=\frac{\lambda}{2}\left\|A-XY\right\|_{F}^{2}$,
$f(X)=\sum_{i=1}^{r}\left\|X_{i}\right\|_{0}+\iota_{X\geq 0}(X)$,
$g(Y)=\iota_{Y\geq 0}(Y)$, $\iota_{C}$ is the indicator function on $C$.

We use the extended Yale-B dataset and the ORL dataset, which are standard
facial recognition benchmarks consisting of human face images111
http://www.cad.zju.edu.cn/home/dengcai/Data/FaceData.html.. The extended
Yale-B dataset contains 2414 cropped images of size $32\times 32$, while the
ORL dataset contains 400 images sized $64\times 64$, see Figure 1. For solving
this sparse NMF problem (4.1), [22, 23] gave the details on how to solve the
X-subproblems and Y-subproblems. In the experiment for the Yale dataset, we
extract 49 sparse basis images for the dataset. For the ORL dataset we extract
25 sparse basis images. We take parameter $\lambda=0.5$. In Algorithm 3.1, let
$\phi_{1}(X)=\frac{\mu_{1}}{2}\left\|X\right\|_{2}^{2}$,
$\phi_{2}(Y)=\frac{\mu_{2}}{2}\left\|Y\right\|_{2}^{2}$. In numerical
experiment, we computed $\mu_{1}$ and $\mu_{2}$ by computing the largest
eigenvalues of $\lambda YY^{T}$ and $\lambda X^{T}X$ at $k$-th iteration,
respectively. We use the algorithm PALM, iPALM, GiPALM and Algorithm 3.1 for a
specific sparsity setting $s=0.25$. We choose $\alpha_{1k}=\beta_{1k}=0.2$,
$\alpha_{2k}=\beta_{2k}=0.3$ in Algorithm 3.1 and $\alpha_{1k}=\beta_{1k}=0.5$
in iPALM and GiPALM. The related numerical results of objective values are
reported in Figure 2 and Figure 4 under Yale-B dataset and ORL dataset,
respectively. And we also take extrapolation parameter
$\alpha_{1k}=\alpha_{2k}=\beta_{1k}=\beta_{2k}=\frac{k-1}{k+2}$ in Algorithm
3.1 and $\alpha_{1k}=\beta_{1k}=\frac{k-1}{k+2}$ in iPALM and GiPALM. In
Figure 3 and Figure 5, we report the numerical results under Yale-B dataset
and ORL dataset, respectively. One can observe from these four figures that
the proposed Algorithm 3.1 can get slightly lower values than the other three
algorithms within almost the same computation time.

We also compare Algorithm 3.1 with PALM, iPALM and GiPALM for different
sparsity settings (the value of $s$). We take extrapolation parameter
$\alpha_{1k}=\alpha_{2k}=\beta_{1k}=\beta_{2k}=\frac{k-1}{k+2}$ in Algorithm
3.1 and $\alpha_{1k}=\beta_{1k}=\frac{k-1}{k+2}$ in iPALM and GiPALM. The
results of the basis images are shown in Figure 6. One can observe from Figure
6 that for smaller values of $s$, the four algorithms lead to more compact
representations. This might improve the generalization capabilities of the
representation.

(a)  Figure 1: ORL face database which includes 400 normalized cropped frontal
faces which we used in our S-NMF example.

(a) Epoch counts comparison

(b) Wall-clock time comparison

Figure 2: Objective decrease comparison of sparse-NMF with $s=25\%$ on Yale
dataset.

(a) Epoch counts comparison

(b) Wall-clock time comparison

Figure 3: Objective decrease comparison of sparse-NMF with $s=25\%$ on Yale
dataset.

(a) Epoch counts comparison

(b) Wall-clock time comparison

Figure 4: Objective decrease comparison of sparse-NMF with $s=25\%$ on ORL
dataset.

(a) Epoch counts comparison

(b) Wall-clock time comparison

Figure 5: Objective decrease comparison of sparse-NMF with $s=25\%$ on ORL
dataset.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 6: The results for 25 basis faces using different sparsity settings.
From left column to right column are the results of PALM, iPALM, GiPALM and
Algorithm 3.1, respectively. From top row to bottom row are the result of
$s=25\%$ (sparsity of $s=25\%$ means that each basis face contains $75\%$
nonzero pixels), $s=33\%$ and $s=50\%$, respectively. Clearly, stronger
sparsity leads to a more compact representation.

###  4.2 Signal recovery

The sparse signal recovery problem has been studied extensively. Supposed $x$
is an unknown vector in $\mathbb{R}^{m}$ (a signal), given a matrix
$A\in\mathbb{R}^{n\times m}$, and an observation $b\in\mathbb{R}^{n}$, we plan
to recover $x$ from observation $b$ such that $x$ is of the sparsest structure
(that is, $x$ has the fewest nonzero components). In this section, the sparse
signal recovery problem can be modulated by the following $L_{0}$-problem

$\min_{x}\\\ {\rm}\ \ \left\|x\right\|_{0}$ $s.t.\\\ {\rm}\ \ Ax=b.$

The $L_{0}$-problem in the form of $L_{1/2}$ regularization can be described
as

$\min_{x}\\\ {\rm}\ \
\frac{1}{2}\left\|Ax-b\right\|_{2}^{2}+\eta\left\|x\right\|_{1/2}^{1/2},$
$None$

where $\eta>0$ is used to balance regularization and data fitting.
$\left\|x\right\|_{1/2}$ is the $L_{1/2}$ quasi-norm of $\mathbb{R}^{m}$,
defined by
$\left\|x\right\|_{1/2}=({\textstyle\sum_{i=1}^{m}}\left|x_{i}\right|^{1/2})^{2}$.
Now we illustrate how to implement the Algorithm 3.1 for solving the above
model.
By introducing an auxiliary variable $y\in\mathbb{R}^{m}$. The model (4.2) can
be reformulated as

$\min_{x,y}\\\ {\rm}\ \
\frac{1}{2}\left\|Ax-b\right\|_{2}^{2}+\eta\left\|y\right\|_{1/2}^{1/2},$
$s.t.\\\ {\rm}\ \ x=y,$

which indeed can be described as

$\min_{x,y}\\\ {\rm}\ \
\frac{1}{2}\left\|Ax-b\right\|_{2}^{2}+\eta\left\|y\right\|_{1/2}^{1/2}+\frac{\gamma}{2}\left\|x-y\right\|_{2}^{2}.$
$None$

When we set
$f(x)=\frac{1}{2}\left\|Ax-b\right\|_{2}^{2},Q(x,y)=\frac{\gamma}{2}\left\|x-y\right\|_{2}^{2},g(y)=\eta\left\|y\right\|_{1/2}^{1/2}$
, then (4.3) becomes nonconvex minimization problem (1.1). Let
$\phi_{1}(x)=\langle
x,Mx\rangle,\phi_{2}(y)=\frac{\lambda}{2}\left\|y\right\|_{2}^{2}$, then we
can obtain colsed form of solution for solving the $x$-subproblem and
$y$-subproblem of Algorithm 3.1 respectively as follows.
Setting $M=\mu I-A^{T}A$, the $x$-subproblem corresponds to the following
optimization problem

$\displaystyle x_{k+1}$ $\displaystyle\in\arg\min_{x\in\mathbb{R}^{l}}$
$\displaystyle\left\\{\frac{1}{2}\left\|Ax-b\right\|_{2}^{2}+\langle
x,\gamma(x_{k}-y_{k})\rangle+\frac{1}{2}\left\|x-x_{k}\right\|_{M}^{2}+\alpha_{1k}\langle
x,x_{k-1}-x_{k}\rangle\right.$ $\displaystyle\left.+\alpha_{2k}\langle
x,x_{k-2}-x_{k-1}\rangle\right\\}$
$\displaystyle=\arg\min_{x\in\mathbb{R}^{l}}$
$\displaystyle\left\\{\frac{1}{2}\left\|Ax\right\|_{2}^{2}-\left\langle
Ax,b\right\rangle+\langle x,\gamma(x_{k}-y_{k})\rangle+\frac{1}{2}\left\langle
x-x_{k},(\mu I-A^{T}A)(x-x_{k})\right\rangle\right.$
$\displaystyle\left.+\alpha_{1k}\langle
x,x_{k-1}-x_{k}\rangle+\alpha_{2k}\langle x,x_{k-2}-x_{k-1}\rangle\right\\}$
$\displaystyle=\arg\min_{x\in\mathbb{R}^{l}}$
$\displaystyle\left\\{\frac{\mu}{2}\left\|x\right\|_{2}^{2}-\langle x,\mu
x_{k}-A^{T}Ax_{k}\rangle-\left\langle x,A^{T}b\right\rangle+\langle
x,\gamma(x_{k}-y_{k})\rangle\right.$ $\displaystyle\left.+\alpha_{1k}\langle
x,x_{k-1}-x_{k}\rangle+\alpha_{2k}\langle x,x_{k-2}-x_{k-1}\rangle\right\\},$

which has an explicit expression

$x_{k+1}=\frac{1}{\mu}\left[\mu
x_{k}-A^{T}Ax_{k}+A^{T}b-\gamma(x_{k}-y_{k})+\alpha_{1k}(x_{k}-x_{k-1})+\alpha_{2k}(x_{k-1}-x_{k-2})\right].$

The $y$-subproblem corresponds to the following optimization problem

$\displaystyle y_{k+1}$ $\displaystyle\in\arg\min_{y\in\mathbb{R}^{m}}$
$\displaystyle\left\\{\eta\left\|y\right\|_{1/2}^{1/2}+\langle
y,\gamma(y_{k}-x_{k+1})\rangle+\frac{\lambda}{2}\left\|y-y_{k}\right\|_{2}^{2}+\beta_{1k}\langle
y,y_{k-1}-y_{k}\rangle\right.$ $\displaystyle\left.+\beta_{2k}\langle
y,y_{k-2}-y_{k-1}\rangle\right\\}$
$\displaystyle=\arg\min_{y\in\mathbb{R}^{m}}$
$\displaystyle\left\\{\eta\left\|y\right\|_{1/2}^{1/2}+\frac{\lambda}{2}\left\|y-\frac{1}{\lambda}\left[\lambda
y_{k}+\gamma(x_{k+1}-y_{k})+\beta_{1k}(y_{k}-y_{k-1})\right.\right.\right.$
$\displaystyle\left.\left.\left.+\beta_{2k}(y_{k-1}-y_{k-2})\right]{}\right\|^{2}\right\\}$
$\displaystyle=\mathcal{H}$
$\displaystyle\left(y_{k}+\frac{1}{\lambda}\left[\gamma(x_{k+1}-y_{k})+\beta_{1k}(y_{k}-y_{k-1})+\beta_{2k}(y_{k-1}-y_{k-2})\right],\frac{\eta}{\lambda}\right),$

where, for any $\kappa>0$, $\mathcal{H}(\cdot,\kappa)$ is called the half
shrinkage operator [37] defined as

${\mathcal{H}\rm}(a;\kappa)=\left\\{h_{\kappa}(a_{1}),h_{\kappa}(a_{2}),\dots,h_{\kappa}(a_{n})\right\\}^{T}$

with

$h_{\kappa}(a_{i{\tiny}})=\left\\{\begin{matrix}\frac{2a_{i}}{3}(1+cos(\frac{2\pi}{3}-\frac{2}{3}\varphi(a_{i}))),&\left|a_{i}\right|>\frac{3}{4}\kappa^{2/3},\\\
0,&otherwise\end{matrix}\right.$

and
$\varphi(a_{i})=$arccos$\left(\frac{\kappa}{8}\left(\frac{\left|a_{i}\right|}{3}\right)^{-3/2}\right).$

In this experiment, each entry of $A$ is drawn from the standard normal
distribution, and then all columns of $A$ are normalized according to
$L_{1/2}$ quasi-norm, noted that $\left\|A\right\|\leq 1$; we genergte a
random sparse vector as $x$; the noise vector $\omega\sim
N\left(0,10^{-3}I\right)$ and the vector $b=Ax+\omega$; the regularization
parameter $\eta=0.001\left\|A^{T}b\right\|_{\infty}$ and the penalization
parameter is set as $\gamma=0.2$. By Lemma 3.1, we take
$\alpha_{1}>0,\alpha_{2}>0$ such that $2(\alpha_{1}+\alpha_{2})<\rho$, where
$\rho=\min\\{\mu-\left\|A\right\|^{2}-\gamma,\lambda-\gamma\\}$,
$\mu=2,\lambda=1.5$. In Algorithm 3.1 and TiBAM with two-step inertial
extrapolation, we select inertial parameters
$\alpha_{1k}=\alpha_{2k}=\beta_{1k}=\beta_{2k}=0.99\rho/4$; in iBPALM with
one-step inertial extrapolation, we select inertial parameters
$\alpha_{1k}=\beta_{1k}=0.99\rho/2$; BPALM is the special case of Algorithm
3.1 without inertial extrapolation, that is,
$\alpha_{1k}=\alpha_{2k}=\beta_{1k}=\beta_{2k}=0$; We take the origin point as
the initial point for all algorithms and use

$E_{k}=\|x_{k+1}-x_{k}\|+\|y_{k+1}-y_{k}\|<10^{-4}$

as the stopping criterion. In the following, we compare Algorithm 3.1 with
TiBAM, iBPALM and BPALM for matrix $A$ of different dimensions with
$(n,m)=(40,200)$ and $(n,m)=(100,500)$.

Table 1: Compare Algorithm 3.1 with TiBAM, iBPALM and BPALM for $n=40,m=200$.  Algorithm |  | $b=Ax$ |  |  |  | $b=Ax+\omega$ |
---|---|---|---|---|---|---|---
n=40,m=200 | Iter. | Time(s) | $\left\|x_{k}-y_{k}\right\|$ |  | Iter. | Time(s) | $\left\|x_{k}-y_{k}\right\|$
Alg. 3.1 | 713 | 4.1934 | $5.0602\times 10^{-5}$ |  | 810 | 6.5243 | $6.5298\times 10^{-5}$
TiBAM | 1110 | 7.6269 | $8.9963\times 10^{-5}$ |  | 1230 | 9.9036 | $9.1096\times 10^{-5}$
iBPALM | 1378 | 9.8594 | $9.3004\times 10^{-5}$ |  | 1577 | 12.0625 | $1.0239\times 10^{-4}$
BPALM | 2033 | 12.8073 | $1.7996\times 10^{-4}$ |  | 2276 | 14.7630 | $1.9939\times 10^{-4}$

(a) $b=Ax$

(b) $b=Ax+\omega$

Figure 7: The value of $\|x_{k+1}-x_{k}\|+\|y_{k+1}-y_{k}\|$ versus the
iteration numbers for 4.2 with $n=40,m=200$.

Table 1 and Table 2 report the number of iteration, CPU time and the 2-norm of
$x_{k}-y_{k}$ for $(n,m)=(40,200)$ and $(n,m)=(100,500)$, respectively. In the
numerical results, “Iter.” denotes the number of iterations. “Time” denotes
the CPU time. It is obvious that Algorithm 3.1 has an advantage over TiBAM,
iBPALM and BPALM in terms of iteration number and time for solving the above
problem. In Figure 7 and Figure 8, the left picture records the downward trend
of the error function without noise $(b=Ax)$ and the right picture shows the
trend of the error function with noise $(b=Ax+\omega)$ for $(n,m)=(40,200)$
and $(n,m)=(100,500)$, respectively. They show the efficiency and advantage of
Algorithm 3.1.

Table 2: Compare Algorithm 3.1 with TiBAM, iBPALM and BPALM for $n=100,m=500$. Algorithm |  | $b=Ax$ |  |  |  | $b=Ax+\omega$ |
---|---|---|---|---|---|---|---
n=100,m=500 | Iter. | Time(s) | $\left\|x_{k}-y_{k}\right\|$ |  | Iter. | Time(s) | $\left\|x_{k}-y_{k}\right\|$
Alg. 3.1 | 1610 | 40.3756 | $9.9978\times 10^{-5}$ |  | 1920 | 46.5243 | $1.5527\times 10^{-4}$
TiBAM | 2108 | 48.5674 | $1.3827\times 10^{-4}$ |  | 2467 | 52.6250 | $1.8418\times 10^{-4}$
iBPALM | 2732 | 55.7853 | $2.0018\times 10^{-4}$ |  | 3196 | 59.1563 | $2.1329\times 10^{-4}$
BPALM | 3731 | 65.3309 | $2.9998\times 10^{-4}$ |  | 4023 | 69.7543 | $3.0698\times 10^{-4}$

(a) $b=Ax$

(b) $b=Ax+\omega$

Figure 8: The value of $\|x_{k+1}-x_{k}\|+\|y_{k+1}-y_{k}\|$ versus the
iteration numbers for 4.2 with $n=100,m=500$.

###  4.3 Nonconvex quadratic fractional programming

In this subsection, we provide some numerical experiments which we carried out
in order to illustrate the numerical convergence of Algorithm 3.1 with
different Bregman distances. The following list are various functions with its
Bregman distances:

(i) Define the function $\varphi_{1}(x)=\mu\sum_{i=1}^{m}x_{i}\ln x_{i}$ with
domain

$\text{dom}\varphi_{1}=\\{x=(x_{1},x_{2},\cdots,x_{m})^{T}\in\mathbb{R}^{m}:x_{i}>0,i=1,2,\cdots,m\\}$

and range ran$\varphi_{1}=(-\infty,+\infty)$. Then

$\nabla\varphi_{1}(x)=\mu(1+\ln(x_{1}),1+\ln(x_{2}),\cdots,1+\ln(x_{m}))^{T}$

and the Bregman distance (the Kullback-Leibler distance) with respect to
$\varphi_{1}$ is

$D_{\varphi_{1}}(x,y)=\mu\sum_{i=1}^{m}\big{(}x_{i}\ln\big{(}\frac{x_{i}}{y_{i}}\big{)}+y_{i}-x_{i}\big{)},\
\ \forall x,y\in\mathbb{R}_{++}^{m}.$

(ii) Define the function $\varphi_{2}(x)=-\mu\sum_{i=1}^{m}\ln x_{i}$ with
domain

$\text{dom}\varphi_{2}=\\{x=(x_{1},x_{2},\cdots,x_{m})^{T}\in\mathbb{R}^{m}:x_{i}>0,i=1,2,\cdots,m\\}$

and range ran$\varphi_{2}=(-\infty,+\infty)$. Then

$\nabla\varphi_{2}(x)=\mu(-\frac{1}{x_{1}},-\frac{1}{x_{2}},\cdots,-\frac{1}{x_{m}})^{T}$

and the Bregman distance (the Itakura-Saito distance) with respect to
$\varphi_{2}$ is

$D_{\varphi_{2}}(x,y)=\mu\sum_{i=1}^{m}\big{(}\frac{x_{i}}{y_{i}}-\ln\big{(}\frac{x_{i}}{y_{i}}\big{)}-1\big{)},\
\ \forall x,y\in\mathbb{R}_{++}^{m}.$

(iii) Define the function $\varphi_{3}(x)=\frac{\mu}{2}\|x\|^{2}$ with domain
$\text{dom}\varphi_{3}=\mathbb{R}^{m}$ and range ran$\varphi_{3}=[0,+\infty)$.
Then $\nabla\varphi_{3}(x)=x$ and the Bregman distance (the squared Euclidean
distance) with respect to $\varphi_{3}$ is

$D_{\varphi_{3}}(x,y)=\frac{\mu}{2}\|x-y\|^{2},\ \ \forall
x,y\in\mathbb{R}^{m}.$

It is clear that $\varphi_{i}$ is $1$-strongly convex ($i=1,2,3$).

We consider the quadratic fractional programming problem

$\min_{x\in C}\\\ f(x):=\frac{x^{T}Mx+a^{T}x+c}{b^{T}x+d}$

with

$C=\left\\{x\in\mathbb{R}^{m}:1\leq x_{i}\leq 3,i=1,2,\cdots,m\right\\},$

where $M:\mathbb{R}^{m}\to\mathbb{R}^{m}$ is a bounded linear operator,
$a\in\mathbb{R}^{m}$, $b\in\mathbb{R}^{m}$, $c=-2$ and $d=20$. By [5], we know
$f$ is pseudo-convex on the open set
$E=\left\\{x\in\mathbb{R}^{m}:b^{T}x+d>0\right\\}$, if $C\subseteq E$, then
$f$ is nonconvex.
The quadratic fractional programming problem can be rewritten as (1.1):

$\min\\{f(x)+\frac{\gamma}{2}\|x-y\|^{2}_{2}+\iota_{C}(y):x,\
y\in\mathbb{R}^{m}\\},$

where $\iota_{C}$ is the indicator function on $C$.
We now elaborate the $x$-subproblem and $y$-subproblem of Algorithm 3.1
respectively as follows.

$\begin{cases}\begin{aligned}
x_{k+1}\in\arg\min_{x\in\mathbb{R}^{m}}\\{&f(x)+\langle
x,\gamma(x_{k}-y_{k})\rangle+D_{\phi_{1}}(x,x_{k})+\alpha_{1k}\langle
x,x_{k-1}-x_{k}\rangle\\\ &+\alpha_{2k}\langle x,x_{k-2}-x_{k-1}\rangle\\},\\\
y_{k+1}\in\arg\min_{y\in\mathbb{R}^{m}}\\{&\iota_{C}(y)+\langle
y,\gamma(y_{k}-x_{k+1})\rangle+D_{\phi_{2}}(y,y_{k})+\beta_{1k}\langle
y,y_{k-1}-y_{k}\rangle\\\ &+\beta_{2k}\langle
y,y_{k-2}-y_{k-1}\rangle\\}.\end{aligned}\end{cases}$

In this experiment, we take $\alpha_{1k}=\beta_{1k}=0.2$,
$\alpha_{2k}=\beta_{2k}=0.3,\gamma=10,\mu=36$. We perform the numerical tests
of Algorithms 3.1 with different Bregman distance and use

$E_{k}=\|x_{k+1}-x_{k}\|+\|y_{k+1}-y_{k}\|<10^{-4}$

as the stopping criterion. We use “Alg(ij)” to denote Algorithm 3.1 with
$\phi_{1}(x)=\varphi_{i}(x)$ and $\phi_{2}(x)=\varphi_{j}(x)$ $(1\leq i\leq
3,1\leq j\leq 3)$. For the above quadratic fractional programming problem, we
will give the numerical results of Algorithm 3.1 for different matrix $M$ and
dimensions. We randomly selected the starting point, carried out 30 times
randomly, and obtained the averaged iteration number and the averaged CPU
time. The numerical results for the performance of Algorithm 3.1 with
different Bregman distance are shown in Table 3 and Table 4. In the numerical
results, “Iter.” denotes the number of outer iterations. “InIterx.” and
“InItery.” denote, respectively, the number of inner iterations for solving
$x-$subproblem and $y-$subproblem. “Time” denotes the CPU time.

Problem 1. We consider a fixed matrix $M=\begin{bmatrix}5&-1&2&0&2\\\
-1&6&-1&3&0\\\ 2&-1&3&0&1\\\ 0&3&0&5&0\\\ 2&0&1&0&4\end{bmatrix}$, and give
the vectors $a=(1,2,-1,-2,1)^{T},b=(1,0,-1,0,1)^{T}$, in this case,
$C\subseteq E$. For different bregman distances, the numerical results are
given for one-step inertial and two-step inertial extrapolation in Table 3
below, where the one-step inertial extrapolation $\alpha_{1k}=\beta_{1k}=0.5$.

Table 3: Numerical results of Algorithm 3.1 for fixed matrix $M$ | Alg.(11) | Alg.(12) | Alg.(13) | Alg.(21) | Alg.(22) | Alg.(23) | Alg.(31) | Alg.(32) | Alg.(33)
---|---|---|---|---|---|---|---|---|---
$one-step$ |  |  |  |  |  |  |  |  |
Iter. | 662 | 675 | 657 | 307 | 308 | 305 | 3649 | 3497 | 3918
InIterx. | 992 | 1003 | 988 | 493 | 497 | 490 | 3812 | 3724 | 4090
InItery. | 704 | 719 | 699 | 307 | 310 | 305 | 3907 | 3731 | 4237
Time(s) | 0.4063 | 0.3750 | 0.2031 | 0.1205 | 0.2188 | 0.0718 | 0.9188 | 0.8563 | 1.1938
$two-step$ |  |  |  |  |  |  |  |  |
Iter. | 529 | 534 | 523 | 294 | 295 | 293 | 2984 | 2828 | 3513
InIterx. | 830 | 834 | 825 | 493 | 495 | 491 | 3237 | 3114 | 3698
InItery. | 559 | 565 | 552 | 294 | 296 | 292 | 3193 | 3019 | 3765
Time(s) | 0.2813 | 0.3281 | 0.1563 | 0.0938 | 0.1938 | 0.0681 | 0.8063 | 0.7501 | 1.0531

It can be seen that the Kullback-Leibler distance and the Itakura-Saito
distance have computational advantage than the squared Euclidean distance for
our proposed Algorithm 3.1 in terms of number of iteration and CPU time. For
one-step inertial, the computation result shows that Alg(23)
($\phi_{1}(x)=\varphi_{2}(x)$, $\phi_{2}(x)=\varphi_{3}(x)$) has the best
performance while Alg(33) ($\phi_{1}(x)=\varphi_{3}(x)$,
$\phi_{2}(x)=\varphi_{3}(x)$) has the poorest performance for solving the
quadratic fractional programming problem. For two-step inertial, the
computation result shows that Alg(23) ($\phi_{1}(x)=\varphi_{2}(x)$,
$\phi_{2}(x)=\varphi_{3}(x)$) has the best performance while Alg(33)
($\phi_{1}(x)=\varphi_{3}(x)$, $\phi_{2}(x)=\varphi_{3}(x)$) has the poorest
performance for solving the quadratic fractional programming problem.

Problem 2. We randomly select matrix $M\in\mathbb{R}^{m\times m}$, and vectors
$a\in\mathbb{R}^{m},b\in\mathbb{R}^{m}$. For $m=5,20,50,100$, we give their
numerical results for Algorithm 3.1 with different Bregman distances in Table
4.

Table 4: Numerical results of Algorithm 3.1 for random matrix $M$ |  | Alg.(11) |  | Alg.(12) |  | Alg.(13) |  | Alg.(21) |  | Alg.(22) |  | Alg.(23) |  | Alg.(31) |  | Alg.(32) |  | Alg.(33)
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$m=5$ |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
Iter. |  | 431 |  | 432 |  | 430 |  | 106 |  | 111 |  | 105 |  | 2069 |  | 1760 |  | 2625
InIterx. |  | 812 |  | 812 |  | 811 |  | 406 |  | 410 |  | 405 |  | 2348 |  | 2078 |  | 2851
InItery. |  | 469 |  | 471 |  | 467 |  | 298 |  | 301 |  | 297 |  | 2286 |  | 1927 |  | 2935
Time(s) |  | 0.3125 |  | 0.2438 |  | 0.1813 |  | 0.1250 |  | 0.1156 |  | 0.0575 |  | 0.9844 |  | 0.7969 |  | 1.0156
$m=20$ |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
Iter. |  | 445 |  | 445 |  | 445 |  | 158 |  | 156 |  | 155 |  | 2921 |  | 2437 |  | 3654
InIterx. |  | 872 |  | 872 |  | 872 |  | 427 |  | 429 |  | 410 |  | 3223 |  | 2770 |  | 3916
InItery. |  | 471 |  | 470 |  | 470 |  | 309 |  | 308 |  | 307 |  | 3115 |  | 2588 |  | 3924
Time(s) |  | 0.9375 |  | 0.2813 |  | 0.1938 |  | 0.1388 |  | 0.1406 |  | 0.0781 |  | 1.1875 |  | 0.9531 |  | 1.2188
$m=50$ |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
Iter. |  | 468 |  | 468 |  | 468 |  | 315 |  | 312 |  | 310 |  | 3488 |  | 2812 |  | 4400
InIterx. |  | 907 |  | 907 |  | 907 |  | 483 |  | 478 |  | 475 |  | 3788 |  | 3149 |  | 4656
InItery. |  | 485 |  | 485 |  | 485 |  | 313 |  | 311 |  | 310 |  | 3704 |  | 2978 |  | 4695
Time(s) |  | 1.1250 |  | 0.3750 |  | 0.6250 |  | 0.1563 |  | 0.1750 |  | 0.1094 |  | 1.9063 |  | 1.8906 |  | 2.9219
$m=100$ |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
Iter. |  | 651 |  | 638 |  | 657 |  | 382 |  | 379 |  | 377 |  | 4692 |  | 3719 |  | 5529
InIterx. |  | 960 |  | 947 |  | 966 |  | 649 |  | 643 |  | 641 |  | 4903 |  | 3976 |  | 5712
InItery. |  | 694 |  | 679 |  | 700 |  | 385 |  | 380 |  | 377 |  | 4944 |  | 3911 |  | 5842
Time(s) |  | 2.8750 |  | 1.3906 |  | 2.0313 |  | 1.4375 |  | 0.7656 |  | 0.5938 |  | 3.6094 |  | 2.2500 |  | 4.6875

From Table 4, it can be seen from that the numerical results show that Alg(23)
($\phi_{1}(x)=\varphi_{2}(x)$, $\phi_{2}(x)=\varphi_{3}(x)$) performs best
when the dimension of the matrix $M$ is low, and Alg(12)
($\phi_{1}(x)=\varphi_{1}(x)$, $\phi_{2}(x)=\varphi_{2}(x)$) performs a little
better when the dimension of the matrix $M$ is slightly higher.

##  5 Conclusion

Based on the proximal alternating linearized minimization algorithm, two-step
inertial Bregman proximal alternating linearized minimization algorithm is
proposed to solve nonconvex and nonsmooth nonseparable optimization problems.
We construct appropriate benefit function and obtain that the generated
sequence is globally convergent to a critical point, under the assumptions
that the objective function satisfies the Kurdyka–Łojasiewicz inequality and
the parameters satisfy certain conditions. In numerical experiments, we choose
appropriate Bregman distance such that the solution of subproblems have closed
form for solving the sparse signal recovery problem. We also apply different
Bregman distance to solve quadratic fractional programming problem. Numerical
results are reported to show the effectiveness of the proposed algorithm.

## References

  * [1] Nikolova M., Ng M.K., Zhang S.Q., Ching W.K., Eﬃcient reconstruction of piecewise constant images using nonsmooth nonconvex minimization, SIAM J. Imaging Sci., 2008, 1(1), 2-25.
  * [2] Gu S.H., Zhang L., Zuo W.M., Feng X.C., Weighted nuclear norm minimization with application to image denoising, In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2014, 2862-2869.
  * [3] Bian W., Chen X.J., Linearly constrained non-Lipschitz optimization for image restoration, SIAM J. Imaging Sci., 2015, 8(4), 2294-2322.
  * [4] Bolte J., Sabach S., Teboulle M., Proximal alternating linearized minimization for nonconvex and nonsmooth problems, Math. Program., 2014, 146, 459-494.
  * [5] Bot R.I., Csetnek E.R., Vuong P.T., The forward-backward-forward method from continuous and discrete perspective for pseudo-monotone variational inequalities in Hilbert spaces, Eur. J. Oper. Res., 2020, 49-60.
  * [6] Attouch H., Bolte J., Svaiter B.F., Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized Guass-Seidel methods, Math. Program., 2013, 137, 91-129.
  * [7] Donoho D.L., Compressed sensing, IEEE Trans. Inform. Theory, 2006, 4, 1289-1306.
  * [8] Boyd S., Parikh N., Chu E., Peleato B., Eckstein J., Distributed optimization and statistical learning via the alternating direction method of multipliers, Foundations and Trends in Machine Learning, 2011, 3(1), 1-122.
  * [9] Bertsekas D.P., Tsitsiklis J.N., Parallel and Distributed Computation: Numerical Methods, Prentice hall, Englewood Cliffs, NJ, 1989.
  * [10] Bertsekas D.P., Nonlinear programming, J. Oper. Res. Soc., 1977, 48, 334-334.
  * [11] Beck A., Tetruashvili L., On the convergence of block coordinate descent type methods, SIAM J. Optim., 2013, 23, 2037-2060.
  * [12] Auslender A., Asymptotic properties of the Fenchel dual functional and applications to decomposition problems, J. Optim. Theory Appl., 1992, 73(3), 427-449.
  * [13] Attouch H., Bolte J., Redont P., Soubeyran A., Proximal alternating minimization and projection methods for nonconvex problems: an approach based on the Kurdyka–Łojasiewicz inequality, Math. Oper. Res., 2010, 35, 438-457.
  * [14] Nikolova M., Tan P., Alternating structure-adapted proximal gradient descent for nonconvex block-regularised problems, SIAM J. Optim., 2019, 29(3), 2053-2078.
  * [15] Ochs P., Chen Y., Brox T., and Pock T., iPiano: Inertial proximal algorithm for nonconvex optimization, SIAM J. Imaging Sci., 2014, 7(2), 1388-1419.
  * [16] Boţ R.I., Csetnek E.R., An inertial Tseng’s type proximal algorithm for nonsmooth and nonconvex optimization problems, J. Optim. Theory Appl., 2016, 171(2), 600-616.
  * [17] Polyak B.T., Some methods of speeding up the convergence of iteration methods, USSR Comput. Math. Math. Phys., 1964, 4, 1-17.
  * [18] Le T.K.H., Nicolas G., Panagiotis P., Inertial block proximal methods for non-convex non-smooth optimization, Proceedings of the 37th International Conference on Machine Learning, PMLR, 2020, 119, 5671-5681.
  * [19] Feng J.K., Zhang H.B., Zhang K.L., Zhao P.F., An inertial Douglas–Rachford splitting algorithm for nonconvex and nonsmooth problems, Concurrency and Computation Practice and Experience, 2021, https://doi.org/10.1002/cpe.6343.
  * [20] Zhang Y.X., He S.N., Inertial proximal alternating minimization for nonconvex and nonsmooth problems, J. Inequal. Appl., 2017, 2017, 232.
  * [21] Pock T., Sabach S., Inertial proximal alternating linearized minimization (iPALM) for nonconvex and nonsmooth problems, SIAM J. Imaging Sci., 2017, 9, 1756-1787.
  * [22] Gao X., Cai X.J., Han D.R., A Gauss-Seidel type inertial proximal alternating linearized minimization for a class of nonconvex optimization problems, J. Glob. Optim., 2020, 76, 863-887.
  * [23] Wang Q. X., Han D. R., A generalized inertial proximal alternating linearized minimization method for nonconvex nonsmooth problems, Appl. Numer. Math., 2023, 189, 66-87.
  * [24] Yang X., Xu L.L., Some accelerated alternating proximal gradient algorithms for a class of nonconvex nonsmooth problems, J. Glob. Optim., 2022, https://doi.org/10.1007/s10898-022-01214-3.
  * [25] Zhao J., Dong Q.L., Michael Th.R., Wang F.H., Two-step inertial Bregman alternating minimization algorithm for nonconvex and nonsmooth problems, J. Glob. Optim., 2022, 84, 941-966.
  * [26] Chao M.T., Nong F.F., Zhao M.Y., An inertial alternating minimization with Bregman distance for a class of nonconvex and nonsmooth problems, J. Appl. Math. Comput., 2023, 69, 1559-1581.
  * [27] Mordukhovich B., Variational Analysis and Generalized Differentiation, I: Basic Theory. Grundlehren der Mathematischen Wissenschaften, Vol. 330. Springer-Verlag, Berlin, 2006.
  * [28] Bolte J., Daniilidis A., Lewis A., The Łojasiewicz inequality for nonsmooth subanalytic functions with applications to subgradient dynamical systems, SIAM J. Optim., 2006, 17(4), 1205-1223.
  * [29] Rockafellar R.T., Wets R., Variational analysis, Grundlehren der Mathematischen Wissenschaften, Vol. 317. Springer, Berlin, 1998.
  * [30] Wang F., Cao W., Xu Z., Convergence of multi-block Bregman ADMM for nonconvex composite problems, Sci. China Inf. Sci., 2018, 61, 122101.
  * [31] Xu Y., Yin W., A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion, SIAM J. Imaging Sci., 2013, 6, 1758-1789.
  * [32] Hsieh Y.P., Kao Y.C., Mahabadi R.K., Yurtsever A., Kyrillidis A., Cevher V., A non-Euclidean gradient descent framework for nonconvex matrix factorization, IEEE Trans. Signal Process., 2018, 66, 5917-5926.
  * [33] Lee D.D, Seung H.S., Learning the parts of objects by non-negative matrix factorization, Nature, 1999, 788-791.
  * [34] Pan J., Gillis N., Generalized separable nonnegative matrix factorization, IEEE Trans. Pattern Anal. Mach. Intell., 2021, 43(5), 1546-1561.
  * [35] Rousset F., Peyrin F., Ducros N., A semi nonnegative matrix factorization technique for pattern generalization in single-pixel imaging, IEEE Trans. Comput. Imaging, 2018, 4(2), 284-294.
  * [36] Peharz R., Pernkopf F., Sparse nonnegative matrix factorization with $l_{0}$-constraints, Neurocomputing, 2012, 80, 38-46.
  * [37] Xu Z.B., Chang X.Y., Xu F.M., Zhang H.: $L_{1/2}$ regularization: a thresholding representation theory and a fast solver, IEEE Trans. Neural Netw. Learning Syst., 2012, 23, 1013-1027.