Inertial Proximal Alternating Linearized Minimization (iPALM) for Nonconvex and Nonsmooth Problems

TL;DR

This paper introduces an inertial version of the PALM algorithm for nonconvex, nonsmooth optimization problems with semi-algebraic data, proving its global convergence and demonstrating its effectiveness through numerical experiments.

Contribution

It develops and analyzes an inertial PALM algorithm with proven global convergence for complex nonconvex, nonsmooth problems, extending existing methods.

Findings

Proves global convergence of inertial PALM to critical points.

Demonstrates effectiveness in blind image deconvolution.

Shows improved performance in sparse matrix factorization.

Abstract

In this paper we study nonconvex and nonsmooth optimization problems with semi-algebraic data, where the variables vector is split into several blocks of variables. The problem consists of one smooth function of the entire variables vector and the sum of nonsmooth functions for each block separately. We analyze an inertial version of the Proximal Alternating Linearized Minimization (PALM) algorithm and prove its global convergence to a critical point of the objective function at hand. We illustrate our theoretical findings by presenting numerical experiments on blind image deconvolution, on sparse non-negative matrix factorization and on dictionary learning, which demonstrate the viability and effectiveness of the proposed method.

Tables

Table 1. Table 1: Values of the objective function of the sparse NMF problem after K 𝐾 K iterations, using different settings for the inertial parameters α i subscript 𝛼 𝑖 \alpha_{i} and β i subscript 𝛽 𝑖 \beta_{i} , i = 1 , 2 𝑖 1 2 i=1,2 and using computation of the exact Lipschitz constant.

	$K=100$	$K=500$	$K=1000$	$K=5000$	time (s)
${\alpha }_{1,2}={\beta }_{1,2}=0.0$	12968.17	7297.70	5640.11	4088.22	196.63
${\alpha }_{1,2}={\beta }_{1,2}=0.2$	12096.91	8453.29	6810.63	4482.00	190.47
${\alpha }_{1,2}={\beta }_{1,2}=0.4$	12342.81	11496.57	9277.11	5617.02	189.27
${\alpha }_1={\beta }_1=0.2,{\alpha }_2={\beta }_2=0.4$	12111.55	8488.35	6822.95	4465.59	201.05
${\alpha }_1={\beta }_1=0.4,{\alpha }_2={\beta }_2=0.8$	12358.00	11576.72	9350.37	5593.84	200.11
${\alpha }_{1,2}^k={\beta }_{1,2}^k=\left(k-1\right)/\left(k+2\right)$	5768.63	3877.41	3870.98	3870.81	186.62

Table 2. Table 2: Values of the objective function of the sparse NMF problem after K 𝐾 K iterations, using different settings for the inertial parameters α i subscript 𝛼 𝑖 \alpha_{i} and β i subscript 𝛽 𝑖 \beta_{i} , i = 1 , 2 𝑖 1 2 i=1,2 and using backtracking to estimate the Lipschitz constant.

	$K=100$	$K=500$	$K=1000$	$K=5000$	time (s)
${\alpha }_{1,2}={\beta }_{1,2}=0.0$	8926.23	5037.89	4356.65	4005.53	347.17
${\alpha }_{1,2}={\beta }_{1,2}=0.2$	8192.71	4776.64	4181.40	4000.41	349.42
${\alpha }_{1,2}={\beta }_{1,2}=0.4$	8667.62	4696.64	4249.57	4060.95	351.78
${\alpha }_1={\beta }_1=0.2,{\alpha }_2={\beta }_2=0.4$	8078.14	4860.74	4274.46	3951.28	353.53
${\alpha }_1={\beta }_1=0.4,{\alpha }_2={\beta }_2=0.8$	8269.27	4733.76	4243.29	4066.63	357.35
${\alpha }_{1,2}^k={\beta }_{1,2}^k=\left(k-1\right)/\left(k+2\right)$	5071.71	3902.91	3896.40	3869.13	347.90
iPiano ($\beta =0.4$)	14564.65	12200.78	11910.65	7116.22	258.42

Table 3. Table 3: Values of the objective function of the BID problem, after K 𝐾 K iterations and using different settings for the inertial parameters α i subscript 𝛼 𝑖 \alpha_{i} and β i subscript 𝛽 𝑖 \beta_{i} for i = 1 , 2 𝑖 1 2 i=1,2 .

	$K=100$	$K=500$	$K=1000$	$K=5000$	time (s)
${\alpha }_{1,2}={\beta }_{1,2}=0.0$	2969668.92	1177462.72	1031575.57	847268.70	1882.63
${\alpha }_{1,2}={\beta }_{1,2}=0.4$	5335748.90	1402080.44	1160510.16	719295.30	1895.61
${\alpha }_{1,2}={\beta }_{1,2}=0.8$	5950073.38	1921105.31	1447739.06	780109.56	1888.25
${\alpha }_{1,2}^k={\beta }_{1,2}^k=\left(k-1\right)/\left(k+2\right)$	2014059.03	978234.23	683694.72	678090.51	1867.19

Table 4. Table 4: Values of the objective function for the convolutional LASSO model using different settings of the inertial parameters.

	$K=100$	$K=200$	$K=500$	$K=1000$	time (s)
${\alpha }_{1,2}={\beta }_{1,2}=0.0$	336.13	328.21	322.91	321.12	3274.97
${\alpha }_{1,2}={\beta }_{1,2}=0.4$	329.20	324.62	321.51	319.85	3185.04
${\alpha }_{1,2}={\beta }_{1,2}=0.8$	325.19	321.38	319.79	319.54	3137.09
${\alpha }_{1,2}^k={\beta }_{1,2}^k=\left(k-1\right)/\left(k+2\right)$	323.23	319.88	318.64	318.44	3325.37

Equations

$$\underset{\mathbf{x}=\left(x_1,x_2\right)}{\min }F\left(\mathbf{x}\right):=f_1\left(x_1\right)+f_2\left(x_2\right)+H\left(\mathbf{x}\right),$$

$$\underset{B,C}{\min }\frac{1}{2}{\Vert A-BC\Vert }_F^2,\text{s.t.}B\ge 0,C\ge 0,$$

$$A=B\ast K+E,$$

$$\underset{U,K}{\min }\mathcal{R}\left(U\right)+\frac{1}{2}{\Vert B\ast K-A\Vert }_F^2,\text{s.t.}0\le U\le 1,K\in {\Delta }_{mn}.$$

$$x^{k+1}=x^k-\tau \nabla f\left(x^k\right)+\beta \left(x^k-x^{k-1}\right),$$

$$\ddot{x}\left(t\right)+c\dot{x}\left(t\right)+g\left(x\left(t\right)\right)=0,$$

$$x^{k+1}=x^k-\sum _{i=0}^k{\alpha }_i\nabla f\left(x^i\right).$$

$$\text{mbox}minimizeF\left(\mathbf{x}\right):=f_1\left(x_1\right)+f_2\left(x_2\right)+H\left(\mathbf{x}\right)\text{mbox}overall\mathbf{x}=\left(x_1,x_2\right)\in {\mathbb{R}}^{n_1}\times {\mathbb{R}}^{n_2},$$

$${\mathrm{prox}}_t^{\sigma }\left(p\right):=\mathrm{argmin}\left\{\sigma \left(q\right)+\frac{t}{2}{\Vert q-p\Vert }^2:q\in {\mathbb{R}}^d\right\},\left(t>0\right).$$

$$\Vert {\nabla }_{x_1}H\left(u,x_2\right)-{\nabla }_{x_1}H\left(v,x_2\right)\Vert \le L_1\left(x_2\right)\Vert u-v\Vert ,\forall u,v\in {\mathbb{R}}^{n_1}.$$

$$\inf \left\{L_1\left(x_2\right):x_2\in B_2\right\}$$

$$\sup \left\{L_1\left(x_2\right):x_2\in B_2\right\}$$

$$\Vert \left({\nabla }_{x_1}H\left(x_1,x_2\right)-{\nabla }_{x_1}H\left(y_1,y_2\right),{\nabla }_{x_2}H\left(x_1,x_2\right)-{\nabla }_{x_2}H\left(y_1,y_2\right)\right)\Vert$$

$$\le M\Vert \left(x_1-y_1,x_2-y_2\right)\Vert .$$

$$y_1^k$$

$$z_1^k$$

$$x_1^{k+1}$$

$$y_2^k$$

$$z_2^k$$

$$x_2^{k+1}$$

$$\mathrm{crit}\sigma =\left\{u\in \mathrm{dom}\sigma :0\in \partial \sigma \left(u\right)\right\}.$$

$${\Phi }_{\eta }\equiv \left\{\phi \in C\left[\left[0,\eta \right),{\mathbb{R}}_{+}\right]\text{mbox}suchthat\{\begin{array}{ll}\phi \left(0\right)=0& \\ \phi \in C^1& \text{mbox}on\left(0,\eta \right)\\ {\phi }^{\prime }\left(s\right)>0& \text{mbox}foralls\in \left(0,\eta \right)\end{array}\}\right\}.$$

$$u\in U\cap [\sigma (\overline{u})<\sigma (u)<\sigma (\overline{u})+\eta ],$$

$$\phi \left(\sigma \left(u\right)-\sigma \left(\overline{u}\right)\right)\mathrm{dist}\left(0,\partial \sigma \left(u\right)\right)\ge 1,$$

$$\mathrm{dist}\left(x,S\right):=\inf \left\{\Vert y-x\Vert :y\in S\right\}.$$

$$\left\{u\in {\mathbb{R}}^d:\mathrm{dist}\left(u,\Omega \right)<\epsilon \right\}\cap \left[\sigma \left(\overline{u}\right)<\sigma \left(u\right)<\sigma \left(\overline{u}\right)+\eta \right],$$

$${\phi }^{\prime }\left(\sigma \left(u\right)-\sigma \left(\overline{u}\right)\right)\mathrm{dist}\left(0,\partial \sigma \left(u\right)\right)\ge 1.$$

$$u^{+}\in {\mathrm{prox}}_t^{\sigma }\left(v-\frac{1}{t}\nabla h\left(w\right)\right),t>0,$$

