Inertial nonconvex alternating minimizations for the image deblurring

TL;DR

This paper introduces a novel inertial nonconvex ADMM algorithm for image deblurring that demonstrates improved efficiency and convergence properties over traditional convex methods, validated through numerical experiments.

Contribution

It extends the classical ADMM with an inertial technique to handle nonconvex TV regularization, providing convergence guarantees under certain conditions.

Findings

The proposed IADMM converges under specific parameter assumptions.

Numerical simulations show improved performance over traditional ADMM.

The method effectively reconstructs blurred images with better numerical stability.

Abstract

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-{\L}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.

Tables

Table 1. TABLE I: Basic notations, where x , y 𝑥 𝑦 x,y stand for points, C 𝐶 C stands for a set, A 𝐴 A stands for matrix

$\parallel \cdot \parallel$ stands for $\parallel \cdot {\parallel }_{L^2\left(\Omega \right)}$ or $\parallel \cdot {\parallel }_{L^2\left({\Omega }^{\ast }\right)}$ ($L_2$ norm)
dist$(x,C):={\min }_{y\in C}\Vert x-y\Vert$	${\Vert A\Vert }_2:={\max }_{\Vert x\Vert =1}\Vert Ax\Vert$
$\parallel \cdot {\parallel }_{\phi }:=\parallel \cdot {\parallel }_{\phi \left({\Omega }^{\ast }\right)}$	$\otimes$ stands for the Kronecker product
${𝒞}^1$ stands for the function class whose derivatives are continuous
for a matrix A, rank(A) stands for its rank, $\text{Null}\left(A\right):=\{x\mid Ax=\text{0}\}$
${\lambda }_{\min }\left(A\right)$ stands for the minimum eigenvalue of $A$

Table 2. TABLE II: Deblurring results for a small value of δ 𝛿 \delta ( δ = 0.001 𝛿 0.001 \delta=0.001 ) using different values of the tolerance ε 𝜀 \varepsilon . Iteration number ( I1, I2, I3 ), Real error ( ERROR ), SNR , Residual ( RES ) and efficiency rates I2/I1 and I3/I1 using the IADMM ( α = 0.5 𝛼 0.5 \alpha=0.5 and 0.2 0.2 0.2 ) and the ADMM (TV1) methods applied to all the 12 test images ( IMG ).

	IADMM ($\alpha =0.5$)				IADMM ($\alpha =0.2$)					ADMM(TV1)
IMG	I1	ERROR	SNR	RES	I2	ERROR	SNR	RES	I2/I1	I3	ERROR	SNR	RES	I3/I1
L0	4	17.4	11.1	4.59e-03	6	17.4	11.1	4.62e-03	1.50	7	17.6	11.0	4.80e-03	1.75
L1	5	16.7	10.0	4.67e-03	8	16.5	10.1	4.60e-03	1.60	9	16.8	9.9	4.95e-03	1.80
L2	2	19.4	8.7	4.36e-03	2	19.9	8.5	4.69e-03	1.00	2	20.2	8.3	4.88e-03	1.00
L3	2	16.6	13.2	3.01e-03	2	17.1	13.0	3.25e-03	1.00	2	17.9	12.6	3.67e-03	1.00
M0	8	30.6	13.8	4.50e-03	13	30.4	13.8	4.67e-03	1.62	16	30.6	13.8	4.88e-03	2.00
M1	9	28.5	13.7	4.33e-03	14	28.7	13.6	4.86e-03	1.56	18	28.5	13.7	4.78e-03	2.00
M2	6	36.6	12.8	3.48e-03	7	40.1	12.0	4.82e-03	1.17	9	39.9	12.1	4.86e-03	1.50
M3	7	35.5	12.5	4.34e-03	11	35.6	12.4	4.50e-03	1.57	13	36.3	12.3	4.97e-03	1.86
H0	8	71.4	10.3	4.50e-03	13	71.0	10.3	4.65e-03	1.62	16	71.4	10.3	4.84e-03	2.00
H1	9	69.2	6.3	4.47e-03	14	69.7	6.2	4.96e-03	1.56	18	69.2	6.3	4.89e-03	2.00
H2	7	74.2	3.3	4.16e-03	10	76.2	3.0	4.87e-03	1.43	13	75.5	3.1	4.77e-03	1.86
H3	6	77.4	1.5	4.08e-03	8	80.8	1.1	4.87e-03	1.33	11	78.9	1.3	4.48e-03	1.83

Table 3. TABLE III: Deblurring results for a large value of δ 𝛿 \delta ( δ = 10 𝛿 10 \delta=10 ) using the tolerance ε = 0.01 𝜀 0.01 \varepsilon=0.01 . Iteration number ( I1, I2, I3 ), Real error ( ERROR ), SNR , Residual ( RES ) and efficiency rates I2/I1 and I3/I1 using the IADMM ( α = 0.5 𝛼 0.5 \alpha=0.5 and 0.2 0.2 0.2 ) and the ADMM (TV1) methods applied to all the 12 test images ( IMG ).

	IADMM ($\alpha =0.5$)				IADMM ($\alpha =0.2$)					ADMM(TV1)
IMG	I1	ERROR	SNR	RES	I2	ERROR	SNR	RES	I2/I1	I3	ERROR	SNR	RES	I3/I1
L0	5	3.0	26.5	4.72e-01	5	2.8	27.1	4.31e-01	1.00	10	6.2	20.0	9.56e-02	2.00
L1	3	3.3	24.1	5.38e-01	5	3.2	24.3	4.62e-01	1.67	10	5.8	19.2	9.48e-02	3.33
L2	1	3.8	22.9	8.22e-01	3	3.7	23.0	5.13e-01	3.00	10	9.3	15.0	9.68e-02	10.00
L3	6	1.4	34.6	2.41e-01	7	1.3	35.5	2.00e-01	1.17	10	7.0	20.7	9.86e-02	1.67
M0	3	6.2	27.6	5.79e-01	5	6.0	28.0	4.98e-01	1.67	10	7.0	26.5	9.74e-02	3.33
M1	7	4.3	30.2	4.40e-01	7	3.9	31.0	4.12e-01	1.00	10	5.4	28.1	9.78e-02	1.43
M2	4	3.0	34.5	1.50e-01	4	3.0	34.7	1.00e-01	1.00	10	20.6	17.8	9.90e-02	2.50
M3	1	14.6	20.2	8.94e-01	1	14.5	20.2	8.68e-01	1.00	10	31.7	13.4	9.64e-02	10.00
H0	5	13.8	24.5	4.96e-01	5	12.6	25.3	4.44e-01	1.00	10	15.2	23.7	9.52e-02	2.00
H1	3	16.6	18.7	5.73e-01	5	16.2	18.9	4.98e-01	1.67	10	20.4	16.9	9.47e-02	3.33
H2	5	15.7	16.8	5.11e-01	5	14.6	17.4	4.57e-01	1.00	10	20.7	14.4	9.29e-02	2.00

