Inertial nonconvex alternating minimizations for the image deblurring
Tao Sun Roberto Barrio Marcos Rodriguez Hao Jiang
math.OC cs.CV cs.MM eess.IV





Inertial nonconvex alternating minimizations for the image deblurring
    Tao Sun, Roberto Barrio, Marcos Rodríguez, and Hao Jiang 
This work was
supported by the National Key Research and Development Program of China (2017YFB0202003), and National Natural Science Foundation of Hunan Province in China (2018JJ3616), and by the Spanish Research projects MTM2015-64095-P, PGC2018-096026-B-I00, and European Regional Development Fund and Diputación General de Aragón (E24-17R).

Tao Sun and Hao Jiang are with the
College of Computer, National University of Defense Technology,
Changsha, 410073, Hunan,  China (e-mail: , ).

Roberto Barrio is with Departamento de Matemática Aplicada, IUMA and CoDy group, Universidad de Zaragoza, Pedro Cerbuna, 12, Zaragoza, 50009, Spain. (e-mail: ).

Marcos Rodríguez is with CoDy group, Zaragoza, Universidad de Zaragoza, Pedro Cerbuna, 12, Zaragoza, 50009, Spain.
(e-mail: ).



    Received: date / Accepted: date
============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================






In image processing, Total Variation (TV) regularization models are commonly used to recover blurred images. One of the most efficient and popular methods to solve the convex TV problem is the Alternating Direction Method of Multipliers (ADMM) algorithm, recently extended using the inertial proximal point method. Although all the classical studies focus on only a convex formulation, recent articles are paying increasing attention to the nonconvex methodology  due to its good numerical performance and properties. In this paper, we propose to extend the classical formulation with a novel nonconvex Alternating Direction Method of Multipliers with the Inertial technique (IADMM). Under certain assumptions on the parameters, we prove the convergence  of the algorithm with the help of the Kurdyka-Łojasiewicz property.
We also present numerical simulations on classical TV image reconstruction problems  to illustrate the efficiency  of the new algorithm  and its behavior compared with the well established ADMM method.



Inertial algorithms, nonconvex method, Kurdyka-Łojasiewicz property, image deblurring, ADMM, inertial proximal ADMM.








§ INTRODUCTION

Denoising and deblurring have numerous applications in communications, control, machine learning, and many other fields of engineering and
science. Restoration of distorted images is, from the theoretical, as well as from the practical point of view, one of the most interesting and important problems of image processing. One special case is the blurring, due, for instance, to incorrect focus and/or blurring due to movement, or added Gaussian noise (a Gaussian blur).

A mathematical model for the process of  blurred images can be expressed as follows. Let Ω≡ℝ^N^2 be a two-dimensional index set representing the image domain, ũ∈Ω
be the original image, f̃ be the observed image, and K̃ be a linear blurring operator. Then, the blurred image can be written <cit.> as

    f̃=K̃(ũ)+e,

where e is an unknown additive noise vector. In this paper, the  blurring operator K̃ is assumed to be known, otherwise, one will deal with   the blind image deblurring problem <cit.>, in which  K̃   also needs to be solved.


In image processing one typically aims at recovering an image from noisy data while still
keeping edges in the image, and this goal is the main reason of the tremendous success of the Total Variation (TV) regularization <cit.> for solving the deblurring problem (although other methods are also used). The TV method can be presented as

    min_ũ(1/2K̃ũ-f̃^2_L^2(Ω)+σ∇ũ_L^1(Ω^*)),

being σ>0 a parameter and where ∇  is the gradient operator, Ω^*:=∇(Ω), and the norms ṽ_L^1(Ω^*):=∑_i,j∈Ω^*|ṽ_i,j|, ũ_L^2(Ω):=√(∑_i,j∈Ω|ũ_i,j|^2).

In most situations, rather than directly minimizing the
support of the image, one is interested in minimizing the support of the gradient of
the recovered image. In most references, the convex methodology is considered <cit.>, but in recent years, some nonconvex methods have been developed <cit.>. The use of a suitable
 nonconvex and nondifferentiable
function allows
possibly a smaller number of measurements than the convex one in compressed sensing <cit.>. In <cit.> the authors showed that nonconvex regularization terms in
total variation-based image restoration yields even better edge preservation when compared
to the convex-type regularization. Moreover, they showed that it seems to be
also more robust with respect to noise.
Nonconvex regularization in image restoration poses significant challenges on the existence of solutions
of associated minimization problems and on the development of efficient solution algorithms.

The main  difference  between the convex and nonconvex methods is replacing the l_1-norm of the variational term by the nonconvex and nondifferentiable
function  v_φ(Ω^*):=∑_i,j∈Ω^*φ(v_i,j) that uses the nonconvex regulation function φ, and that we refer to as the ℓ^q semi-norm (0<q<1). Therefore, the general nonconvex deblurring model is presented as

    min_ũ(1/2K̃ũ-f̃^2_L^2(Ω)+σ∇ũ_φ(Ω^*)).



Many efficient numerical
algorithms have been developed for solving the TV regularization problem.
One of the most efficient methods for the convex problem (<ref>) is the Alternating Direction Method of Multipliers (ADMM) algorithm  <cit.>. In the general case (convex and nonconvex cases depend of the function φ), the method is constructed by introducing an auxiliary variable ṽ, which actually represents ∇ũ,  to reformulate (<ref>)  into a composite optimization problem with linear constraints.
The augmented
Lagrange dual  function is   then

    ℒ̃_δ^φ(ũ,ṽ,p̃):=1/2K̃ũ-f̃^2_L^2(Ω)+σṽ_φ(Ω^*)
               -⟨p̃,∇ũ-ṽ⟩+δ/2∇ũ-ṽ_L^2(Ω^*)^2,