$$g\left(u^{+}\right)\le g\left(u\right)+\frac{L_h+s}{2}{\Vert u^{+}-u\Vert }^2+\frac{t}{2}{\Vert u-v\Vert }^2-\frac{t}{2}{\Vert u^{+}-v\Vert }^2+\frac{L_h^2}{2s}{\Vert u-w\Vert }^2,$$

$$u^{+}\in {\mathrm{argmin}}_{\xi \in {\mathbb{R}}^d}\left\{⟨\xi -v,\nabla h\left(w\right)⟩+\frac{t}{2}{\Vert \xi -v\Vert }^2+\sigma \left(\xi \right)\right\},$$

$$⟨u^{+}-v,\nabla h\left(w\right)⟩+\frac{t}{2}{\Vert u^{+}-v\Vert }^2+\sigma \left(u^{+}\right)\le ⟨u-v,\nabla h\left(w\right)⟩+\frac{t}{2}{\Vert u-v\Vert }^2+\sigma \left(u\right).$$

$$\sigma \left(u^{+}\right)\le ⟨u-u^{+},\nabla h\left(w\right)⟩+\frac{t}{2}{\Vert u-v\Vert }^2-\frac{t}{2}{\Vert u^{+}-v\Vert }^2+\sigma \left(u\right).$$

$$h\left(u^{+}\right)+\sigma \left(u^{+}\right)$$

Full text

Key words: Alternating minimization, blind image deconvolution, block coordinate descent, heavy-ball method, Kurdyka-Łojasiewicz property, nonconvex and nonsmooth minimization, sparse non-negative matrix factorization, dictionary learning.

1 Introduction

In the last decades advances in convex optimization have significantly influenced scientific fields such as image processing and machine learning, which are dominated by computational approaches. However, it is also known that the framework of convexity is often too restrictive to provide good models for many practical problems. Several basic problems such as blind image deconvolution are inherently nonconvex and hence there is a vital interest in the development of efficient and simple algorithms for tackling nonconvex optimization problems.

A large part of the optimization community has also been devoted to the development of general purpose solvers [23], but in the age of big data, such algorithms often come to their limits since they cannot efficiently exploit the structure of the problem at hand. One notable exception is the general purpose limited quasi Newton method [18] which has been published more than 25 years ago but still remains a competitive method.

A promising approach to tackle nonconvex problems is to consider a very rich the class of problems which share certain structure that allows the development of efficient algorithms. One such class of nonconvex optimization problems is given by the sum of three functions:

[TABLE]

where $f_{1}$ and $f_{2}$ are assumed to be general nonsmooth and nonconvex functions with efficiently computable proximal mappings (see exact definition in the next section) and $H$ is a smooth coupling function which is required to have only partial Lipschitz continuous gradients $\nabla_{x_{1}}H$ and $\nabla_{x_{2}}H$ (it should be noted that $\nabla_{{\bf x}}H$ might not be a Lipschitz continuous).

Many practical problems frequently used in the machine learning and image processing communities fall into this class of problems. Let us briefly mention two classical examples (another example will be discussed in Section 5).

The first example is Non-negative Matrix Factorization (NMF) [26, 15]. Given a non-negative data matrix $A\in\mathbb{R}_{+}^{m\times n}$ and an integer $r>0$, the idea is to approximate the matrix $A$ by a product of again non-negative matrices $BC$, where $B\in\mathbb{R}_{+}^{m\times r}$ and $C\in\mathbb{R}_{+}^{r\times n}$. It should be noted that the dimension $r$ is usually much smaller than $\min\{m,n\}$. Clearly this problem is very difficult to solve and hence several algorithms have been developed (see, for example, [33]). One possibility to solve this problem is by finding a solution for the non-negative least squares model given by

[TABLE]

where the non-negativity constraint is understood pointwise and $\left\|{\cdot}\right\|_{F}$ denotes the classical Frobenius norm. The NMF has important applications in image processing (face recognition) and bioinformatics (clustering of gene expressions). Observe that the gradient of the objective function is not Lipschitz continuous but it is partially Lipschitz continuous which enables the application of alternating minimization based methods (see [10]). Additionally, it is popular to impose sparsity constraints on one or both of the unknowns, e.g., $\left\|{C}\right\|_{0}\leq c$, to promote sparsity in the representation. See [10] for the first globally convergent algorithm for solving the sparse NMF problem. As we will see, the complicated sparse NMF can be also simply handled by our proposed algorithm which seems to produce better performances (see Section 5).

The second example we would like to mention is the important but ever challenging problem of blind image deconvolution (BID) [16]. Let $A\in\left[0,1\right]^{M\times N}$ be the observed blurred image of size $M\times N$, and let $B\in\left[0,1\right]^{M\times N}$ be the unknown sharp image of the same size. Furthermore, let $K\in\Delta_{mn}$ denote a small unknown blur kernel (point spread function) of size $m\times n$, where $\Delta_{mn}$ denotes the $mn$-dimensional standard unit simplex. We further assume that the observed blurred image has been formed by the following linear image formation model:

[TABLE]

where $\ast$ denotes a two dimensional discrete convolution operation and $E$ denotes a small additive Gaussian noise. A typical variational formulation of the blind deconvolution problem is given by:

[TABLE]

In the above variational model, $\mathcal{R}$ is an image regularization term, typically a function, that imposes sparsity on the image gradient and hence favoring sharp images over blurred images.

We will come back to both examples in Section 5 where we will show how the proposed algorithm can be applied to efficiently solve these problems.

In [10], the authors proposed a proximal alternating linearized minimization method (PALM) that efficiently exploits the structure of problem (1.1). PALM can be understood as a blockwise application of the well-known proximal forward-backward algorithm [17, 11] in the nonconvex setting. In the case that the objective function $F$ satisfy the so-called Kurdyka-Łojasiewicz (KL) property (the exact definition will be given in Section 3), the whole sequence of the algorithm is guaranteed to converge to a critical point of the problem.

In this paper, we propose an inertial version of the PALM algorithm and show convergence of the whole sequence in case the objective function $F$ satisfy the KL property. The inertial term is motivated from the Heavy Ball method of Polyak [29] which in its most simple version applied to minimizing a smooth function $f$ and can be written as the iterative scheme

[TABLE]

where $\beta$ and $\tau$ are suitable parameters that ensure convergence of the algorithm. The heavy ball method differs from the usual gradient method by the additional inertial term $\beta\left(x^{k}-x^{k-1}\right)$, which adds part of the old direction to the new direction of the algorithm. Therefore for $\beta=0$, we completely recover the classical algorithm of unconstrained optimization, the Gradient Method. The heavy ball method can be motivated from basically three view points.

First, the heavy ball method can be seen as an explicit finite differences discretization of the heavy ball with friction dynamical system (see [2]):

[TABLE]

where $x\left(t\right)$ is a time continuous trajectory, $\ddot{x}\left(t\right)$ is the acceleration, $c\dot{x}\left(t\right)$ for $c>0$ is the friction (damping), which is proportional to the velocity $\dot{x}\left(t\right)$, and $g\left(x\left(t\right)\right)$ is an external gravitational field. In the case that $g=\nabla f$ the trajectory $x\left(t\right)$ is running down the “energy landscape” described by the objective function $f$ until a critical point ($\nabla f=0$) is reached. Due to the presence of the inertial term, it can also overcome spurious critical points of $f$, e.g., saddle points.

Second, the heavy ball method can be seen as a special case of the so-called multi-step algorithms where each step of the algorithm is given as a linear combination of all previously computed gradients [12], that is, algorithm of the following form

[TABLE]

Let us note that in the case that the objective function $f$ is quadratic, the parameters $\alpha_{i}$ can be chosen in a way such that the objective function is minimized at each step. This approach eventually leads to the Conjugate Gradient (CG) method, pointing out a close relationship to inertial based methods.

Third, accelerated gradient methods, as pioneered by Nesterov (see [21] for an overview), are based on a variant of the heavy ball method that use the extrapolated point (based on the inertial force) also for evaluating the gradient in the current step. It turns out that, in the convex setting, these methods improve the worst convergence rate from $\mathcal{O}(1/k)$ to $\mathcal{O}(1/k^{2})$, while leaving the computational complexity of each step basically the same.

In [35], the heavy ball method has been analyzed for the first time in the setting of nonconvex problems. It is shown that the heavy ball method is attracted by the connected components of critical points. The proof is based on considering a suitable Lyapunov function that allows to consider the two-step algorithm as a one-step algorithm. As we will see later, our convergence proof is also based on rewriting the algorithm as a one-step method.

In [24], the authors developed an inertial proximal gradient algorithm (iPiano). The algorithm falls into the class of forward-backward splitting algorithms [11], as it performs an explicit forward (steepest descent) step with respect to the smooth (nonconvex) function followed by a (proximal) backward step with respect to the nonsmooth (convex) function. Motivated by the heavy ball algorithm mentioned before, the iPiano algorithm makes use of an inertial force which empirically shows to improve the convergence speed of the algorithm. A related method based on general Bregman proximal-like distance functions has been recently proposed in [14].

Very recently, a randomized proximal linearization method has been proposed in [34]. The method is closely related to our proposed algorithm but convergence is proven only in the case that the function values are strictly decreasing. This is true only if the inertial force is set to be zero or the algorithm is restarted whenever the function values are not decreasing. In this paper, however, we overcome this major drawback and prove convergence of the algorithm without any assumption on the monotonicity of the objective function.

The remainder of the paper is organized as follows. In Section 2 we give an exact definition of the problem and the proposed algorithm. In Section 3 we state few technical results that will be necessary for the convergence analysis, which will be presented in Section 4. In Section 5 we present some numerical results and analyze the practical performance of the algorithm in dependence of its inertial parameters.

2 Problem Formulation and Algorithm

In this paper we follow [10] and consider the broad class of nonconvex and nonsmooth problems of the following form

[TABLE]