Table 4. TABLE IV: Deblurring results for several values of parameter δ 𝛿 \delta using the tolerance ε = 0.01 𝜀 0.01 \varepsilon=0.01 on M0 image. Iteration number ( I1, I2, I3 ), Real error ( ERROR ), SNR , Residual ( RES ) and efficiency rates I2/I1 and I3/I1 using the IADMM ( α = 0.5 𝛼 0.5 \alpha=0.5 and 0.2 0.2 0.2 ) and the ADMM (TV1) methods.

	IADMM ($\alpha =0.5$)				IADMM ($\alpha =0.2$)					ADMM(TV1)
$\delta$	I1	ERROR	SNR	RES	I2	ERROR	SNR	RES	I2/I1	I3	ERROR	SNR	RES	I3/I1
0.001	6	33.5	13.0	7.40e-03	8	35.3	12.5	8.70e-03	1.33	9	34.0	12.8	7.50e-03	1.50
0.01	4	19.1	17.9	8.88e-03	4	21.2	16.9	9.82e-03	1.00	5	20.0	17.4	7.68e-03	1.25
0.10	3	11.2	22.5	2.39e-02	4	11.4	22.3	2.31e-02	1.33	5	10.8	22.8	2.32e-02	1.67
1.00	1	8.0	25.4	2.25e-01	2	7.5	26.0	2.03e-01	2.00	3	7.3	26.3	2.20e-01	3.00
10.00	3	6.2	27.6	5.79e-01	5	6.0	28.0	4.98e-01	1.67	10	7.0	26.5	9.74e-02	3.33

Equations

$$\tilde{f}=\tilde{K}(\tilde{u})+e,$$

$$\underset{\tilde{u}}{\min }\left(\frac{1}{2}\Vert \tilde{K}\tilde{u}-\tilde{f}{\Vert }_{L^2(\Omega )}^2+\sigma \Vert \nabla \tilde{u}{\Vert }_{L^1({\Omega }^{\ast })}\right),$$

$$\underset{\tilde{u}}{\min }\left(\frac{1}{2}\Vert \tilde{K}\tilde{u}-\tilde{f}{\Vert }_{L^2(\Omega )}^2+\sigma \Vert \nabla \tilde{u}{\Vert }_{\phi ({\Omega }^{\ast })}\right).$$

$${\tilde{\mathcal{L}}}_{\delta }^{\phi }(\tilde{u},\tilde{v},\tilde{p}):=\frac{1}{2}\Vert \tilde{K}\tilde{u}-\tilde{f}{\Vert }_{L^2(\Omega )}^2+\sigma \Vert \tilde{v}{\Vert }_{\phi ({\Omega }^{\ast })}$$

$$-⟨\tilde{p},\nabla \tilde{u}-\tilde{v}⟩+\frac{\delta }{2}\Vert \nabla \tilde{u}-\tilde{v}{\Vert }_{L^2({\Omega }^{\ast })}^2,$$