where δ>0 is a parameter, and the norm X_L^2(Ω^*):=√(∑_i,j∈Ω^*|X_i,j|^2). If φ(·)=|·|_1, we use the notation ℒ̃_δ^1(ũ,ṽ,p̃) for representing ℒ̃_δ^φ(ũ,ṽ,p̃). Now, the standard convex ADMM method (φ(·)=|·|_1) for the deblurring problem can be presented as
 
    {[                        ṽ^k+1                            =   argmin_ṽℒ̃_δ^1(ũ^k,ṽ,p̃^k);                        ũ^k+1                            = argmin_ũℒ̃_δ^1(ũ,ṽ^k+1,p̃^k);                       p̃^k+1                            =         p̃^k-δ(∇ũ^k+1-ṽ^k+1);                              ].


The earlier analyses of convergence and performance of the ADMM algorithms directly depended on the existing results of ADMM framework <cit.>. More recently, motivated by acceleration techniques proposed in <cit.>, inertial algorithms have been proposed for many areas such as (distributed) optimization and imaging sciences in references  <cit.>. The ideas of the inertial strategy have been also applied to ADMM in <cit.>; and under several  assumptions in convex case, some convergence results are proved in those articles. As the nonconvex penalty functions perform more efficiently in some applications, as above commented, nonconvex ADMM has been also developed and studied <cit.>.  The main goal of this paper is to propose a new algorithm that combines the nonconvex methods and the inertial strategy organically.

 In this paper, when φ is nonconvex, we consider a new  inertial scheme for the image deblurring model (<ref>). One of the main differences (and new difficulties) with the convex ADMM, is that in order to properly define the nonconvex ADMM some extra assumptions are needed to prove the convergence. First, at least one of the objective functions has to be smooth. And more,  matrix corresponding to the smooth function is required to be injective, i.e., reversible. Thus, a direct application of the ADMM scheme to the image deblurring model cannot guarantee the convergence because the operator ∇  fails to be injective (although the numerical performance may be  good in some cases). Considering this, we first modify the model (<ref>), and then we develop the new nonconvex inertial ADMM.
By using the  Kurdyka-Łojasiewicz property, we prove the convergence of the new algorithm under several requirements on the parameters. In opposite to the convex case, selecting a suitable parameter δ is crucial to obtain the convergence of the new algorithm. In order to make the method more useful, we
provide a probabilistic strategy for selecting a suitable δ.


The rest of the paper is organized as follows. In Section <ref> we collect some mathematical preliminaries needed for the convergence analysis. Section <ref> presents the details for the new algorithm (inertial alternating minimization algorithm, IADMM) including the schemes and parameters. In Section <ref>, we prove the convergence of the new algorithm.  Section <ref> reports the numerical results and compares the algorithm with convex and nonconvex ADMM. Section <ref> gives some conclusions. Finally, we provide in the Appendixes all the detailed proofs of the proposed results.




§ MATHEMATICAL TOOLS


In this section we present the definitions and basic properties of the subdifferentials
and the Kurdyka-Łojasiewicz functions used later in the convergence analysis.
The basic notations  used in this paper are detailed in Table <ref>.



 §.§ Subdifferentials

We collect several definitions as well as some useful properties in variational and convex
analysis (see the monographes <cit.>).
 For any matrix A, we define A^* to be the adjoint of A.

 Let  J: ℝ^N → (-∞, +∞] be a proper and lower semicontinuous function. The (limiting) subdifferential, or simply the subdifferential, of J at x∈ℝ^N, written as ∂ J(x), is defined  as

    ∂ J(x):={u∈ℝ^N: ∃   x^k→ x, u^k→ u, such  that
    lim_y≠ x^kinf_y→ x^kJ(y)-J(x^k)-⟨ u^k, y-x^k⟩/y-x^k_2≥ 0 }.



It is easy to verify that the Fréchet subdifferential is convex and closed while the subdifferential is closed. When J is convex,  the definition agrees with the subgradient in convex analysis <cit.>.
We call J is strongly convex with constant a ∈ℝ^+ if for any x,y∈dom(J) and any v∈∂ J(x), it holds that
J(y)≥ J(x)+⟨ v,y-x⟩+a/2y-x_2^2.
And J is called as  L-gradient continuous (Lipschitz) with constant L>0 if
∇ J(x)-∇ J(y)_2≤ Lx-y_2.
Noting the closedness of the subdifferential, we have the following simple proposition.




If v^k∈∂ J(x^k), lim_kv^k=v and lim_kx^k=x, then we have
v∈∂ J(x).



A necessary condition for x∈ℝ^N to be a minimizer of J(x) is

    0∈∂ J(x),

which is also sufficient
when J is convex.
A point that satisfies (<ref>) is called (limiting) critical point. The set of critical points of J(x) is denoted by crit(J).



With these basics, we can easily obtain the following proposition.


If  (u̅,v̅,p̅) is a critical point of ℒ^φ_δ whose definition given in (<ref>), it must hold that

    {[         -p̅           ∈     σ∂v̅_φ,;       ∇^*p̅           = K^*(Ku̅-f),;         ∇u̅           =         v̅. ].





Finally, the proximal map of J is defined as

    𝒮_J(x)∈argmin_y{J(y)+1/2x-y_2^2}.

Note that J can be nonconvex. If J is convex, 𝒮_J is a point-to-point operator; otherwise, it may be  point-to-set.



 §.§ Kurdyka-Łojasiewicz property


In this paper the convergence analysis is based on the Kurdyka-Łojasiewicz functions, originated in the seminal works of Łojasiewicz <cit.> and Kurdyka <cit.>.  This kind of functions has played a key role in several recent convergence results on nonconvex minimization problems and they are ubiquitous in applications.



(a) The function J: ℝ^N → (-∞, +∞] is said to have the Kurdyka-Łojasiewicz property at x∈dom(∂ J) if there
 exist η∈ (0, +∞], a neighborhood U of x and a continuous concave function ρ: [0, η)→ℝ^+ such that

  
  * ρ(0)=0.
  
  * ρ is 𝒞^1 on (0, η).
  
  * For all s∈(0, η), ρ'(s)>0.
  
  * For all x in U⋂{x|J(x)<J(x)<J(x)+η}, the Kurdyka-Łojasiewicz inequality holds

    ρ'(J(x)-J(x))·dist(0,∂ J(x))≥ 1.



(b) Proper lower semicontinuous functions which satisfy the Kurdyka-Łojasiewicz inequality at each point of dom(∂ J) are called KŁ functions.



 There are large sets of functions that are KŁ functions <cit.>.



Let J: ℝ^N→ℝ be a proper lower semi-continuous function and Π be a compact set. If J is a constant on Π and J satisfies the KŁ property at each point on Π, then there exists a concave function ρ satisfying the four properties given in Definition <ref>, and constants η,ε>0 such that for any x∈Π and any x satisfying that dist(x,Π)<ε and f(x)<f(x)<f(x)+η, it holds that

    ρ'(J(x)-J(x))·dist(0,∂ J(x))≥ 1.





§ NONCONVEX IADMM ALGORITHM


In this section we introduce the new extended Inertial Alternating Direction Method of Multipliers (ADMM) algorithm for nonconvex functions.



In this  paper, we consider Ω=[1,2,…,N]× [1,2,…,N] (equivalent to space ℝ^N^2) as the two-dimensional index set representing the image domain.  In this case, the image variable constrained on Ω is actually a N× N matrix.  We use the symbol ũ to present its vectorization (a vector of all the columns of the image variable).
And then the original total variation operator then enjoys the following  form

    ∇ (ũ)=(
                     [ I_N⊗ D; D⊗ I_N;        ])·ũ,

where  I_N the identity matrix and D the banded matrix 
    D=(
                                                              [  1 -1      ;     1  ⋱      ;        ⋱ -1   ;           1 -1;    ])_(N-1)× N.

If we directly apply  the inertial ADMM, the convergence is hard to be proved as ∇ fails to be injective. Therefore, we need to modify the image deblurring model (<ref>). To that goal, we define

    𝒯:=(
                         [ I_N⊗ D       ;        D⊗ I_N;        ]),    u:=(
                                                                    [  u_1∈ℝ^N^2; u_2 ∈ℝ^N^2;            ]).

Obviously, we have 𝒯∈ℝ^2N(N-1)× 2N^2.
Noting that

    [             rank(𝒯)=rank(I_N⊗ D)+rank(D⊗ I_N); =rank(I_N)·rank(D)+rank(D)·rank(I_N)=2N(N-1), ]

and thus, 𝒯 is injective.
The following technical lemma focuses on giving a lower bound for the operator 𝒯^*.



For any x∈ℝ^2N(N-1), it holds that

    𝒯^*x≥x/θ,

where θ=1/(2sin(π/2N)).



Then, the image deblurring model (<ref>) is equivalent to
 
    min_u_1=u_2(1/2K̃u_1-f̃^2+σ𝒯 u_φ).

 Instead, we consider its extended penalty form
  
    min_u(1/2K̃u_1-f̃^2+σ𝒯 u_φ+β^2/2u_1-u_2^2),

 where β>0 is a large weight parameter. Therefore, we apply the nonconvex  inertial ADMM to
   
    min_u∈ℝ^2N^2(1/2Ku-f^2+σ𝒯 u_φ),

 where K=(
            [  K̃   0;  β𝕀 -β𝕀;     ]) and f=(
            [ f̃;  0;    ]).
 This leads us to define the function

    ℒ_δ^φ(u,v,p):=1/2Ku-f^2
           +σv_φ-⟨ p,𝒯 u-v⟩+δ/2𝒯 u-v^2.


Inertial methods have witnessed  great success in convex ADMM and nonconvex first-order algorithms. In the nonconvex optimization community, the inertial style ADMM has never been proposed and analyzed.  The convex inertial  ADMM has been proposed in <cit.>, in which one  first uses the “inertial method" to refresh the current sequence with last iteration, and then performs the ADMM scheme with the updated variables. However, the direct extension of convex ADMM   is not allowed   in the nonconvex settings.  This is because without convexity, several descents are heavily dependent on the continuities of the functions, which φ may fail to obey.  And the difference of function values at two  different   points is hard to estimate, which leads to  troubles in the convergence proof. Thus, in the updating of v^k+1, we used u^k rather than the updated one. And the nonconvex IADMM scheme proposed in this paper is defined as follows

    {[                      (u^k, v^k, p^k)                                    =                      (u^k, v^k, p^k);                                      +α(u^k-u^k-1, v^k-v^k-1, p^k-p^k-1),;                                v^k+1                                    =            argmin_vℒ_δ^φ(u^k,v,p^k),;                                u^k+1                                    =          argmin_uℒ_δ^φ(u,v^k+1,p^k),;                                p^k+1                                    =                p^k-δ(𝒯 u^k+1-v^k+1), ].

where α>0 is a free parameter chosen by the user.   Actually, if α=0, the algorithm then will reduce to basic ADMM.






Now we can focus on rewriting the inertial scheme for the image deblurring model (<ref>).
First, we rearrange the minimization for v^k+1,

    v^k+1=𝒮_σ/δ·_φ(Ω^*)(𝒯 u^k-p^k/δ),

being S_J the proximal map of J (<ref>).
For a matrix v, and indices (i,j)∈Ω^*,

    [𝒮_σ/δ·_φ(v)]_i,j=
    𝒮_σ/δφ(v_i,j).

The scheme for updating u^k+1 can be rewritten as

    K^*(Ku^k+1-f)-𝒯^* p^k+δ𝒯^*(𝒯 u^k+1-v^k+1)=0.

That is also

    u^k+1=(K^*K+δ𝒯^*𝒯)^-1(K^* f+δ𝒯^*v^k+1+𝒯^* p^k).


Taking into account (<ref>), (<ref>) and (<ref>) we propose a nonconvex inertial version of the IADMM algorithm (Algorithm <ref>).




Assumption 1:  We assume that Null(K_0)⋂Null(𝒯)={0}, where K_0:=(
            [  K̃   0;   𝕀 - 𝕀;     ]). And the minimum single value is given as ν. 

          

This hypothesis also indicates that the matrix  (
     [   𝒯; K_0;     ]) is reversible.
Note that the rank of 𝒯 is 2N(N-1).   Then, the assumed hypothesis is easy to be satisfied.


We remark that u^k+1 is the minimizer of ℒ_δ^φ(u,v^k+1,p^k), which is strongly convex with λ_min(K^*K+δ𝒯^*𝒯). If we set δ,ρ>1, then we have

    ℒ_δ^φ(u^k,v^k+1,p^k)- ℒ_δ^φ(u^k+1,v^k+1,p^k)
       ≥λ_min(K^*K+δ𝒯^*𝒯)/2u^k+1-u^k^2≥ν/2u^k+1-u^k^2,

where we used the fact λ_min(K^*K+δ𝒯^*𝒯)≥λ_min(K^*_0K_0+𝒯^*𝒯) when δ,β≥ 1.

The following problem is what ν exactly is. In a real situation, the dimensions of 𝒯 are large, and so, a direct calculation leads to a large computational cost. Therefore, we provide a probabilistic method to estimate a suitable value of θ. If K^*K+𝒯^*𝒯 is reversible, it is easy to see that
1/ν=(K^*K+𝒯^*𝒯)^-1_2. Then, if we obtain a bound  (K^*K+𝒯^*𝒯)^-1_2≤ c, we then have ν≥1/c. To that goal we employ a Lemma proposed in <cit.>:

For a fixed positive integer M and a
real number b>1, and given an independent family {w_i}_i=1,2,…,M of standard
Gaussian vectors, we have that

    (K^*K+𝒯^*𝒯)^-1_2
    ≤ b√(2/π)max_i=1,2,…,M(K^*K+𝒯^*𝒯)^-1w_i_2

with probability at least 1-b^-M.



Note that for computing (K^*K+𝒯^*𝒯)^-1w_i we just need several FFT and inverse FFT. Therefore, its  computational cost is low (𝒪(N log N), and the estimation of ν is very fast.







§ CONVERGENCE ANALYSIS


This section consists of two parts and provides a complete analysis of the convergence problem of the nonconvex IADMM algorithm.    The first subsection contains the main convergence results, the proof sketch, the difficulties in the proof and theoretical  contributions. While the second subsection introduces the necessary technical lemmas.
Assumption 1 holds through this section.



 §.§ Main results


Assume that the free parameter δ satisfies the condition

    δ>max{1,  6θ^2K_2^4/ν+7α^2θ^2K_2^4/ν}

with θ=1/(2sin(π/2N)), then
any cluster point (u^*,v^*,p^*) is also a critical point of ℒ^φ_δ.



Theorem <ref> describes the stationary point convergence result for the IADMM method, which is free of using the KŁ property of the functions.  If the KŁ property is further assumed, the sequence convergence can be proved giving us the Theorem <ref>.


 Let   condition (<ref>) hold, and the auxiliary function F (<ref>) be a KŁ function. Then the sequence {(u^k,v^k,p^k)}_k=0,1,2,…, generated by Algorithm <ref>, converges to a critical point of ℒ^φ_δ.




The proof  can be divided into two parts, and in order to help the reader we first give a brief sketch of the proof:

I. In the first part we introduce  an auxiliary sequence {w^k}_k=0,1,2,…, where  w^k  are composite points from {(u^k,v^k,p^k)}_k=0,1,2,….  An auxiliary function F is also proposed. In Lemma <ref>, we prove a “sufficient descent condition” of the values of F at
{w^k}_k=0,1,2,…, i.e.,

    F(w^k)-F(w^k+1)≥hu^k+1-u^k^2,

where h>0 is a positive constant.

II. We prove a “relative error condition" of
{w^k}_k=0,1,2,…, i.e., there exists v^k+1∈∂ F(w^k+1) such that

    v^k+1≤γ(u^k+1-u^k+u^k-u^k-1),

where γ>0 is a positive constant.  Note that this condition is different from the “real” relative error condition proposed in <cit.>.

The major  difficulty in  deriving these two conditions is the use of inertial terms, with which  the descent values are lower bounded  by  v^k+1-v^k^2 and p^k+1-p^k^2 rather than  v^k+1-v^k^2 and p^k+1-p^k^2. Similarly, the  relative error is bounded by v^k+1-v^k and p^k+1-p^k. The relative error can be expanded by triangle inequalities, which is relatively proved. However, for the sufficient descent, the use of the triangle inequalities is much more difficult and technical for the lower boundedness.

The theoretical contributions in this paper are two-fold. The first one is, of course, dealing with the difference caused by the inertial terms. This part also includes how to design the scheme of the algorithm, whose details have been presented in previous section.   The second one is to determine the parameters for IADMM  applied to the image deblurring.

The main results  can be proved with the following  lemmas.




 §.§ Technical lemmas



Let the sequence {(u^k,v^k,p^k)}_k=0,1,2,… be generated by Algorithm <ref> to solve problem (<ref>), then

    p^k-p^k+1≤θK_2^2  u^k+1-u^k,

where θ is given in Lemma <ref> and K_2 denotes the spectral radius of K.



Now we provide the main technical lemma that states the descent condition for a suitable function of the sequences of the Algorithm <ref>.


Let the sequence {(u^k,v^k,p^k)}_k=0,1,2,… be generated by Algorithm <ref> and conditions of Theorem <ref> hold. By defining the auxiliary function F

    F(u,v,p,x):=ℒ_δ^φ(u,v,p)+
        7α^2θ^2K_2^4/2δu-x^2,

and

    w^k:=(u^k,v^k,p^k,u^k-1),

we have that

    F(w^k)-F(w^k+1)≥hu^k+1-u^k^2,

where

    h:=ν/2-( 3θ^2K_2^4/δ+7α^2θ^2K_2^4/2δ)>0.

If the sequence  {(u^k,v^k,p^k)}_k=0,1,2,… is bounded, then, it holds that

    lim_ku^k+1-u^k=lim_kv^k+1-v^k=lim_kp^k+1-p^k=0.





Based on Lemma <ref>, it is important to guarantee that the condition (<ref>) can be satisfied.
This fact can be reached if δ is large enough.  Fortunately, in the Algorithm <ref> the parameter δ can be fixed by the user. Thus, the parameter δ shall be chosen large enough to guarantee the convergence considering condition (<ref>).




If the nonconvex regulation function φ is coercive and

    δ>θ^2K_2^4,

the sequence {w^k}_k=0,1,2,… is bounded.





By combining the  conditions (<ref>) and (<ref>), we just need
that

    δ>max{1, 6θ^2K_2^4/ν+
        7α^2θ^2K_2^4/ν, θ^2K_2^4}.





Let  the sequence {(u^k,v^k,p^k)}_k=0,1,2,… be generated by Algorithm <ref>. Then, for any k, there exists γ>0 and s^k+1∈∂ F(w^k+1) such that

    s^k+1≤γ(u^k+1-u^k+u^k-u^k-1).




Now, we recall a  definition about the  limit point set ℳ  introduced in <cit.>, which
denotes the set of all the stationary points generated by the nonconvex IADMM. The specific mathematical definition of ℳ is given as follows.




Let {w^k}_k=0,1,2,… be generated by the nonconvex IADMM .
We define the set ℳ by

    ℳ:={x∈ℝ^N: ∃ an increasing sequence of integers
               {k_j}_j∈ℕ such that x^k_j→ x  as  j→∞}.





Let the sequence {(u^k,v^k,p^k)}_k=0,1,2,… be generated by Algorithm <ref>, F the auxiliary function defined in (<ref>) and suppose that  condition (<ref>) holds. Then, we have the following results.

(1) ℳ is nonempty and ℳ⊆crit(F).

(2) lim_kdist(w^k,ℳ)=0.

(3) The objective function F is finite and constant on ℳ.





§ NUMERICS


In this section, we illustrate the effectiveness of the proposed algorithm on different numerical blurred images with Gaussian blur.





All the programs have been written entirely in , and all the experiments are
implemented under Linux running on a desktop computer with an Intel Core
i5-2400S CPU (2.5 GHz) and 4 GB Memory. The FFT subroutines used in the algorithms are taken from the [] library.
As test problems we have selected  twelve images (see Figure <ref>), which include seven images from USC-SIPI image database[], two classical test images (Lena and cameraman), one text image from “El Quijote" book and two medical images. In order to obtain the blurred images, we use,  as it is common in literature, the blurring operator generated using a convolution with Gaussian kernel ( = 17 × 17,  =0,  = 7) and circular mapping on the edges of the image.






The proposed algorithm IADMM (Algorithm <ref>) is compared with the widely used augmented Lagrangian
methods (ADMM <cit.>) for image deblurring. We mainly consider two models, i.e., φ(·)=|·| and φ(·)=|·|^1/2. We call them as TV1 and TV(1/2) methods, respectively. Note that TV1 is a convex method, while TV(1/2) is nonconvex. In the tests we have considered (unless so indicated) for all the methods the value of the penalty parameter δ = 0.001 (a small one) and/or δ = 10 (a large one) just to see the behaviour of the IADMM algorithm.

The
performance of the deblurring algorithms is quantitatively measured by means of the objective function (Equations <ref> or <ref>), the Real Error as ·_2 of the difference between the original and deblurred images, the signal-to-noise ratio (SNR) <cit.>

    SNR(u,u^*)=10 ×log_10(u-u̅^2/u-u^*^2),

where  u and u^* denote the original image and the restored image, respectively, and u̅ represents the mean of the original image u, and the residual
((k, k+1) and _i(k, k+1) in the standard and inertial versions, respectively) as described in <cit.>:

    [ (k, k+1) :=  (u^k+1,   p^k+1)-(u^k,   p^k)/1+(u^k,   p^k),; _i(k, k+1) :=  (u^k+1,   p^k+1)-(u^k,  p^k)/1+(u^k,  p^k). ]

In the tests we do not provide CPU tests as all the algorithms have a very similar computation cost per iteration (mainly from the FFT routines). Therefore, there is almost no difference between pictures showing iterations or CPU cost, and the consumed CPU is basically proportional to the respective number of iterations.


On our first test we use the brain M0 image with a low value of δ = 0.001 and we show in Figure <ref> that the performances of the ADMM TV1 and TV(1/2) are quite similar in error and SNR. Therefore, in the rest of comparisons we will just consider the TV1 method. For the IADMM method the inertial parameter α in Algorithm <ref> is investigated firstly by using two different values (α = 0.2 and 0.5).
The main difference observed in these tests is that the nonconvex method is more unstable once reached the maximum precision (at this point the ADMM TV1 convex method seems to be the more stable with a quite smooth behaviour). The fastest convergence is observed using the IADMM (with the largest stepsize value α = 0.5) method, but when the maximum precision is obtained unstable behaviour appears. Therefore, using mainly the information provided by the  (Eq. (<ref>)), we provide a stop control criterion (in the same spirit as
<cit.>)
that stops the iterative process at the black points of Figure <ref> (d) in the tests. That is, we stop when the following condition is hold

    (k, k+1) < ε  or  (k-1, k) < (k, k+1).

In the first case the convergence till the desired tolerance error is obtained, while in the second case the algorithm has reached its stability limit and the residual grows, behaving later in an unstable way. Note that the residual is used in the stop criterion, as it uses known data from the iterations (and which does not depend on the original image that is unknown). We remark that the use of stop control techniques avoids the use of unnecessary iterations, and also to stop at the limit accuracy of the method. Also, from the pictures we observe that the IADMM algorithm provides
enough precision in a lower number of iterations. The larger α means faster method, but at the price of a more unstable method as it can be seen on Fig. <ref>(d). On that picture we observe that when the residual begins to behave chaotically, with sudden increases, it means that it is advisable to stop the iterative process as considered in the stop criterion (<ref>) (black dot points on Fig. <ref>(d)). With that criterion, the IADMM method seems to be an interesting option for fast deblurring problems.





To observe more clearly the influence of the parameter α in Algorithm <ref> we perform several tests on the brain M0 image on Figure <ref> for values α = 0.1, 0.2, 0.5, 0.8, 1.0 and 2.0. Note that this parameter plays a role similar to the stepsize (as it also occurs to the parameter δ), as it controls the perturbation at each step. A large value will provide, when the method works, a quite fast method, but on the other hand it makes the method more unstable. In fact, from the plot <ref>(b) we observe that in this case it would be optimal to use the parameter value α = 0.8 in combination with the stop criterion, giving the maximum precision in just 28 iterations. Besides, it is shown that after the values selected by the stop criterion the residual begins to oscillate among values that provides similar error but that generates an unstable behaviour giving rise to an increment of the error in subsequent iterations (this instability is delayed when the parameter α decreases, what is expected because the increment is smaller, as the vertical lines connecting the error and residual plots show).

The influence of the penalty parameter δ is also quite relevant, but a detailed analysis is out of the scope of this paper. On the Figure <ref> we show the evolution of the residual using two values of α = 0.2, 0.5
and several values of the parameter δ = 0.001, 0.01, 0.1, 1 and 10. We observe that low values of the penalty parameter δ gives a lower residual, but the error is lower for large δ providing a faster convergence, and it has a big effect on the empirical performance of the methods as shown in <cit.>, but it remains to study optimal combinations of the parameters δ and α and suitable criteria for an automatic selection (this will be part of the next steps in our study of these methods).



On Figure <ref>, we present the original medium resolution brain M0 image, the blurred one using, as indicated, a convolution with Gaussian kernel, and the results of the IADMM deblurring images using two error tolerances (ε=5 × 10^3, 5 × 10^4) in the stop criterion (<ref>). We can see that in both cases the quality of the recovered image is visually good.



  

A more detailed analysis is shown on Table <ref>, where we present for the IADMM (α=0.5 and 0.2) and the ADMM (TV1) methods (we do not show results for the TV(1/2) as they are quite similar to those of the TV1 ones), the iteration number, real error (), , the residual () and the efficiency rates of the IADMM (α=0.5) method over the two other ones,  and , for different values of the tolerance ε and using a small value of δ=0.001. In all the methods we have used the same values of the parameters and a similar stop criterion (<ref>) with the respective definition of the residual (Eqn. (<ref>)). From the tests, the inertial IADMM nonconvex methods are the fastest, as expected, but also with a more unstable behaviour, also as expected due to the nonconvex version of the methods. Note that from Figure <ref> the ADMM TV1 and TV(1/2) perform similarly for low δ, but if we make simulations small differences appear in medium-high resolution, obtaining for the M1 image the TV1 method 113 iterations and a ratio 2.26, but the TV(1/2)  99 iterations and a ratio  1.98; and for the H0 image the TV1 uses 110 iterations and a ratio  2.75, but the TV(1/2)   81 iterations and a ratio  2.03. That is, the nonconvex TV(1/2) may perform better in some circumstances, but the instability may also appear. In case of using larger values of the parameter δ, as the convergence theorems suggest, all the methods are much faster, giving a
better performance as shown on Table <ref> for δ = 10. Finally, we present on Table <ref> some results obtained by changing the value of the parameter δ and with the fixed value of the tolerance ε = 0.01. Again, the IADMM (α=0.5) method presents the best performance. Therefore, globally, the inertial IADMM version seems to be an interesting option, being the fastest one.

As above commented, how to select optimal combinations of the parameters δ and α, the use on other types of blur, development of suitable criteria for automatic selection of all the parameters and so on, remains a research issue, but these goals are out of the scope of this article.




§ CONCLUSION


In this paper, a more efficient  nonconvex inertial alternating minimization algorithm (IADMM)
is developed for solving the Total Variation model for  image deblurring. The proposed scheme is based on using the inertial proximal strategy on nonconvex ADMM methods. By the using the Kurdyka-Łojasiewicz property, we prove the convergence of the algorithm under several reasonable assumptions.
 Numerical experiments demonstrate that
the proposed algorithm overall outperforms the widely used augmented Lagrangian ADMM methods, being a fast option, although it can be unstable for high precision. This instability is easily controlled via a suitable stop control criteria.



§ APPENDIX A: PROOF OF LEMMA <REF>


Direct computation yields

    𝒯^* 𝒯=(
                      [ (I_N⊗ D)(I_N⊗ D)^*                   ;                    (D⊗ I_N)(D⊗ I_N)^*;                    ]).

Then, we have
 
    λ_min( 𝒯^*)^2=λ_min( 𝒯𝒯^*)=
       min{λ_min( (I_N⊗ D)(I_N⊗ D)^*), λ_min((D⊗ I_N)(D⊗ I_N)^* )}.

With Lemma 7.2 in <cit.>, we obtain

    λ_min( (I_N⊗ D)(I_N⊗ D)^*)
           =λ_min((D⊗ I_N)(D⊗ I_N)^* )=λ_min(DD^*).

Noting DD^* is a symmetric tridiagonal matrix

    DD^*=(
               [  2 -1      ; -1  2  ⋱   ;     ⋱  ⋱ -1;       -1  2;    ]),

and using the result in <cit.>,

    λ_min(DD^*)=2+2cos(π-π/N)
           =4cos^2(π/2-π/2N)=4sin^2(π/2N),

we obtain
λ_min(𝒯^*)≥ 2sin(π/2N).
Now, taking into account that 𝒯^* is full row-rank,
𝒯^* x≥λ_min(𝒯^*)x≥ 2sin(π/2N)x.



§ APPENDIX B: PROOF OF LEMMA <REF>

As the point u^k+1 is the minimization of ℒ_δ^φ(u,v^k+1,p^k), then the optimization condition yields

    K^*(K u^k+1-f)-𝒯^*p^k+δ𝒯^*(𝒯 u^k+1-v^k+1)=0.

Therefore, we can derive

    K^*(K u^k+1-f)-𝒯^*p^k+1=0,

and replacing k+1 by k

    K^*(K u^k-f)-𝒯^*p^k=0.

Subtracting (<ref>) and (<ref>), and using Lemma <ref>,

    p^k-p^k+1≤θ𝒯^*(p^k-p^k+1)
         =θK^*K(u^k+1-u^k≤θK_2^2·u^k+1-u^k.




§ APPENDIX C: PROOF OF LEMMA <REF>

Note that v^k+1 is the minimizer of ℒ_δ^φ(u^k,v,p^k) with respect to the variable v. Then we have:

    ℒ_δ^φ(u^k,v^k+1,p^k)≤ℒ_δ^φ(u^k,v^k,p^k).

In (<ref>), we have obtained

    ℒ_δ^φ(u^k+1,v^k+1,p^k)   +ν/2u^k+1-u^k^2≤ℒ_δ^φ(u^k,v^k+1,p^k).

By direct calculations, we obtain

    ℒ_δ^φ(u^k+1,v^k+1,p^k+1)
       =ℒ_δ^φ(u^k+1,v^k+1,p^k)+⟨p^k-p^k+1,𝒯 u^k+1-v^k+1⟩
       =ℒ_δ^φ(u^k+1,v^k+1,p^k)+1/δp^k-p^k+1_2^2
       ≤ℒ_δ^φ(u^k+1,v^k+1,p^k)+2/δp^k-p^k_2^2+2/δp^k-p^k+1_2^2
       ≤ℒ_δ^φ(u^k+1,v^k+1,p^k)+2α^2/δp^k-p^k-1_2^2+2/δp^k-p^k+1_2^2
       ≤ℒ_δ^φ(u^k+1,v^k+1,p^k)+2θ^2K_2^4/δu^k+1-u^k^2
                + 2α^2θ^2K_2^4/δu^k-u^k-1^2.

and

    ℒ_δ^φ(u^k,v^k,p^k)=ℒ_δ^φ(u^k,v^k,p^k)+⟨ p^k-p^k,𝒯 u^k-v^k⟩
       =ℒ_δ^φ(u^k,v^k,p^k)+1/δ⟨ p^k-p^k,p^k+1-p^k⟩
       ≤ℒ_δ^φ(u^k,v^k,p^k)+p^k-p^k^2/2δ+p^k+1-p^k^2/2δ
       ≤ℒ_δ^φ(u^k,v^k,p^k)+3α^2p^k-p^k-1^2/2δ+p^k+1-p^k^2/δ
       ≤ℒ_δ^φ(u^k,v^k,p^k)+ θ^2K_2^4/δu^k+1-u^k^2
                +3α^2θ^2K_2^4/2δu^k-u^k-1^2.

Combining the above equations, we obtain

    ℒ_δ^φ(x^k+1,y^k+1,p^k+1)+ν/2u^k+1-u^k^2
         ≤ℒ_δ^φ(x^k,y^k,p^k)+ 3θ^2K_2^4/δu^k+1-u^k^2 
               +7α^2θ^2K_2^4/2δu^k-u^k-1^2.

Thus, taking into account the definition of F (<ref>), we have

    F(w^k)-F(w^k+1)
              ≥[ν/2-( 3θ^2K_2^4/δ+7α^2θ^2K_2^4/2δ)]u^k+1-u^k^2.

Now, using the condition (<ref>) and the definition (<ref>) of h, we obtain the descent result given in Eq. (<ref>).

The boundness of {(u^k,v^k,p^k)}_k=0,1,2,… implies the boundness of {w^k}_k=0,1,2,…, and the continuity of F implies that {F(w^k)}_k=0,1,2,… is bounded. From the above proven result Eq. (<ref>),  F(w^k) is decreasing. Thus, the sequence {F(w^k)}_k=0,1,2,… is convergent, i.e., lim_k[F(w^k)-F(w^k+1)]=0.  Therefore, from (<ref>), we have

    lim_ku^k+1-u^k_2≤lim_k√(F(w^k)-F(w^k+1)/h)=0,

which indicates  lim_k p^k+1-p^k=0 due to (<ref>).  The scheme of updating p^k+1 gives us lim_k v^k+1-v^k=0.





§ APPENDIX D: PROOF OF LEMMA <REF>

From (<ref>), we derive

    p^k^2≤θ^2K_2^4·Ku^k-f^2.

By direct calculations, we obtain

    ℒ_δ^φ(u^k,v^k,p^k)= 1/2Ku^k-f^2+σv^k_φ
            -⟨ p^k,𝒯 u^k-v^k⟩+δ/2𝒯 u^k-v^k^2
       =1/2Ku^k-f^2+σv^k_φ +δ/2𝒯 u^k-v^k-p^k/δ^2-p^k^2/2δ
       ≥1/2Ku^k-f^2+σv^k_φ+δ/2𝒯 u^k-v^k-p^k/δ^2 
            -θ^2K_2^4/2δKu^k-f^2
       =(1/2-θ^2K_2^4/2δ)Ku^k-f^2 +δ/2𝒯 u^k-v^k-p^k/δ^2+σv^k_φ.

Thus, we have

    F(w^k)   ≥(1/2-θ^2K_2^4/2δ)Ku^k-f^2
           +δ/2𝒯 u^k-v^k-p^k/δ^2+σv^k_φ.

From Lemma <ref>, {F(w^k)}_k=0,1,2,… is bounded.
This means, using the definition of F(w^k), the boundedness of sequences {v^k_φ}_k=0,1,2,…,  {𝒯 u^k-v^k-p^k/δ}_k=0,1,2,… and {Ku^k-f}_k=0,1,2,…. Taking into account the coercivity of φ, {v^k}_k=0,1,2,… is bounded.  Now, taking into account (<ref>), {p^k}_k=0,1,2,… is bounded. Then, {(
     [ 𝒯; K;   ])
u^k}_k=0,1,2,… is bounded, and using the reversibility of (
     [ 𝒯; K;   ]), we obtain that {u^k}_k=0,1,2,… is bounded. Therefore, all the components of
the sequence {w^k}_k=0,1,2,… are bounded, and so the result is proved.



§ APPENDIX E: PROOF OF LEMMA <REF>

For the v^k+1 component, we  have

    -p^k+δ (𝒯 u^k-v^k+1)∈∂σv^k+1_φ.

Easy computations on the new defined term s^k+1_v give

    s^k+1_v:=(1+δ)(p^k+1-p^k) ∈∂_v F(w^k+1).

Direct calculations hold

    s^k+1_v   ≤(1+δ)p^k+1-p^k+(1+δ)p^k-p^k
       =(1+δ)p^k+1-p^k+α(1+δ)p^k-p^k-1
       ≤(1+δ)θK_2^2u^k+1-u^k  + α(1+δ)θK_2^2u^k-u^k-1.

Thus, we obtain that

    s^k+1_v≤γ_v(u^k+1-u^k+u^k-u^k-1),

where the parameter γ_v is given by

    γ_v=max{(1+δ)θK_2^2,  α(1+δ)θK_2^2}.



For the u^k+1 component, we  have

    𝒯^*p^k-δ𝒯^*(𝒯 u^k+1-v^k+1)=K^*(Ku^k+1-f).

With (<ref>), we have the new term s^k+1_u

    s^k+1_u   :=𝒯^*p^k-𝒯^*p^k+1+
        7α^2θ^2K_2^4/2δ(u^k+1-u^k)∈∂_u F(w^k+1).

It is easy to obtain

    s^k+1_u≤γ_u(u^k+1-u^k+u^k-u^k-1),

where 
    γ_u=max{θ𝒯_2K_2^2+7α^2θ^2K_2^4/2δ,  αθ𝒯_2K_2^2}.

Obviously, it holds

    s^k+1_p:=-𝒯 u^k+1+v^k+1∈∂_p F(w^k+1).

From the scheme, we can easily see

    s^k+1_p   =1/δp^k+1-p^k≤p^k+1-p^k/δ+αp^k-p^k-1/δ
       ≤γ_p(u^k+1-u^k+
        u^k-u^k-1),

where γ_p=max{θK_2^2/δ
,  αθK_2^2/δ}.
Therefore, we have

    ∇_xF(w^k+1)≤7α^2θ^2K_2^4/δ (d^k+1-d^k+u^k-u^k-1).



Letting s^k+1:=(s^k+1_v,s^k+1_u,s^k+1_p,∇_xF(w^k+1)), we have s^k+1∈∂ F(w^k+1). With all the  above equations (<ref>), (<ref>), (<ref>), and (<ref>), we obtain the global bound

    s^k+1   ≤(γ_u+γ_v+γ_p+7α^2θ^2K_2^4/δ)
           ×(u^k+1-u^k+u^k-u^k-1).

Denoting
γ=γ_u+γ_v+γ_p+7α^2θ^2K_2^4/δ,
and then, we finish the proof.




§ APPENDIX F: PROOF OF THEOREM <REF>

For any cluster point (u^*,v^*,p^*), there exists {k_j}_j=0,1,2,… such that lim_j(u^k_j,v^k_j,p^k_j)=(u^*,v^*,p^*). Then, from Lemma <ref>, we have
lim_j(u^k_j+1,v^k_j+1,p^k_j+1)=lim_j(u^k_j,v^k_j,p^k_j)=(u^*,v^*,p^*).
From the scheme of Algorithm <ref>, we have the following conditions

    -p^k_j+δ(𝒯 u^k_j-v^k_j+1)   ∈∂σv^k_j+1_φ,
    𝒯^*p^k_j-δ𝒯^*(𝒯 u^k_j+1-v^k_j+1)   =K^*(Ku^k_j+1-f),
    
     p^k_j+1   =p^k_j-δ(𝒯 u^k_j+1-v^k_j+1).

Letting j→+∞, with Proposition <ref>, we are then led to

    -p^*   ∈   ∂σv^*_φ,
    𝒯^*p^*   =   K^*(Ku^*-f),
    𝒯 u^*   =   v^*.

Finally, from Proposition <ref>, we obtain that (u^*,v^*,p^*) is a critical point of ℒ^φ_δ.






§ APPENDIX G: PROOF OF LEMMA <REF>

(1) From Lemma <ref> the sequences {(u^k,v^k,p^k)}_k=0,1,2,…  and {w^k}_k=0,1,2,… are bounded. Thus, ℳ is nonempty. Assume that w^*∈ℳ, then from the definition there exists
a subsequence w^k_i=(u^k_i,v^k_i,p^k_i,u^k_i-1)→ w^*=(u^*,v^*,p^*,u^∙). In that case, from Lemma <ref>, we also have u^k_i-1→ u^*. That is also u^∙=u^*, i.e., w^*=(u^*,v^*,p^*,u^*).   Besides, from Lemma <ref> and Lemma <ref>, we have s^k_i∈∂ F(w^k_i) and
s^k_i→0.  Finally, Proposition <ref> indicates that 0∈∂ F(w^*), i.e. w^*∈crit (F).

(2) This item follows as a consequence of the definition of the limit point set (Definition <ref>).

(3) In Eq. (<ref>), by taking limits on i, we obtain
limsup_iv^k_i_φ≤v^*_φ.
And with the closedness of φ, we have the other bound
   liminf_iv^k_i_φ≥v^*_φ.
That means
lim_iv^k_i_φ=v^*_φ.
Using the continuity of the rest of functions that defines F, we have
lim_iF(w^k_i)=F(w^*).
 On the other hand, Lemma <ref> indicates {F(w^k)}_k=0,1,2,… is decreasing. Noting the boundedness of {F(w^k)}_k=0,1,2,…, this sequence has a lower bound. Thus, {F(w^k)}_k=0,1,2,… is convergent. And then, we have
lim_kF(w^k)=lim_iF(w^k_i)=F(w^*).




§ APPENDIX H: PROOF OF THEOREM <REF>


From Lemma <ref>, F is constant on ℳ.
Let w^*=(u^*,v^*,p^*,x^*) be a stationary point of {w^k}_k=0,1,2,…. Then, from the definition of w^k we have
u^*=x^*.
Also from Lemma <ref>, we have dist(w^k,ℳ)<ε and F(w^k)<F(w^*)+η for any k>K, for some K. Hence, from Lemma <ref>, we have

    dist(0,∂ F(w^k))ρ'(F(w^k)-F(w^*))≥ 1,

that together with Lemma <ref> give us

    1/ρ'(F(w^k)-F(w^*))≤dist(0,∂ F(w^k))
           ≤γ(u^k+1-u^k+u^k-u^k-1).

Then, the concavity of ρ yields that

    F(w^k)-F(w^k+1)= F(w^k)-F(w^*)-[F(w^k+1)-F(w^*)]
         ≤ρ[F(w^k)-F(w^*)]-ρ[F(w^k+1)-F(w^*)]/ρ'[F(w^k)-F(w^*)]
         ≤γ(u^k+1-u^k+u^k-u^k-1)
         ×{ρ[F(w^k)-F(w^*)]-ρ[F(w^k+1)-F(w^*)]}.

Using Lemma <ref>, we have

    hu^k+1-u^k^2≤γ(u^k+1-u^k+u^k-u^k-1) 
         ×{ρ[F(w^k)-F(w^*)]-ρ[F(w^k+1)-F(w^*)]},

which is equivalent to

    h/γu^k+1-u^k≤2·1/2√(h/γ)√(u^k+1-u^k+u^k-u^k-1)
         ×√(ρ[F(w^k)-F(w^*)]-ρ[F(w^k+1)-F(w^*)]).

Using the Schwartz's inequality, we then derive

    h/γu^k+1-u^k≤h/4γ(u^k+1-u^k+u^k-u^k-1) 
          +{ρ[F(w^k)-F(w^*)]-ρ[F(w^k+1)-F(w^*)]}.

Summing (<ref>) from K to K+j yields that

    [ h/2γ∑_k=K^K+j-2u^k+1-u^k+3h/4γu^K+j+1-u^K+j;  ≤ρ[F(w^K)-F(w^*)]
       -φ[F(w^K+j+1)-F(w^*)];                           +h/4γu^K-1-u^K-2. ]

Letting j→+∞, and applying Lemma <ref>, we have
h/2γ∑_k=K^+∞u^k+1-u^k<+∞.
Then, {u^k}_k=0,1,2,… is convergent. And noting that u^* is a stationary point of {u^k}_k=0,1,2,… we also have that {u^k}_k=0,1,2,…  converges to u^*. With (<ref>), the convergence of {u^k}_k=0,1,2,… indicates the convergence of {p^k}_k=0,1,2,…. And then the sequence {p̂^k}_k=0,1,2,… is also convergent. With the identity p^k+1=p^k-δ(𝒯 u^k+1-v^k+1), {v^k}_k=0,1,2,… is convergent.
That means   (u^k,v^k,p^k) converges to some point (u^*,v^*,p^*). Finally, from Theorem <ref>, we have that (u^*,v^*,p^*) is a critical point of ℒ^φ_δ.







10

Hansen:2006
P. C. Hansen, J. G. Nagy, and D. P. O'Leary, Deblurring Images: Matrices,
  Spectra, and Filtering.
Philadelphia, PA, USA: Society for Industrial and Applied
  Mathematics, 2006.

lucy1974iterative
L. B. Lucy, “An iterative technique for the rectification of observed
  distributions,” The Astronomical Journal, vol. 79, p. 745, 1974.

rudin1992nonlinear
L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise
  removal algorithms,” Physica D: Nonlinear Phenomena, vol. 60, no. 1-4,
  pp. 259–268, 1992.

chen2015inertial
C. Chen, R. H. Chan, S. Ma, and J. Yang, “Inertial proximal ADMM for
  linearly constrained separable convex optimization,” SIAM Journal on
  Imaging Sciences, vol. 8, no. 4, pp. 2239–2267, 2015.

wang2008new
Y. Wang, J. Yang, W. Yin, and Y. Zhang, “A new alternating minimization
  algorithm for total variation image reconstruction,” SIAM Journal on
  Imaging Sciences, vol. 1, no. 3, pp. 248–272, 2008.

chen2018segmenting
F. Chen, H. Zhou, C. Grecos, and P. Ren, “Segmenting oil spills from blurry
  images based on alternating direction method of multipliers,” IEEE
  Journal of Selected Topics in Applied Earth Observations and Remote Sensing,
  vol. 11, no. 6, pp. 1858–1873, 2018.

hintermuler2013nonconvex
M. Hintermüller and T. Wu, “Nonconvex TV^q-models in image
  restoration: Analysis and a trust-region regularization–based
  superlinearly convergent solver,” SIAM Journal on Imaging Sciences,
  vol. 6, no. 3, pp. 1385–1415, 2013.

hintermuller2014smoothing
M. Hintermüller and T. Wu, “A smoothing descent method for nonconvex
  TV^q-models,” in Efficient Algorithms for Global Optimization
  Methods in Computer Vision, pp. 119–133, Springer, 2014.

xu2011image
L. Xu, C. Lu, Y. Xu, and J. Jia, “Image smoothing via l_0 gradient
  minimization,” ACM Transactions on Graphics (TOG), vol. 30, no. 6,
  p. 174, 2011.

Nikolova2008
M. Nikolova, M. K. Ng, S. Zhang, and W.-K. Ching, “Efficient reconstruction of
  piecewise constant images using nonsmooth nonconvex minimization,” SIAM
  Journal on Imaging Sciences, vol. 1, no. 1, pp. 2–25, 2008.

boyd2011distributed
S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed
  optimization and statistical learning via the alternating direction method of
  multipliers,” Foundations and Trends in Machine Learning, vol. 3,
  no. 1, pp. 1–122, 2011.

eckstein1994some
J. Eckstein and M. Fukushima, “Some reformulations and applications of the
  alternating direction method of multipliers,” in Large scale
  optimization, pp. 115–134, Springer, 1994.

gabay1976dual
D. Gabay and B. Mercier, “A dual algorithm for the solution of nonlinear
  variational problems via finite element approximation,” Computers &
  Mathematics with Applications, vol. 2, no. 1, pp. 17–40, 1976.

boley2013local
D. Boley, “Local linear convergence of the alternating direction method of
  multipliers on quadratic or linear programs,” SIAM Journal on
  Optimization, vol. 23, no. 4, pp. 2183–2207, 2013.

deng2012global
W. Deng and W. Yin, “On the global and linear convergence of the generalized
  alternating direction method of multipliers,” Journal of Scientific
  Computing, pp. 1–28, 2012.

he2012non
B. He and X. Yuan, “On non-ergodic convergence rate of Douglas–Rachford
  alternating direction method of multipliers,” Numerische Mathematik,
  vol. 130, no. 3, pp. 567–577, 2012.

he20121
B. He and X. Yuan, “On the o(1/n) convergence rate of the Douglas-Rachford
  alternating direction method,” SIAM Journal on Numerical Analysis,
  vol. 50, no. 2, pp. 700–709, 2012.

hong2012linear
M. Hong and Z.-Q. Luo, “On the linear convergence of the alternating direction
  method of multipliers,” Mathematical Programming, vol. 162,
  pp. 165–199, Mar 2017.

polyak1964some
B. T. Polyak, “Some methods of speeding up the convergence of iteration
  methods,” USSR Computational Mathematics and Mathematical Physics,
  vol. 4, no. 5, pp. 1–17, 1964.

banert2015backward
S. Banert and R. I. Boţ, “Backward penalty schemes for monotone
  inclusion problems,” Journal of Optimization Theory and Applications,
  vol. 166, no. 3, pp. 930–948, 2015.

bot2015penalty
R. I. Bot and E. R. Csetnek, “Penalty schemes with inertial effects for
  monotone inclusion problems,” Optimization, vol. 66, no. 6,
  pp. 965–982, 2017.

chen2015general
C. Chen, S. Ma, and J. Yang, “A general inertial proximal point algorithm for
  mixed variational inequality problem,” SIAM Journal on Optimization,
  vol. 25, no. 4, pp. 2120–2142, 2015.

ochs2014ipiano
P. Ochs, Y. Chen, T. Brox, and T. Pock, “iPiano: Inertial proximal algorithm
  for nonconvex optimization,” SIAM Journal on Imaging Sciences, vol. 7,
  no. 2, pp. 1388–1419, 2014.

sun2018noner
T. Sun, P. Yin, D. Li, C. Huang, L. Guan, and H. Jiang, “Non-ergodic
  convergence analysis of heavy-ball algorithms,” in AAAI Conference on
  Artificial Intelligence, 2019.

bot2014inertial
R. I. Bot and E. R. Csetnek, “An inertial alternating direction method of
  multipliers,” arXiv:1404.4582, 2014.

li2014splitting
G. Li and T. Pong, “Global convergence of splitting methods for nonconvex
  composite optimization,” SIAM Journal on Optimization, vol. 25, no. 4,
  pp. 2434–2460, 2015.

hong2016convergence
M. Hong, Z.-Q. Luo, and M. Razaviyayn, “Convergence analysis of alternating
  direction method of multipliers for a family of nonconvex problems,” SIAM Journal on Optimization, vol. 26, no. 1, pp. 337–364, 2016.

sun2017iteratively
T. Sun, H. Jiang, and L. Cheng, “Iteratively linearized reweighted alternating
  direction method of multipliers for a class of nonconvex problems,” arXiv:1709.00483, 2017.

sun2017alternating
T. Sun, P. Yin, L. Cheng, and H. Jiang, “Alternating direction method of
  multipliers with difference of convex functions,” Advances in
  Computational Mathematics, pp. 1–22, 2017.

wang2014convergence
F. Wang, Z. Xu, and H.-K. Xu, “Convergence of alternating direction method
  with multipliers for non-convex composite problems,” arXiv:1410.8625,
  2014.

wang2015global
Y. Wang, W. Yin, and J. Zeng, “Global convergence of ADMM in nonconvex
  nonsmooth optimization,” Journal of Scientific Computing, vol. 78,
  no. 1, pp. 29–63, 2019.

sun2017convergence
T. Sun, H. Jiang, L. Cheng, and W. Zhu, “A convergence frame for inexact
  nonconvex and nonsmooth algorithms and its applications to several
  iterations,” arXiv preprint arXiv:1709.04072, 2017.

xu2016empirical
Z. Xu, S. De, M. Figueiredo, C. Studer, and T. Goldstein, “An empirical study
  of admm for nonconvex problems,” arXiv preprint arXiv:1612.03349,
  2016.

taylor2016training
G. Taylor, R. Burmeister, Z. Xu, B. Singh, A. Patel, and T. Goldstein,
  “Training neural networks without gradients: A scalable ADMM approach,”
  in International Conference on Machine Learning, pp. 2722–2731, 2016.

mordukhovich2006variational
B. S. Mordukhovich, Variational analysis and generalized differentiation
  I: Basic theory, vol. 330.
Springer Science & Business Media, 2006.

nesterov2004introductory
Y. Nesterov, Introductory lectures on convex optimization, vol. 87.
Springer Science & Business Media, 2004.

rockafellar2009variational
R. T. Rockafellar and R. J.-B. Wets, Variational analysis, vol. 317.
Springer Science & Business Media, 2009.

rockafellar2015convex
R. T. Rockafellar, Convex analysis.
Princeton university press, 2015.

Lojasiewicz63
S. Łojasiewicz, Les Équations aux Dérivées Partielles, ch. Une
  propriété topologique des sous-ensembles analytiques réels, pp. 8–89.
Paris: Éditions du Centre National de la Recherche Scientifique,
  1963.

Kurdyka98
K. Kurdyka, “On gradients of functions definable in o-minimal structures,”
  Annales de l'Institut Fourier, vol. 48, pp. 769–783, 1998.

attouch2013convergence
H. Attouch, J. Bolte, and B. F. Svaiter, “Convergence of descent methods for
  semi-algebraic and tame problems: proximal algorithms, forward–backward
  splitting, and regularized Gauss–Seidel methods,” Mathematical
  Programming, vol. 137, no. 1-2, pp. 91–129, 2013.

bolte2014proximal
J. Bolte, S. Sabach, and M. Teboulle, “Proximal alternating linearized
  minimization for nonconvex and nonsmooth problems,” Mathematical
  Programming, vol. 146, no. 1-2, pp. 459–494, 2014.

halko2011finding
N. Halko, P.-G. Martinsson, and J. A. Tropp, “Finding structure with
  randomness: Probabilistic algorithms for constructing approximate matrix
  decompositions,” SIAM Review, vol. 53, no. 2, pp. 217–288, 2011.

alizadeh1998primal
F. Alizadeh, J.-P. A. Haeberly, and M. L. Overton, “Primal-dual interior-point
  methods for semidefinite programming: convergence rates, stability and
  numerical results,” SIAM Journal on Optimization, vol. 8, no. 3,
  pp. 746–768, 1998.

link
T. I. Lakoba, “The eigenvalues of tridiagonal matrices,” Technical Note
  of the Norwegian University of Science and Technology, Department of
  Mathematical Sciences, NoteTMA4205, 2009.



