where $f_{1}$ and $f_{2}$ are extended valued (i.e., giving the possibility of imposing constraints separately on the blocks $x_{1}$ and $x_{2}$) and $H$ is a smooth coupling function (see below for more precise assumptions on the involved functions). We would like to stress from the beginning that even though all the discussions and results of this paper derived for two blocks of variables $x_{1}$ and $x_{2}$, they hold true for any finite number of blocks. This choice was done only for the sake of simplicity of the presentation of the algorithm and the convergence results.

As we discussed in the introduction, the proposed algorithm can be viewed either as a block version of the recent iPiano algorithm [24] or as an inertial based version of the recent PALM algorithm [10]. Before presenting the algorithm it will be convenient to recall the definition of the Moreau proximal mapping [20]. Given a proper and lower semicontinuous function $\sigma:\mathbb{R}^{d}\rightarrow\left(-\infty,\infty\right]$, the proximal mapping associated with $\sigma$ is defined by

[TABLE]

Following [10], we take the following as our blanket assumption.

Assumption A**.**

• $\rm{(i)}$

$f_{1}:\mathbb{R}^{n_{1}}\rightarrow\left(-\infty,\infty\right]$ and $f_{2}:\mathbb{R}^{n_{2}}\rightarrow\left(-\infty,\infty\right]$ are proper and lower semicontinuous functions such that $\inf_{\mathbb{R}^{n_{1}}}f_{1}>-\infty$ and $\inf_{\mathbb{R}^{n_{2}}}f_{2}>-\infty$.

• $\rm{(ii)}$

$H:\mathbb{R}^{n_{1}}\times\mathbb{R}^{n_{2}}\rightarrow\mathbb{R}$ is differentiable and $\inf_{\mathbb{R}^{n_{1}}\times\mathbb{R}^{n_{2}}}F>-\infty$.

• $\rm{(iii)}$

For any fixed $x_{2}$ the function $x_{1}\rightarrow H\left(x_{1},x_{2}\right)$ is $C^{1,1}_{L_{1}(x_{2})}$, namely the partial gradient $\nabla_{x_{1}}H\left(x_{1},x_{2}\right)$ is globally Lipschitz with moduli $L_{1}\left(x_{2}\right)$, that is,

[TABLE]

Likewise, for any fixed $x_{1}$ the function $x_{2}\rightarrow H\left(x_{1},x_{2}\right)$ is assumed to be $C^{1,1}_{L_{2}(x_{1})}$.

• $\rm{(iv)}$

For $i=1,2$ there exists $\lambda_{i}^{-},\lambda_{i}^{+}>0$ such that

[TABLE]

for any compact set $B_{i}\subseteq\mathbb{R}^{n_{i}}$, $i=1,2$.

• $\rm{(v)}$

$\nabla H$ is Lipschitz continuous on bounded subsets of $\mathbb{R}^{n_{1}}\times\mathbb{R}^{n_{2}}$. In other words, for each bounded subset $B_{1}\times B_{2}$ of $\mathbb{R}^{n_{1}}\times\mathbb{R}^{n_{2}}$ there exists $M>0$ such that:

[TABLE]

We propose now the inertial Proximal Alternating Linearized Minimization (iPALM) algorithm.

iPALM: Inertial Proximal Alternating Linearized Minimization

Initialization: start with any $\left(x_{1}^{0},x_{2}^{0}\right)\in\mathbb{R}^{n_{1}}\times\mathbb{R}^{n_{2}}$.

For each $k=1,2,\ldots$ generate a sequence $\left\{\left(x_{1}^{k},x_{2}^{k}\right)\right\}_{k\in\mathbb{N}}$ as follows:

2.1.

Take $\alpha_{1}^{k},\beta_{1}^{k}\in\left[0,1\right]$ and $\tau_{1}^{k}>0$. Compute

$\displaystyle y_{1}^{k}$ $\displaystyle=x_{1}^{k}+\alpha_{1}^{k}\left(x_{1}^{k}-x_{1}^{k-1}\right),$

(2.5)

$\displaystyle z_{1}^{k}$ $\displaystyle=x_{1}^{k}+\beta_{1}^{k}\left(x_{1}^{k}-x_{1}^{k-1}\right),$

(2.6)

$\displaystyle x_{1}^{k+1}$ $\displaystyle\in\mathrm{prox}_{\tau_{1}^{k}}^{f_{1}}\left(y_{1}^{k}-\frac{1}{\tau_{1}^{k}}\nabla_{x_{1}}H\left(z_{1}^{k},x_{2}^{k}\right)\right).$

(2.7)

2.2.

Take $\alpha_{2}^{k},\beta_{2}^{k}\in\left[0,1\right]$ and $\tau_{2}^{k}>0$. Compute

$\displaystyle y_{2}^{k}$ $\displaystyle=x_{2}^{k}+\alpha_{2}^{k}\left(x_{2}^{k}-x_{2}^{k-1}\right),$

(2.8)

$\displaystyle z_{2}^{k}$ $\displaystyle=x_{2}^{k}+\beta_{2}^{k}\left(x_{2}^{k}-x_{2}^{k-1}\right),$

(2.9)

$\displaystyle x_{2}^{k+1}$ $\displaystyle\in\mathrm{prox}_{\tau_{2}^{k}}^{f_{2}}\left(y_{2}^{k}-\frac{1}{\tau_{2}^{k}}\nabla_{x_{2}}H\left(x_{1}^{k+1},z_{2}^{k}\right)\right).$

(2.10)

The parameters $\tau_{1}^{k}$ and $\tau_{2}^{k}$, $k\in\mathbb{N}$, are discussed in Section 4 but for now, we can say that they are proportional to the respective partial Lipschitz moduli of $H$. The larger the partial Lipschitz moduli the smaller the step-size, and hence the slower the algorithm. As we shall see below, the partial Lipschitz moduli $L_{1}\left(x_{2}\right)$ and $L_{2}\left(x_{1}\right)$ are explicitly available for the examples mentioned at the introduction. However, note that if these are unknown, or still too difficult to compute, then a backtracking scheme [6] can be incorporated and the convergence results developed below remain true, for simplicity of exposition we omit the details.

In Section 5, we will show that the involved functions of the non-negative matrix factorization model (see (1.2)) and of the blind image deconvulation model (see (1.3)) do satisfy Assumption A. For the general setting we point out the following remarks about Assumption A.

• $\rm{(i)}$

The first item of Assumption A is very general and most of the interesting constraints (via their indicator functions) or regularizing functions fulfill these requirements.

• $\rm{(ii)}$

Items (ii)-(v) of Assumption A are beneficially exploited to build the proposed iPALM algorithm. These requirements do not guarantee that the gradient of $H$ is globally Lipschitz (which is the case in our mentioned applications). The fact that $\nabla H$ is not globally Lipschitz reduces the potential of applying the iPiano and PFB methods in concrete applications and therefore highly motivated us to study their block counterparts (PALM in [10] and iPALM in this paper).

• $\rm{(iii)}$

Another advantage of algorithms that exploit block structures inherent in the model at hand is the fact that they achieve better numerical performance (see Section 5) by taking step-sizes which is optimized to each separated block of variables.

• $\rm{(iv)}$

Item (v) of Assumption A holds true, for example, when $H$ is $C^{2}$. In this case the inequalities in (2.4) could be obtained if the sequence, which generated by the algorithm, is bounded.

The iPALM algorithm generalizes few known algorithms for different values of the inertial parameters $\alpha_{i}^{k}$ and $\beta_{i}^{k}$, $k\in\mathbb{N}$ and $i=1,2$. For example, when $\alpha_{i}^{k}=\beta_{i}^{k}=0$, $k\in\mathbb{N}$, we recover the PALM algorithm of [10] which is a block version of the classical Proximal Forward-Backward (PFB) algorithm. When, there is only one block of variables, for instance only $i=1$, we get the iPiano algorithm [24] which is recovered exactly only when $\beta_{1}^{k}=0$, $k\in\mathbb{N}$. It should be also noted that in [24], the authors additionally assume that the function $f_{1}$ is convex (an assumption that is not needed in our case). The iPiano algorithm by itself generalizes two classical and known algorithms, one is the Heavy-Ball method [30] (when $f_{1}\equiv 0$) and again the PFB method (when $\alpha_{1}^{k}=0$, $k\in\mathbb{N}$).

3 Mathematical Preliminaries and Proof Methodology

Throughout this paper we are using standard notations and definitions of nonsmooth analysis which can be found in any classical book, see for instance [31, 19]. We recall here few notations and technical results. Let $\sigma:\mathbb{R}^{d}\rightarrow\left(-\infty,\infty\right]$ be a proper and lower semicontinuous function. Since we are dealing with nonconvex and nonsmooth functions that can have the value $\infty$, we use the notion of limiting subdifferential (or simply subdifferential), see [19], which is denoted by $\partial\sigma$. In what follows, we are interested in finding critical points of the objective function $F$ defined in (2.1). Critical points are those points for which the corresponding subdifferential contains the zero vector ${\bf 0}$. The set of critical points of $\sigma$ is denoted by $\mathrm{crit}\,{\sigma}$, that is,

[TABLE]

An important property of the subdifferential is recorded in the following remark (see [31]).

Remark 3.1**.**

Let $\left\{\left(u^{k},q^{k}\right)\right\}_{k\in\mathbb{N}}$ be a sequence in $\mathrm{graph}\,{\left(\partial\sigma\right)}$ that converges to $\left(u,q\right)$ as $k\rightarrow\infty$. By the definition of $\partial\sigma\left(u\right)$, if $\sigma\left(u^{k}\right)$ converges to $\sigma\left(u\right)$ as $k\rightarrow\infty$, then $\left(u,q\right)\in\mathrm{graph}\,{(\partial\sigma)}$.

The convergence analysis of iPALM is based on the proof methodology which was developed in [5] and more recently extended and simplified in [10]. The main part of the suggested methodology relies on the fact that the objective function of the problem at hand satisfies the Kurdyka-Łojasiewicz (KL) property. Before stating the KL property we will need the definition of the following class of desingularizing functions. For $\eta\in\left(0,\infty\right]$ define