$$\displaystyle\left\{\begin{array}[]{lcl}\tilde{v}^{k+1}&=&\textrm{arg}\min_{\tilde{v}}\mathcal{\tilde{L}}_{\delta}^{1}(\tilde{u}^{k},\tilde{v},\tilde{p}^{k})\\[5.0pt] \tilde{u}^{k+1}&=&\textrm{arg}\min_{\tilde{u}}\mathcal{\tilde{L}}_{\delta}^{1}(\tilde{u},\tilde{v}^{k+1},\tilde{p}^{k})\\[5.0pt] \tilde{p}^{k+1}&=&\tilde{p}^{k}-\delta(\nabla\tilde{u}^{k+1}-\tilde{v}^{k+1})\\ \end{array}\right.$$

$$\partial J(x):=\{u\in {\mathbb{R}}^N:\exists x^k\to x,u^k\to u,\text{such that}$$

$$\underset{y\ne x^k}{\lim }\underset{y\to x^k}{\inf }\frac{J(y)-J(x^k)-⟨u^k,y-x^k⟩}{\Vert y-x^k{\Vert }_2}\ge 0\}.$$

$$\text{0}\in \partial J(x),$$

$$\displaystyle\left\{\begin{array}[]{lcl}-\bar{p}&\in&\sigma\partial\|\bar{v}\|_{\varphi},\\[5.0pt] \nabla^{*}\bar{p}&=&K^{*}(K\bar{u}-f),\\[5.0pt] \nabla\bar{u}&=&\bar{v}.\end{array}\right.$$

$$\displaystyle\mathcal{S}_{J}(x)\in\textrm{arg}\min_{y}\bigg{\{}J(y)+\frac{1}{2}\|x-y\|_{2}^{2}\bigg{\}}.$$

$${\rho }^{\prime }(J(x)-J(\hat{x}))\cdot \text{dist}(0,\partial J(x))\ge 1.$$

$$\displaystyle\nabla(\tilde{u})=\left(\begin{array}[]{c}I_{N}\otimes D\\ D\otimes I_{N}\\ \end{array}\right)\cdot\tilde{u},$$

$$D=\left(\begin{array}[]{ccccc}1&-1&&&\\ &1&\ddots&&\\ &&\ddots&-1&\\ &&&1&-1\\ \end{array}\right)_{(N-1)\times N}.$$

$$\mathcal{T}:=\left(\begin{array}[]{cc}I_{N}\otimes D&\\ &D\otimes I_{N}\\ \end{array}\right),\quad u:=\left(\begin{array}[]{c}u_{1}\in\mathbb{R}^{N^{2}}\\ u_{2}\in\mathbb{R}^{N^{2}}\\ \end{array}\right).$$

$$\begin{array}[]{l}\textrm{rank}(\mathcal{T})=\textrm{rank}(I_{N}\otimes D)+\textrm{rank}(D\otimes I_{N})\\ =\textrm{rank}(I_{N})\cdot\textrm{rank}(D)+\textrm{rank}(D)\cdot\textrm{rank}(I_{N})=2N(N-1),\end{array}$$

$$\Vert {\mathcal{T}}^{\ast }x\Vert \ge \frac{\Vert x\Vert }{\theta },$$

$$\underset{u_1=u_2}{\min }\left(\frac{1}{2}\Vert \tilde{K}{u}_1-\tilde{f}{\Vert }^2+\sigma \Vert \mathcal{T}u{\Vert }_{\phi }\right).$$

$$\underset{u}{\min }\left(\frac{1}{2}\Vert \tilde{K}{u}_1-\tilde{f}{\Vert }^2+\sigma \Vert \mathcal{T}u{\Vert }_{\phi }+\frac{{\beta }^2}{2}\Vert u_1-u_2{\Vert }^2\right),$$

$$\underset{u\in {\mathbb{R}}^{2N^2}}{\min }\left(\frac{1}{2}\Vert Ku-f{\Vert }^2+\sigma \Vert \mathcal{T}u{\Vert }_{\phi }\right),$$

$${\mathcal{L}}_{\delta }^{\phi }(u,v,p):=\frac{1}{2}\Vert Ku-f{\Vert }^2$$

$$+\sigma \Vert v{\Vert }_{\phi }-⟨p,\mathcal{T}u-v⟩+\frac{\delta }{2}\Vert \mathcal{T}u-v{\Vert }^2.$$

$$\displaystyle\left\{\begin{array}[]{lcl}(\widehat{u}^{k},\widehat{v}^{k},\widehat{p}^{k})&=&(u^{k},v^{k},p^{k})\\[4.0pt] &&+\alpha(u^{k}-u^{k-1},v^{k}-v^{k-1},p^{k}-p^{k-1}),\\[4.0pt] v^{k+1}&=&\textrm{arg}\min_{v}\mathcal{L}_{\delta}^{\varphi}(u^{k},v,\widehat{p}^{k}),\\[4.0pt] u^{k+1}&=&\textrm{arg}\min_{u}\mathcal{L}_{\delta}^{\varphi}(u,v^{k+1},\widehat{p}^{k}),\\[4.0pt] p^{k+1}&=&\widehat{p}^{k}-\delta(\mathcal{T}u^{k+1}-v^{k+1}),\end{array}\right.$$

$$\displaystyle v^{k+1}=\mathcal{S}_{\frac{\sigma}{\delta}\|\cdot\|_{\varphi(\Omega^{*})}}\bigg{(}\mathcal{T}u^{k}-\frac{\widehat{p}^{k}}{\delta}\bigg{)},$$

$$[{\mathcal{S}}_{\frac{\sigma }{\delta }\Vert \cdot {\Vert }_{\phi }}(v){]}_{i,j}={\mathcal{S}}_{\frac{\sigma }{\delta }\phi }(v_{i,j}).$$

$$K^{\ast }(K{u}^{k+1}-f)-{\mathcal{T}}^{\ast }{\hat{p}}^k+\delta {\mathcal{T}}^{\ast }(\mathcal{T}{u}^{k+1}-v^{k+1})=\text{0}.$$

$$u^{k+1}=(K^{\ast }K+\delta {\mathcal{T}}^{\ast }\mathcal{T}{)}^{-1}(K^{\ast }f+\delta {\mathcal{T}}^{\ast }{v}^{k+1}+{\mathcal{T}}^{\ast }{\hat{p}}^k).$$

$${\mathcal{L}}_{\delta }^{\phi }(u^k,v^{k+1},{\hat{p}}^k)-{\mathcal{L}}_{\delta }^{\phi }(u^{k+1},v^{k+1},{\hat{p}}^k)$$

$$\ge \frac{{\lambda }_{\min }(K^{\ast }K+\delta {\mathcal{T}}^{\ast }\mathcal{T})}{2}\Vert u^{k+1}-u^k{\Vert }^2\ge \frac{\nu }{2}\Vert u^{k+1}-u^k{\Vert }^2,$$

$$\begin{gathered} \Vert (K^{\ast }K+{\mathcal{T}}^{\ast }\mathcal{T}{)}^{-1}{\Vert }_2 \\ \le b\sqrt{\frac{2}{\pi }}\underset{i=1,2,\dots ,M}{\max }\Vert (K^{\ast }K+{\mathcal{T}}^{\ast }\mathcal{T}{)}^{-1}{w}_i{\Vert }_2 \end{gathered}$$

$$\delta >\max \left\{1,\frac{6{\theta }^2\Vert K{\Vert }_2^4}{\nu }+\frac{7{\alpha }^2{\theta }^2\Vert K{\Vert }_2^4}{\nu }\right\}$$

Full text

I 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 $\Omega\equiv\mathbb{R}^{N^{2}}$ be a two-dimensional index set representing the image domain, $\tilde{u}\in\Omega$ be the original image, $\tilde{f}$ be the observed image, and $\tilde{K}$ be a linear blurring operator. Then, the blurred image can be written [1] as

[TABLE]

where $e$ is an unknown additive noise vector. In this paper, the blurring operator $\tilde{K}$ is assumed to be known, otherwise, one will deal with the blind image deblurring problem [2], in which $\tilde{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 [3] for solving the deblurring problem (although other methods are also used). The TV method can be presented as

[TABLE]

being $\sigma>0$ a parameter and where $\nabla$ is the gradient operator, $\Omega^{*}:=\nabla(\Omega)$, and the norms $\|\tilde{v}\|_{L^{1}(\Omega^{*})}:=\sum_{i,j\in\Omega^{*}}|\tilde{v}_{i,j}|$, $\|\tilde{u}\|_{L^{2}(\Omega)}:=\sqrt{\sum_{i,j\in\Omega}|\tilde{u}_{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 [4, 5, 6], but in recent years, some nonconvex methods have been developed [7, 8, 9]. The use of a suitable nonconvex and nondifferentiable function allows possibly a smaller number of measurements than the convex one in compressed sensing [7]. In [10] 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\|_{\varphi(\Omega^{*})}:=\sum_{i,j\in\Omega^{*}}\varphi(v_{i,j})$ that uses the nonconvex regulation function $\varphi$, and that we refer to as the $\ell^{q}$ semi-norm ($0<q<1$). Therefore, the general nonconvex deblurring model is presented as

[TABLE]

Many efficient numerical algorithms have been developed for solving the TV regularization problem. One of the most efficient methods for the convex problem (2) is the Alternating Direction Method of Multipliers (ADMM) algorithm [11, 12, 13]. In the general case (convex and nonconvex cases depend of the function $\varphi$), the method is constructed by introducing an auxiliary variable $\tilde{v}$, which actually represents $\nabla\tilde{u}$, to reformulate (3) into a composite optimization problem with linear constraints. The augmented Lagrange dual function is then

[TABLE]

where $\delta>0$ is a parameter, and the norm $\|X\|_{L^{2}(\Omega^{*})}:=\sqrt{\sum_{i,j\in\Omega^{*}}|X_{i,j}|^{2}}$. If $\varphi(\cdot)=|\cdot|_{1}$, we use the notation $\mathcal{\tilde{L}}_{\delta}^{1}(\tilde{u},\tilde{v},\tilde{p})$ for representing $\mathcal{\tilde{L}}_{\delta}^{\varphi}(\tilde{u},\tilde{v},\tilde{p})$. Now, the standard convex ADMM method ($\varphi(\cdot)=|\cdot|_{1}$) for the deblurring problem can be presented as

[TABLE]

The earlier analyses of convergence and performance of the ADMM algorithms directly depended on the existing results of ADMM framework [14, 15, 16, 17, 18]. More recently, motivated by acceleration techniques proposed in [19], inertial algorithms have been proposed for many areas such as (distributed) optimization and imaging sciences in references [20, 21, 22, 23, 24]. The ideas of the inertial strategy have been also applied to ADMM in [4, 25]; 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 [26, 27, 28, 29, 30, 31, 32, 33, 34]. 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 $\varphi$ is nonconvex, we consider a new inertial scheme for the image deblurring model (3). 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 $\nabla$ fails to be injective (although the numerical performance may be good in some cases). Considering this, we first modify the model (3), 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 $\delta$ 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 $\delta$.

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

II 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 I.

II-A Subdifferentials

We collect several definitions as well as some useful properties in variational and convex analysis (see the monographes [35, 36, 37, 38]). For any matrix $A$, we define $A^{*}$ to be the adjoint of $A$.

Definition 1

Let $J:\mathbb{R}^{N}\rightarrow(-\infty,+\infty]$ be a proper and lower semicontinuous function. The (limiting) subdifferential, or simply the subdifferential, of $J$ at $x\in\mathbb{R}^{N}$, written as $\partial J(x)$, is defined as

[TABLE]

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 [38]. We call $J$ is strongly convex with constant $a\in\mathbb{R}^{+}$ if for any $x,y\in\textrm{dom}(J)$ and any $v\in\partial J(x)$, it holds that $J(y)\geq J(x)+\langle v,y-x\rangle+\frac{a}{2}\|y-x\|_{2}^{2}.$ And $J$ is called as $L$-gradient continuous (Lipschitz) with constant $L>0$ if $\|\nabla J(x)-\nabla J(y)\|_{2}\leq L\|x-y\|_{2}.$ Noting the closedness of the subdifferential, we have the following simple proposition.

Proposition 1

If $v^{k}\in\partial J(x^{k})$, $\lim_{k}v^{k}=v$ and $\lim_{k}x^{k}=x$, then we have $v\in\partial J(x).$

Definition 2

A necessary condition for $x\in\mathbb{R}^{N}$ to be a minimizer of $J(x)$ is

[TABLE]

which is also sufficient when $J$ is convex. A point that satisfies (9) is called (limiting) critical point. The set of critical points of $J(x)$ is denoted by $\textrm{crit}(J)$.

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

Proposition 2

If $(\bar{u},\bar{v},\bar{p})$ is a critical point of $\mathcal{L}^{\varphi}_{\delta}$ whose definition given in (III), it must hold that

[TABLE]

Finally, the proximal map of $J$ is defined as

[TABLE]

Note that $J$ can be nonconvex. If $J$ is convex, $\mathcal{S}_{J}$ is a point-to-point operator; otherwise, it may be point-to-set.

II-B Kurdyka-Łojasiewicz property

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

Definition 3** ([41, 42])**

(a) The function $J:\mathbb{R}^{N}\rightarrow(-\infty,+\infty]$ is said to have the Kurdyka-Łojasiewicz property at $\widehat{x}\in\textrm{dom}(\partial J)$ if there exist $\eta\in(0,+\infty]$, a neighborhood $U$ of $\widehat{x}$ and a continuous concave function $\rho:[0,\eta)\rightarrow\mathbb{R}^{+}$ such that

$\rho(0)=0$. 2. 2.

$\rho$* is $\mathcal{C}^{1}$ on $(0,\eta)$.* 3. 3.

For all $s\in(0,\eta)$, $\rho^{\prime}(s)>0$. 4. 4.

For all $x$ in $U\bigcap\{x|J(\widehat{x})<J(x)<J(\widehat{x})+\eta\}$, the Kurdyka-Łojasiewicz inequality holds

[TABLE]

(b) Proper lower semicontinuous functions which satisfy the Kurdyka-Łojasiewicz inequality at each point of $\textrm{dom}(\partial J)$ are called KŁ functions.

Remark 1

There are large sets of functions that are KŁ functions [41].

Lemma 1** ([42])**

Let $J:\mathbb{R}^{N}\rightarrow\mathbb{R}$ be a proper lower semi-continuous function and $\Pi$ be a compact set. If $J$ is a constant on $\Pi$ and $J$ satisfies the KŁ property at each point on $\Pi$, then there exists a concave function $\rho$ satisfying the four properties given in Definition 3, and constants $\eta,\varepsilon>0$ such that for any $\widehat{x}\in\Pi$ and any $x$ satisfying that $\textrm{dist}(x,\Pi)<\varepsilon$ and $f(\widehat{x})<f(x)<f(\widehat{x})+\eta$, it holds that

[TABLE]

III 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 $\Omega=[1,2,\ldots,N]\times[1,2,\ldots,N]$ (equivalent to space $\mathbb{R}^{N^{2}}$) as the two-dimensional index set representing the image domain. In this case, the image variable constrained on $\Omega$ is actually a $N\times N$ matrix. We use the symbol $\tilde{u}$ 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

[TABLE]

where $I_{N}$ the identity matrix and $D$ the banded matrix

[TABLE]

If we directly apply the inertial ADMM, the convergence is hard to be proved as $\nabla$ fails to be injective. Therefore, we need to modify the image deblurring model (3). To that goal, we define

[TABLE]

Obviously, we have $\mathcal{T}\in\mathbb{R}^{2N(N-1)\times 2N^{2}}$. Noting that

[TABLE]

and thus, $\mathcal{T}$ is injective. The following technical lemma focuses on giving a lower bound for the operator $\mathcal{T}^{*}$.

Lemma 2

For any $x\in\mathbb{R}^{2N(N-1)}$, it holds that

[TABLE]

where \theta=1/(2\sin(\frac{\pi}{2N})\big{)}.

Then, the image deblurring model (3) is equivalent to

[TABLE]

Instead, we consider its extended penalty form

[TABLE]

where $\beta>0$ is a large weight parameter. Therefore, we apply the nonconvex inertial ADMM to

[TABLE]

where K=\left(\begin{array}[]{cc}\tilde{K}&\textbf{0}\\ \beta\mathbb{I}&-\beta\mathbb{I}\\ \end{array}\right) and f=\left(\begin{array}[]{c}\tilde{f}\\ \textbf{0}\\ \end{array}\right). This leads us to define the function

[TABLE]

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 [4], 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 $\varphi$ 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

[TABLE]

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

Now we can focus on rewriting the inertial scheme for the image deblurring model (30). First, we rearrange the minimization for $v^{k+1}$,

[TABLE]

being $S_{J}$ the proximal map of $J$ (14). For a matrix $v$, and indices $(i,j)\in\Omega^{*}$,

[TABLE]

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

[TABLE]

That is also

[TABLE]

Taking into account (30), (31) and (34) we propose a nonconvex inertial version of the IADMM algorithm (Algorithm 1).

Assumption 1: We assume that $\textrm{Null}(K_{0})\bigcap\textrm{Null}(\mathcal{T})=\{\textbf{0}\}$, where K_{0}:=\left(\begin{array}[]{cc}\tilde{K}&\textbf{0}\\ \mathbb{I}&-\mathbb{I}\\ \end{array}\right). And the minimum single value is given as $\nu$.

This hypothesis also indicates that the matrix \left(\begin{array}[]{c}\mathcal{T}\\ K_{0}\\ \end{array}\right) is reversible. Note that the rank of $\mathcal{T}$ is $2N(N-1)$. Then, the assumed hypothesis is easy to be satisfied.

We remark that $u^{k+1}$ is the minimizer of $\mathcal{L}_{\delta}^{\varphi}(u,v^{k+1},\widehat{p}^{k})$, which is strongly convex with $\lambda_{\min}(K^{*}K+\delta\mathcal{T}^{*}\mathcal{T})$. If we set $\delta,\rho>1$, then we have

[TABLE]

where we used the fact $\lambda_{\min}(K^{*}K+\delta\mathcal{T}^{*}\mathcal{T})\geq\lambda_{\min}(K^{*}_{0}K_{0}+\mathcal{T}^{*}\mathcal{T})$ when $\delta,\beta\geq 1$.

The following problem is what $\nu$ exactly is. In a real situation, the dimensions of $\mathcal{T}$ 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 $\theta$. If $K^{*}K+\mathcal{T}^{*}\mathcal{T}$ is reversible, it is easy to see that ${1}/{\nu}=\|(K^{*}K+\mathcal{T}^{*}\mathcal{T})^{-1}\|_{2}$. Then, if we obtain a bound $\|(K^{*}K+\mathcal{T}^{*}\mathcal{T})^{-1}\|_{2}\leq c$, we then have $\nu\geq{1}/{c}$. To that goal we employ a Lemma proposed in [43]:

Lemma 3** (Lemma 4.1, [43])**

For a fixed positive integer $M$ and a real number $b>1$, and given an independent family $\{w_{i}\}_{i=1,2,\ldots,M}$ of standard Gaussian vectors, we have that

[TABLE]

with probability at least $1-b^{-M}$.

Note that for computing $(K^{*}K+\mathcal{T}^{*}\mathcal{T})^{-1}w_{i}$ we just need several FFT and inverse FFT. Therefore, its computational cost is low ($\mathcal{O}(N\log N$), and the estimation of $\nu$ is very fast.

IV 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.

IV-A Main results

Theorem 1** (Stationary point convergence)**

Assume that the free parameter $\delta$ satisfies the condition

[TABLE]

with \theta=1/(2\sin(\frac{\pi}{2N})\big{)}, then any cluster point $(u^{*},v^{*},p^{*})$ is also a critical point of $\mathcal{L}^{\varphi}_{\delta}$.

Theorem 1 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 2.

Theorem 2** (Sequence convergence)**

Let condition (46) hold, and the auxiliary function $F$ (40) be a KŁ function. Then the sequence $\{(u^{k},v^{k},p^{k})\}_{k=0,1,2,\ldots}$, generated by Algorithm 1, converges to a critical point of $\mathcal{L}^{\varphi}_{\delta}$.

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,\ldots}$, where $w^{k}$ are composite points from $\{(u^{k},v^{k},p^{k})\}_{k=0,1,2,\ldots}$. An auxiliary function $F$ is also proposed. In Lemma 5, we prove a “sufficient descent condition” of the values of $F$ at $\{w^{k}\}_{k=0,1,2,\ldots}$, i.e.,

[TABLE]

where $\widehat{h}>0$ is a positive constant.

II. We prove a “relative error condition” of $\{w^{k}\}_{k=0,1,2,\ldots}$, i.e., there exists $v^{k+1}\in\partial F(w^{k+1})$ such that

[TABLE]

where $\gamma>0$ is a positive constant. Note that this condition is different from the “real” relative error condition proposed in [41].

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}-\widehat{v}^{k}\|^{2}$ and $\|p^{k+1}-\widehat{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}-\widehat{v}^{k}\|$ and $\|p^{k+1}-\widehat{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.

IV-B Technical lemmas

Lemma 4

Let the sequence $\{(u^{k},v^{k},p^{k})\}_{k=0,1,2,\ldots}$ be generated by Algorithm 1 to solve problem (3), then

[TABLE]

where $\theta$ is given in Lemma 2 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 1.

Lemma 5

Let the sequence $\{(u^{k},v^{k},p^{k})\}_{k=0,1,2,\ldots}$ be generated by Algorithm 1 and conditions of Theorem 1 hold. By defining the auxiliary function $F$

[TABLE]

and

[TABLE]

we have that

[TABLE]

where

[TABLE]

If the sequence $\{(u^{k},v^{k},p^{k})\}_{k=0,1,2,\ldots}$ is bounded, then, it holds that

[TABLE]

Remark 2

Based on Lemma 5, it is important to guarantee that the condition (36) can be satisfied. This fact can be reached if $\delta$ is large enough. Fortunately, in the Algorithm 1 the parameter $\delta$ can be fixed by the user. Thus, the parameter $\delta$ shall be chosen large enough to guarantee the convergence considering condition (36).

Lemma 6

If the nonconvex regulation function $\varphi$ is coercive and

[TABLE]

the sequence $\{w^{k}\}_{k=0,1,2,\ldots}$ is bounded.

Remark 3

By combining the conditions (36) and (45), we just need that

[TABLE]

Lemma 7

Let the sequence $\{(u^{k},v^{k},p^{k})\}_{k=0,1,2,\ldots}$ be generated by Algorithm 1. Then, for any $k$, there exists $\gamma>0$ and $s^{k+1}\in\partial F(w^{k+1})$ such that

[TABLE]

Now, we recall a definition about the limit point set $\mathcal{M}$ introduced in [42], which denotes the set of all the stationary points generated by the nonconvex IADMM. The specific mathematical definition of $\mathcal{M}$ is given as follows.

Definition 4

Let $\{w^{k}\}_{k=0,1,2,\ldots}$ be generated by the nonconvex IADMM . We define the set $\mathcal{M}$ by

[TABLE]

Lemma 8

Let the sequence $\{(u^{k},v^{k},p^{k})\}_{k=0,1,2,\ldots}$ be generated by Algorithm 1, $F$ the auxiliary function defined in (40) and suppose that condition (46) holds. Then, we have the following results.

(1) $\mathcal{M}$ is nonempty and $\mathcal{M}\subseteq\textrm{crit}(F)$.

(2) $\lim_{k}\textrm{dist}(w^{k},\mathcal{M})=0$.

(3) The objective function $F$ is finite and constant on $\mathcal{M}$.

V 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 C++, 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 fftw-3111http://www.fftw.org/ library. As test problems we have selected twelve images (see Figure 1), which include seven images from USC-SIPI image database222http://sipi.usc.edu/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 (KernelSize $=17\times 17$, KernelMu $=0$, KernelSigma $=7$) and circular mapping on the edges of the image.

The proposed algorithm IADMM (Algorithm 1) is compared with the widely used augmented Lagrangian methods (ADMM [22]) for image deblurring. We mainly consider two models, i.e., $\varphi(\cdot)=|\cdot|$ and $\varphi(\cdot)=|\cdot|^{\frac{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 $\delta=0.001$ (a small one) and/or $\delta=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 2 or 3), the Real Error as $\|\cdot\|_{2}$ of the difference between the original and deblurred images, the signal-to-noise ratio (SNR) [22]

[TABLE]

where $u$ and $u^{*}$ denote the original image and the restored image, respectively, and $\bar{u}$ represents the mean of the original image $u$, and the residual ($\texttt{res}(k,k+1)$ and ${\texttt{res}_{i}}(k,k+1)$ in the standard and inertial versions, respectively) as described in [22]:

[TABLE]

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 $\delta=0.001$ and we show in Figure 2 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 $\alpha$ in Algorithm 1 is investigated firstly by using two different values ($\alpha=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 $\alpha=0.5$) method, but when the maximum precision is obtained unstable behaviour appears. Therefore, using mainly the information provided by the residual (Eq. (50)), we provide a stop control criterion (in the same spirit as [22]) that stops the iterative process at the black points of Figure 2 (d) in the tests. That is, we stop when the following condition is hold

[TABLE]

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 $\alpha$ means faster method, but at the price of a more unstable method as it can be seen on Fig. 2(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 (51) (black dot points on Fig. 2(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 $\alpha$ in Algorithm 1 we perform several tests on the brain M0 image on Figure 3 for values $\alpha=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 $\delta$), 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 3(b) we observe that in this case it would be optimal to use the parameter value $\alpha=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 $\alpha$ 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 $\delta$ is also quite relevant, but a detailed analysis is out of the scope of this paper. On the Figure 4 we show the evolution of the residual using two values of $\alpha=0.2,0.5$ and several values of the parameter $\delta=0.001,0.01,0.1,1$ and $10$. We observe that low values of the penalty parameter $\delta$ gives a lower residual, but the error is lower for large $\delta$ providing a faster convergence, and it has a big effect on the empirical performance of the methods as shown in [33], but it remains to study optimal combinations of the parameters $\delta$ and $\alpha$ and suitable criteria for an automatic selection (this will be part of the next steps in our study of these methods).

On Figure 5, 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 ($\varepsilon=5\times 10^{3}$, $5\times 10^{4}$) in the stop criterion (51). We can see that in both cases the quality of the recovered image is visually good.

A more detailed analysis is shown on Table II, where we present for the IADMM ($\alpha=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 (ERROR), SNR, the residual (RES) and the efficiency rates of the IADMM ($\alpha=0.5$) method over the two other ones, I2/I1 and I3/I1, for different values of the tolerance $\varepsilon$ and using a small value of $\delta=0.001$. In all the methods we have used the same values of the parameters and a similar stop criterion (51) with the respective definition of the residual (Eqn. (50)). 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 2 the ADMM TV1 and TV(1/2) perform similarly for low $\delta$, 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 $\delta$, as the convergence theorems suggest, all the methods are much faster, giving a better performance as shown on Table III for $\delta=10$. Finally, we present on Table IV some results obtained by changing the value of the parameter $\delta$ and with the fixed value of the tolerance $\varepsilon=0.01$. Again, the IADMM ($\alpha=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 $\delta$ and $\alpha$, 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.

VI 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 2

Direct computation yields

[TABLE]

Then, we have

[TABLE]

With Lemma 7.2 in [44], we obtain

[TABLE]

Noting $DD^{*}$ is a symmetric tridiagonal matrix

[TABLE]

and using the result in [45],

[TABLE]

we obtain \lambda_{\min}(\mathcal{T}^{*})\geq 2\sin\bigg{(}\frac{\pi}{2N}\bigg{)}. Now, taking into account that $\mathcal{T}^{*}$ is full row-rank, \|\mathcal{T}^{*}x\|\geq\lambda_{\min}(\mathcal{T}^{*})\|x\|\geq 2\sin\bigg{(}\frac{\pi}{2N}\bigg{)}\|x\|.

Appendix B: proof of Lemma 4

As the point $u^{k+1}$ is the minimization of $\mathcal{L}_{\delta}^{\varphi}(u,v^{k+1},\widehat{p}^{k})$, then the optimization condition yields

[TABLE]

Therefore, we can derive

[TABLE]

and replacing $k+1$ by $k$

[TABLE]

Subtracting (60) and (61), and using Lemma 2,

[TABLE]

Appendix C: proof of Lemma 5

Note that $v^{k+1}$ is the minimizer of $\mathcal{L}_{\delta}^{\varphi}(u^{k},v,\widehat{p}^{k})$ with respect to the variable $v$. Then we have:

[TABLE]

In (III), we have obtained

[TABLE]

By direct calculations, we obtain

[TABLE]

and

[TABLE]

Combining the above equations, we obtain

[TABLE]

Thus, taking into account the definition of $F$ (40), we have

[TABLE]

Now, using the condition (36) and the definition (43) of $\widehat{h}$, we obtain the descent result given in Eq. (42).

The boundness of $\{(u^{k},v^{k},p^{k})\}_{k=0,1,2,\ldots}$ implies the boundness of $\{w^{k}\}_{k=0,1,2,\ldots}$, and the continuity of $F$ implies that $\{F(w^{k})\}_{k=0,1,2,\ldots}$ is bounded. From the above proven result Eq. (42), $F(w^{k})$ is decreasing. Thus, the sequence $\{F(w^{k})\}_{k=0,1,2,\ldots}$ is convergent, i.e., $\lim_{k}[F(w^{k})-F(w^{k+1})]=0$. Therefore, from (42), we have

[TABLE]

which indicates $\lim_{k}\|p^{k+1}-p^{k}\|=0$ due to (39). The scheme of updating $p^{k+1}$ gives us $\lim_{k}\|v^{k+1}-v^{k}\|=0$.

Appendix D: proof of Lemma 6

From (61), we derive

[TABLE]

By direct calculations, we obtain

[TABLE]

Thus, we have

[TABLE]

From Lemma 5, $\{F(w^{k})\}_{k=0,1,2,\ldots}$ is bounded. This means, using the definition of $F(w^{k})$, the boundedness of sequences $\{\|v^{k}\|_{\varphi}\}_{k=0,1,2,\ldots}$, $\{\mathcal{T}u^{k}-v^{k}-\frac{p^{k}}{\delta}\}_{k=0,1,2,\ldots}$ and $\{Ku^{k}-f\}_{k=0,1,2,\ldots}$. Taking into account the coercivity of $\varphi$, $\{v^{k}\}_{k=0,1,2,\ldots}$ is bounded. Now, taking into account (69), $\{p^{k}\}_{k=0,1,2,\ldots}$ is bounded. Then, \bigg{\{}\left(\begin{array}[]{c}\mathcal{T}\\ K\\ \end{array}\right)u^{k}\bigg{\}}_{k=0,1,2,\ldots} is bounded, and using the reversibility of \left(\begin{array}[]{c}\mathcal{T}\\ K\\ \end{array}\right), we obtain that $\{u^{k}\}_{k=0,1,2,\ldots}$ is bounded. Therefore, all the components of the sequence $\{w^{k}\}_{k=0,1,2,\ldots}$ are bounded, and so the result is proved.

Appendix E: Proof of Lemma 7

For the $v^{k+1}$ component, we have

[TABLE]

Easy computations on the new defined term $s^{k+1}_{v}$ give

[TABLE]

Direct calculations hold

[TABLE]

Thus, we obtain that

[TABLE]

where the parameter $\gamma_{v}$ is given by

[TABLE]

For the $u^{k+1}$ component, we have

[TABLE]

With (73), we have the new term $s^{k+1}_{u}$

[TABLE]

It is easy to obtain

[TABLE]

where

[TABLE]

Obviously, it holds

[TABLE]

From the scheme, we can easily see

[TABLE]

where \gamma_{p}=\max\bigg{\{}\frac{\theta\|K\|_{2}^{2}}{\delta},\,\frac{\alpha\theta\|K\|_{2}^{2}}{\delta}\bigg{\}}. Therefore, we have

[TABLE]

Letting $s^{k+1}:=(s^{k+1}_{v},s^{k+1}_{u},s^{k+1}_{p},\nabla_{x}F(w^{k+1}))$, we have $s^{k+1}\in\partial F(w^{k+1})$. With all the above equations (72), (75), (76), and (77), we obtain the global bound

[TABLE]

Denoting $\gamma=\gamma_{u}+\gamma_{v}+\gamma_{p}+\frac{7\alpha^{2}\theta^{2}\|K\|_{2}^{4}}{\delta},$ and then, we finish the proof.

Appendix F: Proof of Theorem 1

For any cluster point $(u^{*},v^{*},p^{*})$, there exists $\{k_{j}\}_{j=0,1,2,\ldots}$ such that $\lim_{j}(u^{k_{j}},v^{k_{j}},p^{k_{j}})=(u^{*},v^{*},p^{*})$. Then, from Lemma 5, we have $\lim_{j}(u^{k_{j}+1},v^{k_{j}+1},p^{k_{j}+1})=\lim_{j}(\widehat{u}^{k_{j}},\widehat{v}^{k_{j}},\widehat{p}^{k_{j}})=(u^{*},v^{*},p^{*}).$ From the scheme of Algorithm 1, we have the following conditions

[TABLE]

Letting $j\rightarrow+\infty$, with Proposition 1, we are then led to

[TABLE]

Finally, from Proposition 2, we obtain that $(u^{*},v^{*},p^{*})$ is a critical point of $\mathcal{L}^{\varphi}_{\delta}$.

Appendix G: Proof of Lemma 8

(1) From Lemma 6 the sequences $\{(u^{k},v^{k},p^{k})\}_{k=0,1,2,\ldots}$ and $\{w^{k}\}_{k=0,1,2,\ldots}$ are bounded. Thus, $\mathcal{M}$ is nonempty. Assume that $w^{*}\in\mathcal{M}$, then from the definition there exists a subsequence $w^{k_{i}}=(u^{k_{i}},v^{k_{i}},p^{k_{i}},u^{k_{i}-1})\rightarrow w^{*}=(u^{*},v^{*},p^{*},u^{\bullet})$. In that case, from Lemma 5, we also have $u^{k_{i}-1}\rightarrow u^{*}$. That is also $u^{\bullet}=u^{*}$, i.e., $w^{*}=(u^{*},v^{*},p^{*},u^{*})$. Besides, from Lemma 7 and Lemma 5, we have $s^{k_{i}}\in\partial F(w^{k_{i}})$ and $s^{k_{i}}\rightarrow\mathbf{0}$. Finally, Proposition 1 indicates that $\mathbf{0}\in\partial F(w^{*})$, i.e. $w^{*}\in\textrm{crit}(F)$.

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

(3) In Eq. (64), by taking limits on $i$, we obtain $\lim\sup_{i}\|v^{k_{i}}\|_{\varphi}\leq\|v^{*}\|_{\varphi}.$ And with the closedness of $\varphi$, we have the other bound $\lim\inf_{i}\|v^{k_{i}}\|_{\varphi}\geq\|v^{*}\|_{\varphi}$. That means $\lim_{i}\|v^{k_{i}}\|_{\varphi}=\|v^{*}\|_{\varphi}.$ Using the continuity of the rest of functions that defines $F$, we have $\lim_{i}F(w^{k_{i}})=F(w^{*}).$ On the other hand, Lemma 5 indicates $\{F(w^{k})\}_{k=0,1,2,\ldots}$ is decreasing. Noting the boundedness of $\{F(w^{k})\}_{k=0,1,2,\ldots}$, this sequence has a lower bound. Thus, $\{F(w^{k})\}_{k=0,1,2,\ldots}$ is convergent. And then, we have $\lim_{k}F(w^{k})=\lim_{i}F(w^{k_{i}})=F(w^{*}).$

Appendix H: Proof of Theorem 2

From Lemma 8, $F$ is constant on $\mathcal{M}$. Let $w^{*}=(u^{*},v^{*},p^{*},x^{*})$ be a stationary point of $\{w^{k}\}_{k=0,1,2,\ldots}$. Then, from the definition of $w^{k}$ we have $u^{*}=x^{*}.$ Also from Lemma 8, we have $\textrm{dist}(w^{k},\mathcal{M})<\varepsilon$ and $F(w^{k})<F(w^{*})+\eta$ for any $k>K$, for some $K$. Hence, from Lemma 1, we have

[TABLE]

that together with Lemma 7 give us

[TABLE]

Then, the concavity of $\rho$ yields that

[TABLE]

Using Lemma 7, we have

[TABLE]

which is equivalent to

[TABLE]

Using the Schwartz’s inequality, we then derive

[TABLE]

Summing (Appendix H: Proof of Theorem 2) from $K$ to $K+j$ yields that

[TABLE]

Letting $j\rightarrow+\infty$, and applying Lemma 5, we have $\frac{\widehat{h}}{2\gamma}\sum_{k=K}^{+\infty}\|u^{k+1}-u^{k}\|<+\infty.$ Then, $\{u^{k}\}_{k=0,1,2,\ldots}$ is convergent. And noting that $u^{*}$ is a stationary point of $\{u^{k}\}_{k=0,1,2,\ldots}$ we also have that $\{u^{k}\}_{k=0,1,2,\ldots}$ converges to $u^{*}$. With (39), the convergence of $\{u^{k}\}_{k=0,1,2,\ldots}$ indicates the convergence of $\{p^{k}\}_{k=0,1,2,\ldots}$. And then the sequence $\{\hat{p}^{k}\}_{k=0,1,2,\ldots}$ is also convergent. With the identity $p^{k+1}=\widehat{p}^{k}-\delta(\mathcal{T}u^{k+1}-v^{k+1})$, $\{v^{k}\}_{k=0,1,2,\ldots}$ is convergent. That means $(u^{k},v^{k},p^{k})$ converges to some point $(u^{*},v^{*},p^{*})$. Finally, from Theorem 1, we have that $(u^{*},v^{*},p^{*})$ is a critical point of $\mathcal{L}^{\varphi}_{\delta}$.