[TABLE]

The function $\sigma$ is said to have the Kurdyka-Łojasiewicz (KL) property at $\overline{u}\in\mathrm{dom}\,{\partial}\sigma$ if there exist $\eta\in\left(0,\infty\right]$, a neighborhood $U$ of $\overline{u}$ and a function $\varphi\in\Phi_{\eta}$, such that, for all

[TABLE]

the following inequality holds

[TABLE]

where for any subset $S\subset\mathbb{R}^{d}$ and any point $x\in\mathbb{R}^{d}$

[TABLE]

When $S=\emptyset$, we have that $\mathrm{dist}\left(x,S\right)=\infty$ for all $x$. If $\sigma$ satisfies property (3.2) at each point of $\mathrm{dom}\,{\partial\sigma}$, then $\sigma$ is called a KL function.

The convergence analysis presented in the following section is based on the uniformized KL property which was established in [10, Lemma 6, p. 478].

Lemma 3.1**.**

Let $\Omega$ be a compact set and let $\sigma:\mathbb{R}^{d}\rightarrow\left(-\infty,\infty\right]$ be a proper and lower semicontinuous function. Assume that $\sigma$ is constant on $\Omega$ and satisfies the KL property at each point of $\Omega$. Then, there exist $\varepsilon>0$, $\eta>0$ and $\varphi\in\Phi_{\eta}$ such that for all $\overline{u}$ in $\Omega$ and all $u$ in the following intersection

[TABLE]

one has,

[TABLE]

We refer the reader to [9] for a depth study of the class of KL functions. For the important relation between semi-algebraic and KL functions see [8]. In [3, 4, 5, 10], the interested reader can find through catalog of functions which are very common in many applications and satisfy the KL property.

Before concluding the mathematical preliminaries part we would like to mention few important properties of the proximal map (defined in (2.2)). The following result can be found in [31].

Proposition 3.1**.**

Let $\sigma:\mathbb{R}^{d}\rightarrow\left(-\infty,\infty\right]$ be a proper and lower semicontinuous function with $\inf_{\mathbb{R}^{d}}\sigma>-\infty$. Then, for every $t\in\left(0,\infty\right)$ the set $\mathrm{prox}_{t\sigma}\left(u\right)$ is nonempty and compact.

It follows immediately from the definition that $\mathrm{prox}{\sigma}$ is a multi-valued map when $\sigma$ is nonconvex. The multi-valued projection onto a nonempty and closed set $C$ is recovered when $\sigma=\delta_{C}$, which is the indicator function of $C$ that defined to be zero on $C$ and $\infty$ outside.

The main computational effort of iPALM involves a proximal mapping step of a proper and lower semicontinuous but nonconvex function. The following property will be essential in the forthcoming convergence analysis and is a slight modification of [10, Lemma 2, p. 471].

Lemma 3.2** (Proximal inequality).**

Let $h:\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a continuously differentiable function with gradient $\nabla h$ assumed $L_{h}$-Lipschitz continuous and let $\sigma:\mathbb{R}^{d}\rightarrow\left(-\infty,\infty\right]$ be a proper and lower semicontinuous function with $\inf_{\mathbb{R}^{d}}\sigma>-\infty$. Then, for any $v,w\in\mathrm{dom}\,{\sigma}$ and any $u^{+}\in\mathbb{R}^{d}$ defined by

[TABLE]

we have, for any $u\in\mathrm{dom}\,{\sigma}$ and any $s>0$:

[TABLE]

where $g:=h+\sigma$.

Proof.

First, it follows immediately from Proposition 3.1 that $u^{+}$ is well-defined. By the definition of the proximal mapping (see (2.2)) we get that

[TABLE]

and hence in particular, by taking $\xi=u$, we obtain

[TABLE]

Thus

[TABLE]

Invoking first the descent lemma (see [7]) for $h$, and using (3.7), yields

[TABLE]

Now, using the fact that $\left\langle{p,q}\right\rangle\leq\left(s/2\right)\left\|{p}\right\|^{2}+\left(1/2s\right)\left\|{q}\right\|^{2}$ for any two vectors $p,q\in\mathbb{R}^{d}$ and every $s>0$, yields

[TABLE]

where we have used the fact that $\nabla h$ is $L_{h}$-Lipschitz continuous. Thus, combining the last two inequalities proves that (3.6) holds. ∎

Remark 3.2**.**

It should be noted that if the nonsmooth function $\sigma$ is also known to be convex, then we can derive the following tighter upper bound (*cf. *(3.6))

[TABLE]

3.1 Convergence Proof Methodology

In this section we briefly summarize (*cf. *Theorem 3.1 below) the methodology recently proposed in [10] which provides the key elements to obtain an abstract convergence result that can be applied to any algorithm and will be applied here to prove convergence of iPALM. Let $\left\{{{\bf u}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ be a sequence in $\mathbb{R}^{d}$ which was generated from a starting point ${\bf u}^{0}$ by a generic algorithm ${\cal A}$. The set of all limit points of $\left\{{{\bf u}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ is denoted by $\omega\left({\bf u}^{0}\right)$, and defined by

[TABLE]

Theorem 3.1**.**

Let $\Psi:\mathbb{R}^{d}\rightarrow\left(-\infty,\infty\right]$ be a proper, lower semicontinuous and semi-algebraic function with $\inf\Psi>-\infty$. Assume that $\left\{{{\bf u}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ is a bounded sequence generated by a generic algorithm ${\cal A}$ from a starting point ${\bf u}^{0}$, for which the following three conditions hold true for any $k\in\mathbb{N}$.

• $\rm{(C1)}$

There exists a positive scalar $\rho_{1}$ such that

[TABLE]

• $\rm{(C2)}$

There exists a positive scalar $\rho_{2}$ such that for some ${\bf w}^{k}\in\partial\Psi\left({\bf u}^{k}\right)$ we have

[TABLE]

• $\rm{(C3)}$

Each limit point in the set $\omega\left({\bf u}^{0}\right)$ is a critical point of $\Psi$, that is, $\omega\left({\bf u}^{0}\right)\subset\mathrm{crit}\,{\Psi}$.

Then, the sequence $\left\{{{\bf u}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ converges to a critical point ${\bf u}^{\ast}$ of $\Psi$.

4 Convergence Analysis of iPALM

Our aim in this section is to prove that the sequence $\left\{\left(x_{1}^{k},x_{2}^{k}\right)\right\}_{k\in\mathbb{N}}$ which is generated by iPALM converges to a critical point of the objective function $F$ defined in (2.1). To this end we will follow the proof methodology described above in Theorem 3.1. In the case of iPALM, similarly to the iPiano algorithm (see [24]), it is not possible to prove that condition (C1) hold true for the sequence $\left\{\left(x_{1}^{k},x_{2}^{k}\right)\right\}_{k\in\mathbb{N}}$ and the function $F$, namely, this is not a descent algorithm with respect to $F$. Therefore, we first show that conditions (C1), (C2) and (C3) hold true for an auxiliary sequence and auxiliary function (see details below). Then, based on these properties we will show that the original sequence converges to a critical point of the original function $F$.

We first introduce the following notations that simplify the coming expositions. For any $k\in\mathbb{N}$, we define

[TABLE]

it is clear, that using these notations, we have that $\Delta^{k}=\Delta_{1}^{k}+\Delta_{2}^{k}$ for all $k\in\mathbb{N}$. Using these notations we can easily show few basic relations of the sequences $\left\{x_{i}^{k}\right\}_{k\in\mathbb{N}}$, $\left\{y_{i}^{k}\right\}_{k\in\mathbb{N}}$, and $\left\{z_{i}^{k}\right\}_{k\in\mathbb{N}}$, for $i=1,2$, generated by iPALM.

Proposition 4.1**.**

Let $\left\{\left(x_{1}^{k},x_{2}^{k}\right)\right\}_{k\in\mathbb{N}}$ be a sequence generated by iPALM. Then, for any $k\in\mathbb{N}$ and $i=1,2$, we have

• $\rm{(i)}$

$\left\|{x_{i}^{k}-y_{i}^{k}}\right\|^{2}=2\left(\alpha_{i}^{k}\right)^{2}\Delta_{i}^{k}$;

• $\rm{(ii)}$

$\left\|{x_{i}^{k}-z_{i}^{k}}\right\|^{2}=2\left(\beta_{i}^{k}\right)^{2}\Delta_{i}^{k}$;

• $\rm{(iii)}$

$\left\|{x_{i}^{k+1}-y_{i}^{k}}\right\|^{2}\geq 2\left(1-\alpha_{i}^{k}\right)\Delta_{i}^{k+1}+2\alpha_{i}^{k}\left(\alpha_{i}^{k}-1\right)\Delta_{i}^{k}$.

Proof.

The first two items follow immediately from the facts that $x_{i}^{k}-y_{i}^{k}=\alpha_{i}^{k}\left(x_{i}^{k-1}-x_{i}^{k}\right)$ and $x_{i}^{k}-z_{i}^{k}=\beta_{i}^{k}\left(x_{i}^{k-1}-x_{i}^{k}\right)$, for $i=1,2$ (see steps (2.5), (2.6), (2.8) and (2.9)). The last item follows from the following argument

[TABLE]

where we have used the fact that

[TABLE]

that follows from the Cauchy-Schwartz and Young inequalities. This proves item (iii). ∎

Now we prove the following property of the sequence $\left\{\left(x_{1}^{k},x_{2}^{k}\right)\right\}_{k\in\mathbb{N}}$ generated by iPALM.

Proposition 4.2**.**

Suppose that Assumption A holds. Let $\left\{\left(x_{1}^{k},x_{2}^{k}\right)\right\}_{k\in\mathbb{N}}$ be a sequence generated by iPALM, then for all $k\in\mathbb{N}$, we have that

[TABLE]

where $s_{1}^{k}>0$ and $s_{2}^{k}>0$ are arbitrarily chosen, for all $k\in\mathbb{N}$.

Proof.

Fix $k\geq 1$. Under our Assumption A(ii), the function $x_{1}\rightarrow H\left(x_{1},x_{2}\right)$ ($x_{2}$ is fixed) is differentiable and has a Lipschitz continuous gradient with moduli $L_{1}\left(x_{2}\right)$. Using the iterative step (2.7), applying Lemma 3.2 for $h\left(\cdot\right):=H\left(\cdot,x_{2}^{k}\right)$, $\sigma:=f_{1}$ and $t:=\tau_{1}^{k}$ with the points $u=x_{1}^{k}$, $u^{+}=x_{1}^{k+1}$, $v=y_{1}^{k}$ and $w=z_{1}^{k}$ yields that

[TABLE]

where the second inequality follows from Proposition 4.1. Repeating all the arguments above on the iterative step (2.10) yields the following

[TABLE]

By adding (4.3) and (4.4) we get

[TABLE]

This proves the desired result. ∎

Before we proceed and for the sake of simplicity of our developments we would like to chose the parameters $s_{1}^{k}$ and $s_{2}^{k}$ for all $k\in\mathbb{N}$. The best choice can be derived by minimizing the right-hand side of (3.6) with respect to $s$. Simple computations yields that the minimizer should be

[TABLE]

where $u,u^{+},w$ and $L_{h}$ are all in terms of Lemma 3.2. In Proposition 4.2 we have used Lemma 3.2 with the following choices $u=x_{1}^{k}$, $u^{+}=x_{1}^{k+1}$ and $w=z_{1}^{k}$. Thus

[TABLE]

where the last equality follows from step (2.6). Thus, from now on, we will use the following parameters:

[TABLE]

An immediate consequence of this choice of parameters which combined with Proposition 4.2 is recorded now.

Corollary 4.1**.**

Suppose that Assumption A holds. Let $\left\{\left(x_{1}^{k},x_{2}^{k}\right)\right\}_{k\in\mathbb{N}}$ be a sequence generated by iPALM, then for all $k\in\mathbb{N}$, we have that

[TABLE]

Similarly to iPiano, the iPALM algorithm generates a sequence which does not ensure that the function values decrease between two successive elements of the sequence. Thus we can not obtain condition (C1) of Theorem 3.1. Following [24] we construct an auxiliary function which do enjoy the property of function values decreases. Let $\Psi:\mathbb{R}^{n_{1}\times n_{2}}\times\mathbb{R}^{n_{1}\times n_{2}}\rightarrow\left(-\infty,\infty\right]$ be the auxiliary function which is defined as follows

[TABLE]

where $\delta_{1},\delta_{2}>0$, $u_{1}=\left(u_{11},u_{12}\right)\in\mathbb{R}^{n_{1}}\times\mathbb{R}^{n_{2}}$, $u_{2}=\left(u_{21},u_{22}\right)\in\mathbb{R}^{n_{1}}\times\mathbb{R}^{n_{2}}$ and ${\bf u}=\left(u_{1},u_{2}\right)$.

Let $\left\{\left(x_{1}^{k},x_{2}^{k}\right)\right\}_{k\in\mathbb{N}}$ be a sequence generated by iPALM and denote, for all $k\in\mathbb{N}$, $u_{1}^{k}=\left(x_{1}^{k},x_{2}^{k}\right)$, $u_{2}^{k}=\left(x_{1}^{k-1},x_{2}^{k-1}\right)$ and ${\bf u}^{k}=\left(u_{1}^{k},u_{2}^{k}\right)$. We will prove now that the sequence $\left\{{{\bf u}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ and the function $\Psi$ defined above do satisfy conditions (C1), (C2) and (C3) of Theorem 3.1. We begin with proving condition (C1). To this end we will show that there are choices of $\delta_{1}>0$ and $\delta_{2}>0$, such that there exists $\rho_{1}>0$ which satisfies

[TABLE]

It is easy to check that using the notations defined in (4.1), we have, for all $k\in\mathbb{N}$, that

[TABLE]

In order to prove that the sequence $\left\{\Psi\left({\bf u}^{k}\right)\right\}_{k\in\mathbb{N}}$ decreases we will need the following technical result (we provide the proof in Appendix 7).

Lemma 4.1**.**

Consider the functions $g:\mathbb{R}_{+}^{5}\rightarrow\mathbb{R}$ and $h:\mathbb{R}_{+}^{5}\rightarrow\mathbb{R}$ defined as follow

[TABLE]

Let $\varepsilon>0$ and $\bar{\alpha}>0$ be two real numbers for which $0\leq\alpha\leq\bar{\alpha}<0.5\left(1-\varepsilon\right)$. Assume, in addition, that $0\leq L\leq\lambda$ for some $\lambda>0$ and $0\leq\beta\leq\bar{\beta}$ with $\bar{\beta}>0$. If

[TABLE]

then $g\left(\alpha,\beta,\delta_{\ast},\tau_{\ast},L\right)=\varepsilon\delta_{\ast}$ and $h\left(\alpha,\beta,\delta_{\ast},\tau_{\ast},L\right)\geq\varepsilon\delta_{\ast}$.

Based on the mentioned lemma we will set, from now on, the parameters $\tau_{1}^{k}$ and $\tau_{2}^{k}$ for all $k\in\mathbb{N}$, as follow

$\tau_{1}^{k}=\frac{\left(1+\varepsilon\right)\delta_{1}^{k}+\left(1+\beta_{1}^{k}\right)L_{1}(x_{2}^{k})}{1-\alpha_{1}^{k}}\quad\text{and}\quad\tau_{2}^{k}=\frac{\left(1+\varepsilon\right)\delta_{2}^{k}+\left(1+\beta_{2}^{k}\right)L_{2}(x_{1}^{k+1})}{1-\alpha_{2}^{k}}.$

(4.9)

Remark 4.1**.**

If we additionally know that $f_{i}$, $i=1,2$, is convex, then a tighter bound can be used in Proposition 3.2 as described in Remark 3.8. Using the tight bound will improve the possible parameter $\tau_{i}^{k}$ that can be used (*cf. *(4.9)). Indeed, in the convex case, (4.7) and (4.8) are given by

[TABLE]

Thus, the parameters $\tau_{i}^{k}$, $i=1,2$, can be taken in the convex case as follows

[TABLE]

This means that in the convex case, we can take smaller $\tau_{i}^{k}$, $i=1,2$, which means larger step-size in the algorithm. On top of that, in the case that $f_{i}$, $i=1,2$, is convex, it should be noted that a careful analysis shows that in this case the parameters $\alpha_{i}^{k}$, $i=1,2$, can be in the interval $\left[0,1\right)$ and not $\left[0,0.5\right)$ as stated in Lemma 4.1 (see also Assumption B below).

In order to prove condition C1 and according to Lemma 4.1, we will need to restrict the possible values of the parameters $\alpha_{i}^{k}$ and $\beta_{i}^{k}$, $i=1,2$, for all $k\in\mathbb{N}$. The following assumption is essential for our analysis.

Assumption B**.**

Let $\varepsilon>0$ be an arbitrary small number. For all $k\in\mathbb{N}$ and $i=1,2$, there exist $0<{\bar{\alpha}_{i}}<\left(1/2\right)\left(1-\varepsilon\right)$ such that $0\leq\alpha_{i}^{k}\leq{\bar{\alpha}_{i}}$. In addition, $0\leq\beta_{i}^{k}\leq{\bar{\beta}_{i}}$ for some $\bar{\beta_{i}}>0$.

Remark 4.2**.**

It should be noted that using Assumption B, we obtain that $\tau_{1}^{k}\leq\tau_{1}^{+}$ where

[TABLE]

where $\delta_{1}$ is given in (4.7). Similar arguments show that

[TABLE]

Now we will prove a descent property of $\left\{\Psi\left({\bf u}^{k}\right)\right\}_{k\in\mathbb{N}}$.

Proposition 4.3**.**

Let $\left\{{{\bf x}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ be a sequence generated by iPALM which is assumed to be bounded. Suppose that Assumptions A and B hold true. Then, for all $k\in\mathbb{N}$ and $\varepsilon>0$, we have

[TABLE]

where ${\bf u}^{k}=\left({\bf x}^{k},{\bf x}^{k-1}\right)$, $k\in\mathbb{N}$ and $\rho_{1}=\left(\varepsilon/2\right)\min\left\{\delta_{1},\delta_{2}\right\}$ with

[TABLE]

Proof.

From the definition of $\Psi$ (see (4.6)) and Corollary 4.1 we obtain that

[TABLE]

where

[TABLE]

Let $\varepsilon>0$ be an arbitrary. Using (4.9) and (4.12) with the notations of Lemma 4.1 we immediately see that $a_{1}^{k}=g_{1}\left(\alpha_{1}^{k},\beta_{1}^{k},\delta_{1},\tau_{1}^{k},L_{1}(x_{2}^{k})\right)$ and $a_{2}^{k}=g_{2}\left(\alpha_{2}^{k},\beta_{2}^{k},\delta_{2},\tau_{2}^{k},L_{2}(x_{1}^{k+1})\right)$. From Assumptions A and B we get that the requirements of Lemma 4.1 are fulfilled, which means that Lemma 4.1 can be applied. Thus $a_{1}^{k}=\varepsilon\delta_{1}$ and $a_{2}^{k}=\varepsilon\delta_{2}$. Using again the notions of Lemma 4.1, we have that $b_{1}^{k}=h_{1}\left(\alpha_{1}^{k},\beta_{1}^{k},\delta_{1},\tau_{1}^{k},L_{1}(x_{2}^{k})\right)$ and $b_{2}^{k}=h_{2}\left(\alpha_{2}^{k},\beta_{2}^{k},\delta_{2},\tau_{2}^{k},L_{2}(x_{1}^{k+1})\right)$. Thus we obtain from Lemma 4.1 that $b_{1}^{k}\geq\varepsilon\delta_{1}$ and $b_{2}^{k}\geq\varepsilon\delta_{2}$. Hence, for $\rho_{1}=\left(\varepsilon/2\right)\min\left\{\delta_{1},\delta_{2}\right\}$, we have

[TABLE]

where the last equality follows from (4.1). This completes the proof. ∎

Now, we will prove that condition (C2) of Theorem 3.1 holds true for the sequence $\left\{{{\bf u}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ and the function $\Psi$.

Proposition 4.4**.**

Let $\left\{{{\bf x}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ be a sequence generated by iPALM which is assumed to be bounded. Suppose that Assumptions A and B hold true. Assume that ${\bf u}^{k}=\left({\bf x}^{k},{\bf x}^{k-1}\right)$, $k\in\mathbb{N}$. Then, there exists a positive scalar $\rho_{2}$ such that for some ${\bf w}^{k}\in\partial\Psi\left({\bf u}^{k}\right)$ we have

[TABLE]

Proof.

Let $k\geq 2$. By the definition of $\Psi$ (see (4.6)) we have that

[TABLE]

By the definition of $F$ (see (2.1)) and [10, Proposition 1, Page 465] we get that

[TABLE]

From the definition of the proximal mapping (see (2.2)) and the iterative step (2.7) we have

[TABLE]

Writing down the optimality condition yields

[TABLE]

where $\xi_{1}^{k}\in\partial f_{1}\left(x_{1}^{k}\right)$. Hence

[TABLE]

where the last equality follows from step (2.5). By defining

[TABLE]

we obtain from (4.13) that $v_{1}^{k}\in\partial_{x_{1}}F\left({\bf x}^{k}\right)$. Similarly, from the iterative step (2.10), by defining

[TABLE]

we have that $v_{2}^{k}\in\partial_{x_{2}}F\left({\bf x}^{k}\right)$.

Thus, for

[TABLE]

we obtain that

[TABLE]

This means that we have to bound from above the norms of $v_{1}^{k}$ and $v_{2}^{k}$. Since $\nabla H$ is Lipschitz continuous on bounded subsets of $\mathbb{R}^{n_{1}}\times\mathbb{R}^{n_{2}}$ (see Assumption A(v)) and since we assumed that $\left\{{{\bf x}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ is bounded, there exists $M>0$ such that

[TABLE]

where the third inequality follows from (2.6), the fact the sequence $\left\{\tau_{1}^{k}\right\}_{k\in\mathbb{N}}$ is bounded from above by $\tau_{1}^{+}$ (see Remark 4.2) and $\alpha_{1}^{k},\beta_{1}^{k}\leq 1$ for all $k\in\mathbb{N}$. On the other hand, from the Lipschitz continuity of $\nabla_{x_{2}}H\left(x_{1},\cdot\right)$ (see Assumption A(iii)), we have that

[TABLE]

where the third and fourth inequalities follow from (2.9), the fact the sequence $\left\{\tau_{2}^{k}\right\}_{k\in\mathbb{N}}$ is bounded from above by $\tau_{2}^{+}$ (see Remark 4.2) and $\beta_{2}^{k}\leq 1$ for all $k\in\mathbb{N}$. Summing up these estimations, we get from (4.16) that

[TABLE]

where the second inequality follows from the fact that $\left\|{x_{i}^{k}-x_{i}^{k-1}}\right\|\leq\left\|{{\bf x}^{k}-{\bf x}^{k-1}}\right\|$ for $i=1,2$, the third inequality follows from the fact that $\left\|{{\bf x}^{k}-{\bf x}^{k-1}}\right\|\leq\left\|{{\bf u}^{k}-{\bf u}^{k-1}}\right\|$ and the last inequality follows from the fact that $\left\|{{\bf x}^{k}-{\bf x}^{k-1}}\right\|+\left\|{{\bf x}^{k-1}-{\bf x}^{k-2}}\right\|\leq\sqrt{2}\left\|{{\bf u}^{k}-{\bf u}^{k-1}}\right\|$. This completes the proof with $\rho_{2}=\sqrt{2}\left(\tau_{1}^{+}+M+\tau_{2}^{+}+\lambda_{2}^{+}\right)+2\left(\delta_{1}+\delta_{2}\right)$. ∎

So far we have proved that the sequence $\left\{{{\bf u}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ and the function $\Psi$ (see (4.6)) satisfy conditions (C1) and (C2) of Theorem 3.1. Now, in order to get that $\left\{{{\bf u}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ converges to a critical point of $\Psi$, it remains to prove that condition (C3) holds true.

Proposition 4.5**.**

Let $\left\{{{\bf x}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ be a sequence generated by iPALM which is assumed to be bounded. Suppose that Assumptions A and B hold true. Assume that ${\bf u}^{k}=\left({\bf x}^{k},{\bf x}^{k-1}\right)$, $k\in\mathbb{N}$. Then, each limit point in the set $\omega\left({\bf u}^{0}\right)$ is a critical point of $\Psi$.

Proof.

Since $\left\{{{\bf u}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ is assumed to be bounded, the set $\omega\left({\bf u}^{0}\right)$ is nonempty. Thus there exists ${\bf u}^{\ast}=\left(x_{1}^{\ast},x_{2}^{\ast},{\hat{x}_{1}},{\hat{x}_{2}}\right)$ which is a limit point of $\left\{{\bf u}^{k_{l}}\right\}_{l\in\mathbb{N}}$, which is a subsequence of $\left\{{{\bf u}}^{{k}}\right\}_{{k}\in\mathbb{N}}$. We will prove that ${\bf u}^{\ast}$ is a critical point of $\Psi$ (see (4.6)). From condition (C2), for some ${\bf w}^{k}\in\partial\Psi\left({\bf u}^{k}\right)$, we have that

[TABLE]

From Proposition 4.3, it follow that for any $N\in\mathbb{N}$, we have

[TABLE]

Since $F$ is bounded from below (see Assumption A(ii)) and the fact that $\Psi\left(\cdot\right)\geq F\left(\cdot\right)$ we obtain that $\Psi$ is also bounded from below. Thus, letting $N\rightarrow\infty$ in (4.17) yields that

[TABLE]

which means that

[TABLE]

This fact together with condition (C1) implies that $\left\|{{\bf w}^{k}}\right\|\rightarrow 0$ as $k\rightarrow\infty$. Thus, in order to use the closedness property of $\partial\Psi$ (see Remark 3.1) we remain to show that $\left\{\Psi\left({\bf u}^{k}\right)\right\}_{k\in\mathbb{N}}$ converges to $\Psi\left({\bf u}^{\ast}\right)$. Since $f_{1}$ and $f_{2}$ are lower semicontinuous (see Assumption A(i)), we obtain that

[TABLE]

From the iterative step (2.7), we have, for all integer $k$, that

[TABLE]

Thus letting $x_{1}=x_{1}^{\ast}$ in the above, we get

[TABLE]

Choosing $k=k_{l}-1$ and letting $k$ goes to infinity, we obtain

[TABLE]

where we have used the facts that both sequences $\left\{{{\bf x}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ (and therefore $\left\{{{\bf z}}^{{k}}\right\}_{{k}\in\mathbb{N}}$) and $\left\{\tau_{1}^{k}\right\}_{k\in\mathbb{N}}$ (see Remark 4.2) are bounded, $\nabla H$ continuous and that the distance between two successive iterates tends to zero (see (4.19)). For that very reason we also have that $x_{1}^{k_{l}}\rightarrow x_{1}^{\ast}$ as $l\rightarrow\infty$, hence (4.21) reduces to $\limsup_{l\rightarrow\infty}f_{1}\left(x_{1}^{k_{l}}\right)\leq f_{1}\left(x_{1}^{\ast}\right)$. Thus, in view of (4.20), $f_{1}\left(x_{1}^{k_{l}}\right)$ tends to $f_{1}\left(x_{1}^{\ast}\right)$ as $k\rightarrow\infty$. Arguing similarly with $f_{2}$ and $x_{2}^{k+1}$ we thus finally obtain from (4.19) that

[TABLE]

Now, the closedness property of $\partial\Psi$ (see Remark 3.1) implies that ${\bf 0}\in\partial\Psi\left({\bf u}^{\ast}\right)$, which proves that ${\bf u}^{\ast}$ is a critical point of $\Psi$. This proves condition (C3). ∎

Now using the convergence proof methodology of [10] which is summarized in Theorem 3.1 we can obtain the following result.

Corollary 4.2**.**

Let $\left\{{{\bf x}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ be a sequence generated by iPALM which is assumed to be bounded. Suppose that Assumptions A and B hold true. Assume that ${\bf u}^{k}=\left({\bf x}^{k},{\bf x}^{k-1}\right)$, $k\in\mathbb{N}$. If $F$ is a KL function, then the sequence $\left\{{{\bf u}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ converges to a critical point ${\bf u}^{\ast}$ of $\Psi$.

Proof.

The proof follows immediately from Theorem 3.1 since Proposition 4.3 proves that condition (C1) holds true, Proposition 4.4 proves that condition (C2) holds true, and condition (C3) was proved in Proposition 4.5. It is also clear that if $F$ is a KL function then obviously $\Psi$ is a KL function since we just add two quadratic functions. ∎

To conclude the convergence theory of iPALM we have to show that the sequence $\left\{{{\bf x}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ which is generated by iPALM converges to a critical point of $F$ (see (2.1)).

Theorem 4.1** (Convergence of iPALM).**

Let $\left\{{{\bf x}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ be a sequence generated by iPALM which is assumed to be bounded. Suppose that Assumptions A and B hold true. If $F$ is a KL function, then the sequence $\left\{{{\bf x}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ converges to a critical point ${\bf x}^{\ast}$ of $F$.

Proof.

From Corollary 4.2 we have that the sequence $\left\{{{\bf u}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ converges to a critical point ${\bf u}^{\ast}=\left(u_{11}^{\ast},u_{12}^{\ast},u_{21}^{\ast},u_{22}^{\ast}\right)$ of $\Psi$. Therefore, obviously also the sequence $\left\{{{\bf x}}^{{k}}\right\}_{{k}\in\mathbb{N}}$ converges. Let ${\bf x}^{\ast}$ be the limit point of $\left\{{{\bf x}}^{{k}}\right\}_{{k}\in\mathbb{N}}$. Thence $u_{1}^{\ast}=\left(u_{11}^{\ast},u_{12}^{\ast}\right)={\bf x}^{\ast}$ and $u_{2}^{\ast}=\left(u_{21}^{\ast},u_{22}^{\ast}\right)={\bf x}^{\ast}$ (see the discussion on Page 14). We will prove that ${\bf x}^{\ast}$ is a critical point of $F$ (see (2.1)), that is, we have to show that $0\in\partial F\left({\bf x}^{\ast}\right)$. Since ${\bf u}^{\ast}$ is a critical point of $\Psi$, it means that $0\in\partial\Psi\left({\bf u}^{\ast}\right)$. Thus

[TABLE]

which means that

[TABLE]

This proves that ${\bf x}^{\ast}$ is a critical point of $F$. ∎

5 Numerical Results

In this section we consider several important applications in image processing and machine learning to illustrate the numerical performance of the proposed iPALM method. All algorithms have been implemented in Matlab R2013a and executed on a server with Xeon(R) E5-2680 v2 @ 2.80GHz CPUs and running Linux.

5.1 Non-Negative Matrix Factorization

In our first example we consider the problem of using the Non-negative Matrix Factorization (NMF) to decompose a set of facial images into a number of sparse basis faces, such that each face of the database can be approximated by using a small number of those parts. We use the ORL database [32] that consists of $400$ normalized facial images. In order to enforce sparsity in the basis faces, we additionally consider a $\ell_{0}$ sparsity constraint (see, for example, [27]). The sparse NMF problem to be solved is given by

[TABLE]

where $A\in\mathbb{R}^{m\times n}$ is the data matrix, organized in a way that each column of the matrix $A$ corresponds to one face of size $m=64\times 64$ pixels. In total, the matrix holds $n=400$ faces (see Figure 1 for a visualization of the data matrix $A$). The matrix $B\in\mathbb{R}^{m\times r}$ holds the $r$ basis vectors $b_{i}\in\mathbb{R}^{m\times 1}$, where $r$ corresponds to the number of sparse basis faces. The sparsity constraint applied to each basis face requires that the number of non-zero elements in each column vector $b_{i}$, $i=1,2,\ldots,r$ should have less or equal to $s$ non-zero elements. Finally, the matrix $C\in\mathbb{R}^{r\times n}$ corresponds to the coefficient vectors.

The application of the proposed iPALM algorithm to this problem is straight-forward. The first block of variables corresponds to the matrix $B$ and the second block of variables corresponds to the matrix $C$. Hence, the smooth coupling function of both blocks is given by

[TABLE]

The block-gradients and respective block-Lipschitz constants are easily computed via

[TABLE]

The nonsmooth function for the first block, $f_{1}\left(B\right)$, is given by the non-negativity constraint $B\geq 0$ and the $\ell_{0}$ sparsity constraint applied to each column of the matrix $B$, that is,

[TABLE]

Although this function is an indicator function of a nonconvex set, it is shown in [10], that its proximal mapping can be computed very efficiently (in fact in linear time) via

[TABLE]

where $\hat{b_{i}^{+}}=\max\{\hat{b_{i}},0\}$ denotes an elementwise truncation at zero and the operator $T_{s}\left(\hat{b_{i}^{+}}\right)$ corresponds to first sorting the values of $\hat{b_{i}^{+}}$, keeping the $s$ largest values and setting the remaining $m-s$ values to zero.

The nonsmooth function corresponding to the second block is simply the indicator function of the non-negativity constraint of $C$, that is,

[TABLE]

and its proximal mapping is trivially given by

[TABLE]

which is again an elementwise truncation at zero.

In our numerical example we set $r=25$, that is we seek for $25$ sparse basis images. Figure 2 shows the results of the basis faces when running the iPALM algorithm for different sparsity settings. One can see that for smaller values of $s$, the algorithm leads to more compact representations. This might improve the generalization capabilities of the representation.

In order to investigate the properties of iPALM on the inertial parameters, we run the algorithm for a specific sparsity setting ($s=33\%$) using different constant settings of $\alpha_{i}$ and $\beta_{i}$ for $i=1,2$. From our convergence theory (see Proposition 4.3) it follows that the parameters $\delta_{i}$, $i=1,2$, have to be chosen as constants. However, in practice we shall use a varying parameter $\delta_{i}^{k}$, $i=1,2$, and assume that the parameters will become constant after a certain number of iterations. Observe, that the nonsmooth function of the first block is nonconvex ($\ell_{0}$ constraint), while the nonsmooth function of the second block is convex (non-negativity constraint) which will affect the rules to compute the parameters (see Remark 4.1).

We compute the parameters $\tau_{i}^{k}$, $i=1,2$, by invoking (4.7) and (4.8), where we practically choose $\varepsilon=0$. Hence, for the first, completely nonconvex block, we have

[TABLE]

from which it directly follows that $\alpha\in\left[0,0.5\right)$.

For the second block, where the nonsmooth function is convex we follow Remark 4.1 and invoke (4.10) and (4.11) to obtain

[TABLE]

Comparing this parameter to the parameter of the first block we see that now $\alpha\in\left[0,1\right)$ and the value of $\tau$ is smaller by a factor of $2$. Hence, convexity in the nonsmooth function allows for twice larger steps.

In Tables 1 and 2, we report the performance of the iPALM algorithm for different settings of the inertial parameters $\alpha_{i}^{k}$ and $\beta_{i}^{k}$, $i=1,2$ and $k\in\mathbb{N}$. Since we are solving a nonconvex problem, we report the values of the objective function after a certain number of iterations ($100,500,1000$ and $5000$). As already mentioned, setting $\alpha_{i}=\beta_{i}=0$, $i=1,2$, reverts the proposed iPALM algorithm to the PALM algorithm [10]. In order to estimate the (local) Lipschitz constants $L_{1}(x_{2}^{k})$ and $L_{2}(x_{1}^{k+1})$ we computed the exact values of the Lipschitz constants by computing the largest eigenvalues of $C^{k}(C^{k})^{T}$ and $(B^{k})^{T}B^{k}$, respectively. The results are given in table 1. Furthermore, we also implemented a standard backtracking procedure, see for example [6, 24], which makes use of the descent lemma in order to estimate the value of the (local) Lipschitz constant. In terms of iterations of iPALM, the backtracking procedure generally leads to a better overall performance but each iteration also takes more time compared to the exact computation of the Lipschitz constants. The results based on backtracking are shown in table 2.

We tried the following settings for the inertial parameters.

• •

Equal: Here we used the same inertial parameters for both the first and the second block. In case all parameters are set to be zero, we recover PALM. We observe that the use of the inertial parameters can speed up the convergence but for too large inertial parameters we also observe not as good results as in the case that we use PALM, i.e., no inertial is used. Since the problem is highly nonconvex the final value of the objective function can be missleading since it could correspond to a bad stationary point.

• •

Double: Since the nonsmooth function of the second block is convex, we can take twice larger inertial parameters. We observe an additional speedup by taking twice larger inertial parameters. Again, same phenomena occur here, too large inertial parameters yields inferior performances of the algorithm.

• •

Dynamic: We also report the performance of the algorithm in case we choose dynamic inertial parameters similar to accelerated methods in smooth convex optimization [21]. We use $\alpha_{i}^{k}=\beta_{i}^{k}=\left(k-1\right)/\left(k+2\right)$, $i=1,2$, and we set the parameters $\tau_{1}^{k}=L_{1}(x_{2}^{k})$ and $\tau_{2}^{k}=L_{2}(x_{1}^{k+1})$. One can see that this setting outperforms the other settings by a large margin. Although our current convergence analysis does not support this setting, it shows the great potential of using inertial algorithms when tackling nonconvex optimization problems. The investigation of the convergence properties in this setting will be subject to future research.

Finally, we also compare the proposed algorithm to the iPiano algorithm [24], which is similar to the proposed iPALM algorithm but does not make use of the block structure. Note that in theory, iPiano is not applicable since the gradient of the overall problem is not Lipschitz continuous. However, in practice it turns out that using a backtracking procedure to determine the Lipschitz constant of the gradient is working and hence we show comparisons. In the iPiano algorithm, the step-size parameter $\tau$ ($\alpha$ in terms of [24]) was set as

[TABLE]

from which it follows that the inertial parameter $\beta$ can be chosen in the interval $\left[0,0.5\right)$. We used $\beta=0.4$ since it give the best results in our tests. Observe that the performance of iPALM is much better compared to the performance of iPiano. In terms of CPU time, one iteration of iPiano is clearly slower than one iteration of iPALM using the exact computation of the Lipschitz constant, because the backtracking procedure is computationally more demanding than the computation of the largest eigenvalues of $C^{k}(C^{k})^{T}$ and $(B^{k})^{T}B^{k}$. Comparing the iPiano algorithm to the version of iPALM which uses backtracking, one iteration of iPiano is faster since iPALM needs to backtrack the Lipschitz constants for both of the two blocks.

5.2 Blind Image Deconvolution

In our second example, we consider the well-studied (yet challenging) problem of blind image deconvolution (BID). Given a blurry and possibly noisy image $f\in\mathbb{R}^{M}$ of $M=m_{1}\times m_{2}$ pixels, the task is to recover both a sharp image $u$ of the same size and the unknown point spread function $b\in\mathbb{R}^{N}$, which is a small $2D$ blur kernel of size $N=n_{1}\times n_{2}$ pixels. We shall assume that the blur kernel is normalized, that is $b\in\Delta^{N}$, where $\Delta^{N}$ denotes the standard unit simplex defined by

[TABLE]

Furthermore, we shall assume that the pixel intensities of the unknown sharp image $u$ are normalized to the interval $\left[0,1\right]$, that is $u\in U^{M}$, where

[TABLE]

We consider here a classical blind image deconvolution model (see, for example, [28]) defined by

[TABLE]

The first term is a regularization term which favors sharp images and the second term is a data fitting term that ensures that the recovered solution approximates the given blurry image. The parameter $\lambda>0$ is used to balance between regularization and data fitting. The linear operators $\nabla_{p}$ are finite differences approximation to directional image gradients, which in implicit notation are given by

[TABLE]

for $1\leq i\leq m_{1}$ and $1\leq j\leq m_{2}$. We assume natural boundary conditions, that is $\left(\nabla_{p}\right)_{i,j}=0$, whenever the operator references a pixel location that lies outside the domain. The operation $u\ast_{m_{1},m_{2}}b$ denotes the usual 2D modulo - $m_{1},m_{2}$ discrete circular convolution operation defined by (and interpreting the image $u$ and the blur kernel $b$ as 2D arrays)

[TABLE]

For ease the notation, we will rewrite the 2D discrete convolutions as the matrix vector products of the form

[TABLE]

where $K\left(b\right)\in\mathbb{R}^{M\times M}$ is a sparse matrix, where each row holds the values of the blur kernel $b$, and $K\left(u\right)\in\mathbb{R}^{M\times N}$ is a dense matrix, where each column is given by a circularly shifted version of the image $u$. Finally, the function $\phi\left(\cdot\right)$ is a differentiable robust error function, that promotes sparsity in its argument. For a vector $x\in\mathbb{R}^{M}$, the function $\phi$ is defined as

[TABLE]

where $\theta>0$ is a parameter. Here, since the argument of the function are image gradients, the function promotes sparsity in the edges of the image. Hence, we can expect that sharp images result in smaller values of the objective function than blurry images. For images that can be well described by piecewise constant functions (see, for example, the books image in Figure 3), such sparsity promoting function might be well suited to favor sharp images, but we would like to stress that this function could be a bad choice for textured images, since a sharp image usually has much stronger edges than the blurry image. This often leads to the problem that the trivial solution ($b$ being the identity kernel and $u=f$) has a lower energy compared to the true solution.

In order to apply the iPALM algorithm, we identify the following functions

[TABLE]

which is smooth with block Lipschitz continuous gradients given by

[TABLE]

where the operation $\mathrm{vec}\left(\cdot\right)$ denotes the formation of a vector from the values passed to its argument. The nonsmooth function of the first block is given by

[TABLE]

and the proximal map with respect to $f_{1}$ is computed as

[TABLE]

The nonsmooth function of the second block is given by the indicator function of the unit simplex constraint, that is,

[TABLE]

In order to compute the proximal map with respect to $f_{2}$, we use the algorithm proposed in [13] which computes the projection onto the unit simplex in $O(N\log N)$ time.

We applied the BID problem to the books image of size $m_{1}\times m_{2}=495\times 323$ pixels (see Figure 3). We set $\lambda=10^{6}$, $\theta=10^{4}$ and generated the blurry image by convolving it with a s-shaped blur kernel of size $n_{1}\times n_{2}=31\times 31$ pixels. Since the nonsmooth functions of both blocks are convex, we can set the parameters as (compare to the first example)

[TABLE]

To determine the values of the (local) Lipschitz constants, we used again a backtracking scheme [6, 24]. In order to avoid the trivial solution (that is, $u=f$ and $b$ being the identity kernel) we took smaller descent steps in the blur kernel which was realized by multiplying $\tau_{2}^{k}$, $k\in\mathbb{N}$, by a factor of $c=5$. Note that this form of “preconditioning” does not violate any step-size restrictions as we can always take larger values for $\tau$ than the value computed in (5.12).

Table 3 shows an evaluation of the iPALM algorithm using different settings of the inertial parameters $\alpha_{i}$ and $\beta_{i}$ for $i=1,2$. First, we observe that the use of inertial parameters lead to higher values of the objective function after a smaller number of iterations. However, for a larger number of iterations the use of inertial forces leads to significantly lower values. Again, the use of dynamic inertial parameters together with the parameter $\tau_{1}^{k}=L_{1}(x_{2}^{k})$ and $\tau_{2}^{k}=L_{2}(x_{1}^{k+1})$ leads to the best overall performance. In Figure 3 we show the results of the blind deconvolution problem. One can see that the quality of the recovered image as well as the recovered blur kernel is much better using inertial forces. Note that the recovered blur kernel is very close to the true blur kernel but the recovered image appears slightly more piecewise constant than the original image.

5.3 Convolutional LASSO

In our third experiment we address the problem of sparse approximation of an image using dictionary learning. Here, we consider the convolutional LASSO model [36], which is an interesting variant of the well-known patch-based LASSO model [25, 1] for sparse approximations. The convolutional model inherently models the transnational invariance of images, which can be considered as an advantage over the usual patch-based model which treats every patch independently.

The idea of the convolutional LASSO model is to learn a set of small convolution filters $d_{j}\in\mathbb{R}^{l\times l}$, $j=1,2,\ldots,p$, such that a given image $f\in\mathbb{R}^{m\times n}$ can be written as $f\approx\sum_{j=1}^{p}d_{j}\ast_{m,n}v_{j}$, where $\ast_{m,n}$ denotes again the $2D$ modulo $m,n$ discrete circular convolution operation and $v_{j}\in\mathbb{R}^{m\times n}$ are the corresponding coefficient images which are assumed to be sparse. In order to make the convolution filters capture the high-frequency information in the image, we fix the first filter $d_{1}$ to be a Gaussian (low pass) filter $g$ with standard deviation $\sigma_{l}$ and we set the corresponding coefficient image $v_{j}$ equal to the initial image $f$. Furthermore, we assume that the remaining filters $d_{j}$, $j=2,3,\ldots,p$ have zero mean as well as a $\ell_{2}$-norm less or equal to $1$. In order to impose a sparsity prior on the coefficient images $v_{j}$ we make use of the $\ell_{1}$-norm. The corresponding objective function is hence given by

[TABLE]

where the parameter $\lambda>0$ is used to control the degree of sparsity. It is easy to see that the convolutional LASSO model nicely fits to the class of problems that can be solved using the proposed iPALM algorithm. We leave the details to the interested reader. In order to compute the (local) Lipschitz constants we again made use of a backtracking procedure and the parameters $\tau_{i}^{k}$, $i=1,2$ were computed using (5.12).

We applied the convolutional LASSO problem to the Barbara image of size $512\times 512$ pixels, which is shown in Figure 4. The Barbara image contains a lot of stripe-like texture and hence we expect that the learned convolution filters will contain these characteristic structures. In our experiment, we learned a dictionary made of $81$ filter kernels, each of size $9\times 9$ pixels. The regularization parameter $\lambda$ was set to $\lambda=0.2$. From Figure 4 it can be seen that the learned convolution filters indeed contain stripe-like structures of different orientations but also other filters that are necessary to represent the other structures in the image. Table 4 summarizes the performance of the iPALM algorithm for different settings of the inertial parameters. From the results, one can see that larger settings of the inertial parameters lead to a consistent improvement of the convergence speed. Again, using a dynamic choice of the inertial parameters clearly outperforms the other settings.

For illustration purposes we finally applied the learned dictionary to denoise a noisy variant of the same Barbara image. The noisy image has been generated by adding zero-mean Gaussian noise with standard deviation $\sigma=0.1$ to the original image. For denoising we use the previously learned dictionary $d$ and minimizing the convolutional LASSO problem only with respect to the coefficient images $v$. Note that this is a convex problem and hence,it can be efficiently minimized using for example the FISTA algorithm [6]. The denoised image is again shown in Figure 4. Observe that the stripe-like texture is very well preserved in the denoised image.

6 Conclusion

In this paper we proposed iPALM which an inertial variant of the Proximal Alternating Linearized Minimization (PALM) method proposed in [10] for solving a broad class of nonconvex and nonsmoooth optimization problems consisting of block-separable nonsmooth, nonconvex functions with easy to compute proximal mappings and a smooth coupling function with block-Lipschitz continuous gradients. We studied the convergence properties of the algorithm and provide bounds on the inertial and step-size parameters that ensure convergence of the algorithm to a critical point of the problem at hand. In particular, we showed that in case the objective function satisfies the Kurdyka-Łojasiewicz (KL) property, we can obtain finite length property of the generated sequence of iterates. In several numerical experiments we show the advantages of the proposed algorithm to minimize a number of well-studied problems and image processing and machine learning. In our experiments we found that choosing the inertial and step-size parameters dynamically, as pioneered by Nesterov [22], leads to a significant performance boost, both in terms of convergence speed and convergence to a “better” critical point of the problem. Our current convergence theory does not support this choice of parameters but developing a more general convergence theory will be interesting and a subject for future research.

7 Appendix A: Proof of Lemma 4.1

We first recall the result that should be proved.

Lemma 7.1**.**

Consider the functions $g:\mathbb{R}_{+}^{5}\rightarrow\mathbb{R}$ and $h:\mathbb{R}_{+}^{5}\rightarrow\mathbb{R}$ defined as follow

[TABLE]

Let $\varepsilon>0$ and $\bar{\alpha}>0$ be two real numbers for which $0\leq\alpha\leq\bar{\alpha}<0.5\left(1-\varepsilon\right)$. Assume, in addition, that $0\leq L\leq\lambda$ for some $\lambda>0$ and $0\leq\beta\leq\bar{\beta}$ with $\bar{\beta}>0$. If

[TABLE]

then $g\left(\alpha,\beta,\delta_{\ast},\tau_{\ast},L\right)=\varepsilon\delta_{\ast}$ and $h\left(\alpha,\beta,\delta_{\ast},\tau_{\ast},L\right)\geq\varepsilon\delta_{\ast}$.

Proof.

From the definition of $g$ we immediately obtain that

[TABLE]

this proves the first desired result. We next simplify $h\left(\alpha,\beta,\delta_{\ast},\tau_{\ast},L\right)-\varepsilon\delta_{\ast}$ as follows

[TABLE]

Thus, we only remain to show that

[TABLE]

Indeed, simple manipulations yields the following equivalent inequality

[TABLE]

Using now (7.1) and the facts that $\alpha\leq{\bar{\alpha}}$, $\beta\leq\bar{\beta}$ and $0<L\leq\lambda$ we obtain that (7.3) holds true. This completes the proof of the lemma. ∎