L0TV: A Sparse Optimization Method for Impulse Noise Image Restoration

TL;DR

This paper introduces a novel sparse optimization method, TV-PADMM, for impulse noise image restoration that outperforms existing methods by effectively solving a non-convex -norm based TV restoration problem.

Contribution

It proposes a new TV-PADMM algorithm that addresses non-convex -norm TV restoration, with proven convergence, for improved impulse noise removal in images.

Findings

TV-PADMM outperforms state-of-the-art methods in impulse noise removal.

The method effectively handles non-convex -norm optimization problems.

Experimental results show superior denoising and deblurring performance.

Abstract

Total Variation (TV) is an effective and popular prior model in the field of regularization-based image processing. This paper focuses on total variation for removing impulse noise in image restoration. This type of noise frequently arises in data acquisition and transmission due to many reasons, e.g. a faulty sensor or analog-to-digital converter errors. Removing this noise is an important task in image restoration. State-of-the-art methods such as Adaptive Outlier Pursuit(AOP) \cite{yan2013restoration}, which is based on TV with $\ell_{02}$-norm data fidelity, only give sub-optimal performance. In this paper, we propose a new sparse optimization method, called $\ell_0TV$-PADMM, which solves the TV-based restoration problem with $\ell_0$-norm data fidelity. To effectively deal with the resulting non-convex non-smooth optimization problem, we first reformulate it as an equivalent biconvex Mathematical Program with Equilibrium Constraints (MPEC), and then solve it using a proximal Alternating Direction Method of Multipliers (PADMM). Our $\ell_0TV$-PADMM method finds a desirable solution to the original $\ell_0$-norm optimization problem and is proven to be convergent under mild conditions. We apply $\ell_0TV$-PADMM to the problems of image denoising and deblurring in the presence of impulse noise. Our extensive experiments demonstrate that $\ell_0TV$-PADMM outperforms state-of-the-art image restoration methods.

Figures

Tables

Table 1. Table I: Data Fidelity Models

Data Fidelity Function	Noise and References
${\ell }_2(\mathrm{𝐊𝐮},𝐛)={\Vert \mathrm{𝐊𝐮}-𝐛\Vert }_2^2$	add. Gaussian noise [47, 14]
${\ell }_1(\mathrm{𝐊𝐮},𝐛)={\Vert \mathrm{𝐊𝐮}-𝐛\Vert }_1$	add. Laplace noise [60, 23]
${\ell }_{\infty }(\mathrm{𝐊𝐮},𝐛)={\Vert \mathrm{𝐊𝐮}-𝐛\Vert }_{\infty }$	add. uniform noise [22, 51]
${\ell }_p(\mathrm{𝐊𝐮},𝐛)=⟨\mathrm{𝐊𝐮}-𝐛\odot \log \left(\mathrm{𝐊𝐮}\right),\mathrm{𝟏}⟩$	mul. Poisson noise [36, 49]
${\ell }_g(\mathrm{𝐊𝐮},𝐛)=⟨\log \left(\mathrm{𝐊𝐮}\right)+𝐛\odot \frac{1}{\mathrm{𝐊𝐮}},\mathrm{𝟏}⟩$	mul. Gamma noise [3, 53]
${\ell }_r(\mathrm{𝐊𝐮},𝐛)=⟨\log \left(\mathrm{𝐊𝐮}\right)+𝐛\odot 𝐛\odot \frac{1}{\mathrm{𝟐}𝐊𝐮},\mathrm{𝟏}⟩$	mul. Rayleigh noise [48, 2]
${\ell }_{02}(\mathrm{𝐊𝐮},𝐛)={\Vert \mathrm{𝐊𝐮}-𝐛+𝐳\Vert }_2^2,s.t.{\Vert 𝐳\Vert }_0\le k$	mixed Gaussian impulse noise [59]
${\ell }_0(\mathrm{𝐊𝐮},𝐛)={\Vert \mathrm{𝐊𝐮}-𝐛\Vert }_0$	add./mul. impulse noise [ours]

Table 2. Table II: Regularization Models

Regularization Function	Description and References
${\Omega }_{\text{tik}}(𝐠,𝐡)={\sum }_{i=1}^n{𝐠}_i^2+{𝐡}_i^2$	Tikhonov-like [1]
${\Omega }_{{\text{tv}}_2}(𝐠,𝐡)={\sum }_{i=1}^n{\left({𝐠}_i^2+{𝐡}_i^2\right)}^{\frac{1}{2}}$	Isotropic [47, 53]
${\Omega }_{{\text{tv}}_1}(𝐠,𝐡)={\sum }_{i=1}^n\left|{𝐠}_i\right|+\left|{𝐡}_i\right|$	Anisotropic [50, 60]
${\Omega }_{\text{stv}}(𝐠,𝐡)={\sum }_{i=1}^n{\left({𝐠}_i^2+{𝐡}_i^2+{\epsilon }^2\right)}^{\frac{1}{2}}$	smooth TV [18, 51]
${\Omega }_{\text{pot}}(𝐠,𝐡)={\sum }_{i=1}^n{\left|{𝐠}_i\right|}_0+{\left|{𝐡}_i\right|}_0$	Potts model [56, 57]
$\begin{array}{c}{\Omega }_{\text{hub}}(𝐠,𝐡)={\sum }_{i=1}^n\phi ({𝐠}_i;{𝐡}_i),\\ \phi ({𝐠}_i;{𝐡}_i)=\{\begin{array}{cc}\epsilon \parallel {𝐠}_i;{𝐡}_i{\parallel }_2^2/2;\parallel {𝐠}_i;{𝐡}_i{\parallel }_2\le 1/\epsilon & \\ \parallel {𝐠}_i;{𝐡}_i{\parallel }_2-\epsilon /2;\text{otherwise}& \end{array}\end{array}$	Huber-Like [46]

Table 3. Table III: General denoising problems. The results separated by ‘/’ are S ​ N ​ R 0 𝑆 𝑁 subscript 𝑅 0 SNR_{0} , S ​ N ​ R 1 𝑆 𝑁 subscript 𝑅 1 SNR_{1} and S ​ N ​ R 2 𝑆 𝑁 subscript 𝑅 2 SNR_{2} , respectively. The 1 s ​ t superscript 1 𝑠 𝑡 1^{st} , 2 n ​ d superscript 2 𝑛 𝑑 2^{nd} , and 3 r ​ d superscript 3 𝑟 𝑑 3^{rd} best results are colored with red , blue and green , respectively.

	BM3D	${\ell }_{1}TV$-$SBM$	$MFM$	$TSM$	${\ell }_{02}TV$-$AOP$	${\ell }_{P}TV$-$PADMM$	${\ell }_{0}TV$-$PDA$	${\ell }_{0}TV$-$PADM$
Random-Value Impulse Noise
walkbridge+10%	93/7.1/11.0	95/12.3/15.6	92/7.7/12.3	95/11.8/12.9	96/12.8/16.6	95/12.1/13.8	97/14.1/16.9	97/13.8/15.9
walkbridge+30%	76/3.7/7.1	89/8.6/11.0	82/6.1/10.3	85/5.8/7.8	89/8.4/12.1	89/7.8/11.5	91/9.6/12.8	91/9.5/11.9
walkbridge+50%	59/2.2/4.3	76/4.9/5.7	67/4.1/7.0	69/2.7/4.8	76/5.4/8.1	79/5.4/8.7	84/7.0/10.1	85/7.0/9.2
walkbridge+70%	42/1.0/1.9	56/2.0/1.7	45/2.0/3.3	50/1.3/2.2	53/2.5/4.0	59/3.0/5.0	65/4.0/6.2	76/5.1/7.0
walkbridge+90%	26/-0.1/-0.1	32/-0.2/-1.1	28/0.3/0.5	30/0.0/-0.0	31/0.4/0.8	30/0.4/0.8	34/0.7/1.3	57/2.7/3.9
pepper+10%	67/5.0/9.9	99/19.1/21.5	99/15.0/22.2	97/13.5/15.8	74/5.4/11.3	99/13.6/20.3	100/20.2/24.6	99/18.0/21.0
pepper+30%	55/3.7/7.0	96/12.3/13.6	96/11.4/16.3	87/6.3/9.5	72/5.2/10.7	98/12.0/16.8	98/15.1/19.7	98/14.6/18.3
pepper+50%	44/2.4/4.5	85/6.7/6.7	85/7.0/9.7	71/3.5/5.5	65/4.5/8.9	94/9.7/13.1	96/11.8/15.7	96/11.6/14.4
pepper+70%	33/1.2/2.1	63/2.8/2.1	59/3.1/4.4	52/1.6/2.4	51/2.7/4.7	79/5.2/6.2	84/6.8/8.9	93/9.0/11.4
pepper+90%	24/0.2/0.1	35/0.1/-1.0	30/0.6/0.6	31/0.3/0.1	28/0.7/1.1	35/0.9/1.0	39/1.3/1.7	76/4.2/4.8
mandrill+10%	74/3.3/6.0	89/8.1/9.0	92/6.9/6.9	93/9.6/9.6	84/3.7/7.4	93/9.6/9.6	95/11.1/11.5	95/10.8/10.3
mandrill+30%	63/2.0/3.6	83/5.9/6.6	76/3.8/5.9	83/4.7/4.9	73/3.0/5.5	85/5.8/6.8	87/6.8/7.4	86/6.4/6.5
mandrill+50%	50/1.1/2.2	73/3.6/3.7	65/2.9/4.6	69/2.0/3.4	61/2.2/4.0	74/3.6/5.0	77/4.6/5.6	78/4.4/4.6
mandrill+70%	36/0.4/0.8	57/1.4/0.6	51/1.5/2.4	52/0.9/1.5	47/1.2/2.2	62/2.3/3.4	64/2.9/3.9	70/3.1/3.5
mandrill+90%	28/-0.3/-0.6	36/-0.6/-1.9	37/0.2/0.4	34/-0.1/-0.4	33/0.1/0.3	39/0.5/0.9	42/0.8/1.2	58/1.9/2.5
lake+10%	92/6.9/12.5	98/16.9/21.3	96/11.3/17.7	97/14.0/15.0	97/8.7/16.1	98/14.3/19.2	98/17.2/21.1	98/16.7/19.5
lake+30%	75/4.3/8.1	93/11.3/13.9	91/9.3/14.4	86/7.1/10.0	92/7.9/13.9	95/10.5/15.0	95/12.7/16.7	95/12.0/14.3
lake+50%	58/2.6/4.9	79/6.5/7.2	71/5.9/9.4	69/3.7/5.9	78/6.2/10.2	88/8.3/11.7	91/10.0/13.7	90/9.5/11.5
lake+70%	41/1.3/2.3	54/2.9/2.6	42/2.5/4.1	47/1.8/2.8	43/2.8/4.6	60/4.7/7.0	68/5.8/8.6	84/7.4/9.0
lake+90%	24/0.3/0.3	26/0.5/-0.4	25/0.6/0.8	26/0.5/0.4	24/0.6/1.0	13/0.7/1.1	26/1.1/1.7	62/4.2/5.3
jetplane+10%	39/2.5/6.1	99/17.5/21.0	98/11.5/17.5	98/12.8/13.3	39/3.4/8.3	99/13.1/19.1	99/17.0/20.0	98/15.6/17.0
jetplane+30%	32/0.7/2.6	95/10.3/11.5	94/9.0/13.3	87/5.0/7.3	38/3.2/7.5	97/10.4/15.0	97/12.4/15.7	97/11.5/12.6
jetplane+50%	27/-0.6/-0.1	80/4.5/4.0	75/4.2/6.7	69/1.5/2.8	34/2.4/5.2	92/7.9/10.6	94/9.3/12.2	94/9.0/10.0
jetplane+70%	22/-1.7/-2.4	53/0.6/-0.7	42/0.2/0.9	47/-0.5/-0.5	23/-0.6/-0.3	67/3.2/4.8	74/4.4/6.4	90/6.7/7.4
jetplane+90%	18/-2.5/-4.1	25/-1.8/-3.6	25/-1.7/-2.5	26/-1.8/-2.9	18/-2.3/-3.4	14/-1.6/-2.2	26/-1.2/-1.5	74/3.4/3.7
Salt-and-Pepper Impulse Noise
walkbridge+10%	90/5.4/9.9	96/12.9/17.3	90/7.6/12.4	98/15.8/19.9	98/16.3/20.7	98/15.8/19.9	99/17.2/22.7	99/17.5/23.2
walkbridge+30%	71/3.0/4.5	94/10.4/14.3	83/6.3/9.8	96/11.7/16.4	94/10.5/15.2	96/11.7/16.4	96/12.0/17.1	97/12.3/17.5
walkbridge+50%	51/-0.1/-1.7	89/8.1/11.4	71/4.0/5.4	92/9.3/14.0	88/7.8/11.8	92/9.3/13.9	92/9.2/13.8	93/9.5/14.3
walkbridge+70%	32/-2.0/-4.6	82/6.1/8.7	49/1.4/2.7	87/7.3/11.5	69/4.4/6.9	87/7.3/11.5	85/6.9/11.0	87/7.4/11.6
walkbridge+90%	15/-3.2/-6.2	67/3.7/5.1	26/0.2/0.6	73/4.8/7.8	36/0.9/1.6	73/4.8/7.7	56/3.3/5.8	74/4.8/7.8
pepper+10%	68/4.9/9.6	99/14.8/20.1	99/15.0/21.8	100/20.5/24.9	74/5.4/11.4	100/20.5/24.9	100/23.2/30.5	100/23.9/31.0
pepper+30%	52/3.1/4.8	98/14.6/18.3	95/10.8/13.6	99/16.8/22.9	73/5.4/11.2	99/16.8/22.9	99/17.7/24.8	100/18.5/25.6
pepper+50%	38/0.3/-1.1	97/12.9/16.1	84/6.1/7.0	99/14.9/21.5	71/5.2/10.6	99/14.8/21.5	99/14.5/21.1	99/15.4/22.4
pepper+70%	25/-1.5/-3.9	95/10.6/13.3	57/2.1/3.4	98/12.5/18.5	61/3.9/7.4	98/12.5/18.5	96/11.4/16.9	98/12.7/18.7
pepper+90%	14/-2.7/-5.5	89/7.2/8.5	27/0.4/0.6	93/8.8/12.7	32/1.2/1.9	93/8.8/12.5	75/4.8/7.9	93/9.0/12.9
mandrill+10%	77/2.7/4.9	93/9.8/11.3	90/4.5/6.9	97/13.1/14.3	87/4.2/9.2	97/13.1/14.3	98/14.4/17.1	98/14.5/17.2
mandrill+30%	61/1.5/2.3	90/7.8/9.0	75/4.0/5.9	92/8.9/10.7	79/3.6/7.2	92/8.9/10.7	93/9.3/11.8	93/9.4/11.9
mandrill+50%	44/-0.9/-2.8	84/5.7/6.6	67/2.7/3.3	87/6.6/8.5	68/2.8/5.2	87/6.6/8.5	87/6.7/8.8	88/6.8/8.8
mandrill+70%	27/-2.7/-5.6	76/3.8/4.3	48/1.1/1.9	80/4.9/6.5	54/2.0/3.6	80/4.9/6.5	79/4.8/6.6	80/4.9/6.5
mandrill+90%	10/-3.8/-7.2	63/2.0/1.9	36/0.3/0.6	69/3.1/4.3	35/0.4/0.8	69/3.1/4.3	59/2.4/3.8	69/3.1/4.4
lake+10%	91/6.6/11.9	99/16.4/22.9	96/11.3/17.6	99/19.6/25.9	99/9.0/17.2	99/19.6/25.7	100/20.3/27.5	100/20.6/27.9
lake+30%	71/3.9/5.6	97/13.6/18.7	90/9.1/12.8	98/15.0/21.4	97/8.6/16.0	98/15.0/21.3	98/15.1/21.7	99/15.4/22.3
lake+50%	52/1.2/-0.4	94/11.2/15.3	76/5.7/6.8	97/12.5/18.3	91/7.7/13.6	97/12.5/18.2	96/12.2/17.9	97/12.7/18.6
lake+70%	33/-0.5/-3.0	90/9.0/12.1	52/2.4/3.7	93/10.4/15.2	63/5.0/8.2	93/10.4/15.2	91/9.7/14.4	94/10.4/15.2
lake+90%	18/-1.6/-4.5	80/6.2/7.5	26/0.5/0.9	84/7.3/10.1	25/1.1/1.9	83/7.3/10.1	51/4.3/7.3	84/7.4/10.2
jetplane+10%	49/2.5/6.0	100/17.0/23.4	98/11.6/17.3	100/20.4/26.8	39/3.4/8.5	100/20.4/26.8	100/20.7/28.0	100/21.3/29.2
jetplane+30%	39/0.6/1.2	98/13.6/17.9	93/8.3/10.4	99/15.5/21.9	40/3.4/8.3	99/15.5/21.9	99/15.3/21.6	99/15.9/22.7
jetplane+50%	33/-1.4/-4.1	96/10.9/14.1	79/4.0/5.1	98/12.7/18.4	39/3.1/7.2	98/12.7/18.4	98/12.1/17.3	98/12.9/18.5
jetplane+70%	30/-2.8/-6.4	93/8.5/10.5	53/0.3/1.2	96/10.2/14.6	32/1.2/3.0	96/10.2/14.6	94/9.2/13.3	96/10.3/14.6
jetplane+90%	28/-3.7/-7.9	87/5.6/6.0	26/-1.7/-2.1	89/6.6/8.6	29/-1.9/-2.8	89/6.6/8.6	54/2.4/4.8	89/6.8/8.7
Mixed Impulse Noise (Half Random-Value Noise and Half Salt-and-Pepper Noise)
walkbridge+10%	91/6.1/10.1	93/10.6/14.7	91/7.5/12.3	96/12.6/13.3	96/12.5/16.0	96/12.6/13.3	98/14.8/17.8	98/15.1/17.9
walkbridge+30%	73/3.6/6.7	90/8.4/11.8	83/6.3/10.3	88/6.6/8.3	89/8.6/12.2	92/8.6/12.2	93/10.2/13.5	93/10.2/12.9
walkbridge+50%	55/1.5/1.9	81/5.7/7.0	70/4.3/6.8	76/3.5/5.7	78/5.7/8.7	85/6.3/10.0	86/7.6/10.8	87/7.6/10.1
walkbridge+70%	37/-0.5/-1.8	63/2.4/1.9	50/2.0/2.9	58/1.9/3.2	56/2.8/4.9	72/4.4/7.2	74/5.1/7.9	80/5.7/7.9
walkbridge+90%	21/-1.9/-4.0	34/-0.6/-2.1	30/0.1/0.4	34/0.3/0.5	31/0.6/1.3	38/1.2/2.0	40/1.3/2.3	63/3.3/4.9
pepper+10%	68/5.0/9.7	98/13.9/19.5	99/15.0/22.0	98/14.3/16.0	74/5.4/11.3	99/14.4/19.9	100/21.0/25.6	99/19.9/23.4
pepper+30%	54/3.7/6.8	97/12.7/16.0	96/11.4/15.4	91/7.5/10.8	72/5.3/10.8	98/12.8/18.5	99/15.8/20.7	98/14.9/18.4
pepper+50%	41/1.8/2.3	92/8.5/8.6	86/7.0/8.9	80/4.5/7.0	68/4.8/9.5	97/11.2/16.1	97/12.6/17.1	97/12.6/15.7
pepper+70%	29/-0.1/-1.2	73/3.6/2.4	62/3.0/3.6	63/2.5/3.8	54/3.3/5.9	90/8.1/10.7	92/9.1/12.5	94/10.1/12.8
pepper+90%	19/-1.4/-3.4	39/-0.2/-2.0	33/0.4/0.5	37/0.6/0.7	31/1.0/1.5	53/2.1/2.5	49/2.2/2.9	82/5.6/6.6
mandrill+10%	76/3.0/5.3	86/6.8/8.3	91/5.5/6.8	95/10.4/10.1	83/3.6/7.3	95/10.5/10.3	96/12.1/12.4	96/11.7/11.2
mandrill+30%	63/1.8/3.4	82/5.4/6.6	74/3.9/6.0	85/5.3/5.1	73/2.9/5.3	88/6.5/7.4	89/7.3/8.1	89/7.3/7.5
mandrill+50%	47/0.6/0.6	75/3.7/4.0	67/3.0/4.4	74/2.5/3.8	61/2.2/3.9	78/4.4/5.6	80/5.0/5.9	81/5.0/5.3
mandrill+70%	32/-1.0/-2.6	60/1.3/0.2	53/1.5/1.8	58/1.3/2.1	48/1.4/2.6	68/2.9/4.2	69/3.3/4.3	73/3.5/3.9
mandrill+90%	20/-2.4/-4.9	35/-1.2/-3.3	36/0.3/0.5	37/0.2/0.1	33/0.3/0.6	46/1.1/1.8	45/1.0/1.3	62/2.2/2.8
lake+10%	91/6.8/12.0	98/14.6/20.5	96/11.3/17.7	98/15.0/15.5	97/8.7/16.1	98/15.0/19.4	99/18.0/22.2	99/17.9/21.2
lake+30%	73/4.3/7.6	95/11.7/15.7	91/9.3/13.7	90/8.0/10.5	92/7.9/13.8	96/11.0/16.4	96/13.1/17.2	96/12.8/15.6
lake+50%	55/2.3/2.7	87/7.9/9.0	75/6.1/8.6	78/4.8/7.2	82/6.6/11.0	92/9.3/13.1	92/10.4/14.1	92/10.0/12.2
lake+70%	37/0.6/-0.6	66/3.7/3.1	44/2.8/3.7	58/2.6/4.1	48/3.7/6.2	82/7.0/9.8	83/7.7/10.8	87/7.9/9.4
lake+90%	22/-0.6/-2.7	34/0.4/-1.1	20/0.6/0.7	30/0.8/1.0	24/0.8/1.5	22/1.5/2.4	33/2.0/3.1	74/5.3/6.0
jetplane+10%	44/2.6/6.0	99/15.4/20.8	98/11.6/17.5	99/13.9/13.3	39/3.4/8.3	99/13.9/19.3	99/17.6/20.8	99/16.8/18.5
jetplane+30%	36/0.8/2.5	97/11.6/14.2	94/8.8/12.3	91/6.3/8.2	38/3.2/7.7	98/11.0/16.6	98/13.1/16.4	98/12.6/14.1
jetplane+50%	30/-0.8/-1.6	90/6.6/6.2	79/4.5/5.8	79/2.8/4.2	37/2.8/6.1	95/9.0/12.7	95/10.0/13.0	95/9.7/10.7
jetplane+70%	25/-2.1/-4.5	68/1.7/-0.1	45/0.6/0.5	60/0.4/0.9	25/0.7/1.9	88/6.3/7.9	87/6.6/8.8	91/7.3/8.0
jetplane+90%	22/-3.1/-6.4	34/-1.8/-4.4	19/-1.8/-2.4	30/-1.5/-2.3	16/-2.1/-3.0	19/-0.8/-0.9	32/-0.2/-0.1	79/4.2/4.4

Table 4. Table IV: General deblurring problems. The results separated by ‘/’ are S ​ N ​ R 0 𝑆 𝑁 subscript 𝑅 0 SNR_{0} , S ​ N ​ R 1 𝑆 𝑁 subscript 𝑅 1 SNR_{1} and S ​ N ​ R 2 𝑆 𝑁 subscript 𝑅 2 SNR_{2} , respectively. The 1 s ​ t superscript 1 𝑠 𝑡 1^{st} , 2 n ​ d superscript 2 𝑛 𝑑 2^{nd} , and 3 r ​ d superscript 3 𝑟 𝑑 3^{rd} best results are colored with red , blue and green , respectively.

	Corrupted	${\ell }_{1}TV$-$SBM$	$TSM$	${\ell }_{p}TV$-$PADMM$	${\ell }_{02}TV$-$AOP$	${\ell }_{0}TV$-$PDA$	${\ell }_{0}TV$-$PADM$
Random-Valued Impulse Noise
walkbridge+10%	63/2.9/3.4	74/4.8/8.6	72/4.6/8.2	77/5.1/9.2	81/5.6/10.1	76/5.0/9.0	91/7.0/13.2
walkbridge+30%	52/1.1/0.0	72/4.6/8.1	61/3.7/6.8	75/4.9/8.7	79/5.4/9.7	74/4.8/8.7	86/6.4/11.7
walkbridge+50%	42/-0.2/-1.9	63/3.8/6.9	46/2.4/4.6	71/4.5/8.0	75/4.9/8.6	73/4.7/8.3	84/6.0/11.0
walkbridge+70%	31/-1.2/-3.2	46/2.1/3.8	33/1.1/2.3	55/2.9/5.1	65/3.3/4.8	69/4.3/7.7	81/5.6/10.1
walkbridge+90%	21/-2.0/-4.2	28/0.3/0.8	25/0.2/0.5	31/0.6/1.2	33/0.4/0.6	42/1.7/3.0	67/3.7/5.8
pepper+10%	81/4.9/4.5	94/9.3/14.7	93/8.3/13.6	70/5.1/10.1	96/9.7/15.8	94/9.0/14.7	99/11.1/19.8
pepper+30%	66/2.1/0.3	92/8.5/13.3	82/5.7/9.9	68/4.9/9.7	96/9.7/15.8	93/8.8/14.1	98/10.7/18.8
pepper+50%	52/0.4/-1.8	83/6.4/9.9	58/3.4/6.0	65/4.6/8.9	95/9.3/14.9	92/8.5/13.5	98/10.4/17.8
pepper+70%	37/-0.8/-3.2	58/3.1/4.7	37/1.6/2.9	52/3.0/5.4	82/5.1/5.4	90/7.8/12.1	97/9.8/16.4
pepper+90%	23/-1.8/-4.3	29/0.6/1.0	24/0.4/0.7	29/0.9/1.3	38/0.9/0.7	54/2.5/3.5	85/6.1/7.2
mandrill+10%	59/1.6/1.3	67/2.9/4.7	65/2.7/4.3	54/2.1/3.8	68/3.0/4.5	68/3.1/5.0	78/4.3/7.3
mandrill+30%	50/0.0/-1.7	66/2.9/4.6	60/2.3/3.9	52/2.1/3.7	68/3.0/4.6	67/3.0/4.8	76/4.0/6.8
mandrill+50%	40/-1.1/-3.4	64/2.7/4.3	50/1.6/2.9	51/2.0/3.5	68/2.9/4.5	66/2.9/4.6	73/3.6/6.0
mandrill+70%	30/-2.0/-4.7	53/1.8/3.1	40/0.9/1.7	46/1.6/2.9	64/2.5/3.6	65/2.7/4.4	70/3.3/5.4
mandrill+90%	21/-2.7/-5.6	38/0.5/0.9	36/0.3/0.6	34/0.4/0.7	42/0.6/0.8	49/1.5/2.5	65/2.7/4.2
lake+10%	71/4.8/4.9	84/7.6/11.6	83/7.3/11.3	83/6.7/11.3	89/8.6/13.8	84/7.7/12.1	96/10.0/17.4
lake+30%	59/2.6/1.2	81/7.1/10.8	65/5.2/8.9	80/6.4/10.7	89/8.5/13.2	83/7.4/11.6	94/9.5/15.9
lake+50%	46/1.1/-0.7	68/5.5/8.8	35/3.2/5.6	76/6.0/9.8	86/7.9/11.9	82/7.2/11.1	92/9.1/15.1
lake+70%	34/0.0/-2.1	35/2.6/4.5	22/1.6/2.9	39/3.3/5.6	66/4.3/5.4	79/6.7/10.2	89/8.5/13.8
lake+90%	22/-0.9/-3.1	22/0.6/1.0	16/0.4/0.8	22/0.7/1.3	21/0.6/0.8	31/2.1/3.5	74/5.6/7.2
jetplane+10%	76/3.3/2.1	88/6.7/9.9	88/6.1/9.7	63/2.8/6.5	93/7.9/12.5	89/6.8/10.5	98/9.1/16.6
jetplane+30%	63/0.7/-1.9	86/6.2/9.1	68/3.2/6.3	66/2.7/6.2	93/7.8/12.0	88/6.6/10.0	97/8.8/15.6
jetplane+50%	49/-0.9/-3.9	74/3.9/6.6	34/0.9/2.6	55/2.5/5.6	91/7.0/9.7	87/6.3/9.4	95/8.4/14.2
jetplane+70%	36/-2.1/-5.3	37/0.3/1.3	22/-0.7/-0.3	35/-0.1/0.6	64/1.5/1.9	84/5.8/8.5	93/7.8/12.4
jetplane+90%	23/-3.0/-6.3	23/-1.7/-2.3	14/-1.9/-2.5	16/-2.2/-3.3	20/-1.7/-2.5	30/0.0/0.6	80/4.5/5.1
Salt-and-Pepper Impulse Noise
walkbridge+10%	61/2.0/0.8	73/4.8/8.5	80/5.6/10.1	76/5.1/9.1	80/5.6/10.1	76/5.0/9.0	94/7.4/14.3
walkbridge+30%	48/-0.5/-3.2	71/4.5/7.9	79/5.4/9.7	74/4.8/8.5	79/5.4/9.7	75/4.9/8.8	92/7.2/13.7
walkbridge+50%	35/-2.1/-5.3	67/4.1/7.3	77/5.2/9.3	72/4.5/8.1	77/5.2/9.3	73/4.8/8.5	90/6.8/12.9
walkbridge+70%	22/-3.3/-6.7	53/2.8/5.2	75/5.0/8.8	61/3.5/6.4	75/4.9/8.8	71/4.5/8.1	86/6.4/11.8
walkbridge+90%	8/-4.2/-7.7	31/0.6/1.0	73/4.7/8.3	34/0.9/1.7	73/4.7/8.3	59/3.4/6.3	79/5.4/9.9
pepper+10%	79/3.6/1.3	94/8.9/14.2	96/9.7/15.8	69/5.0/10.0	96/9.6/15.8	94/9.1/14.8	99/11.4/20.3
pepper+30%	62/0.2/-3.2	92/8.5/13.2	96/9.6/15.7	69/4.9/9.6	96/9.6/15.7	94/8.9/14.4	99/11.2/19.7
pepper+50%	45/-1.7/-5.4	87/7.3/11.2	95/9.4/15.4	66/4.7/9.1	95/9.4/15.4	93/8.6/13.8	99/10.9/19.1
pepper+70%	28/-3.0/-6.8	70/4.3/6.5	95/9.2/14.8	56/3.7/6.8	95/9.2/14.9	91/8.3/13.0	98/10.3/18.2
pepper+90%	11/-4.1/-7.9	33/0.8/1.1	94/8.8/14.1	32/1.1/1.8	94/8.8/14.1	79/5.6/8.8	96/9.5/15.8
mandrill+10%	58/0.7/-1.3	67/2.9/4.7	67/2.9/4.4	53/2.1/3.8	67/2.9/4.4	68/3.1/5.0	86/5.2/9.5
mandrill+30%	45/-1.7/-5.2	65/2.8/4.4	67/2.9/4.5	52/2.1/3.6	67/2.9/4.5	68/3.0/4.9	83/4.9/8.7
mandrill+50%	32/-3.2/-7.2	64/2.6/4.2	66/2.8/4.4	51/2.0/3.5	66/2.8/4.4	67/3.0/4.7	80/4.5/7.9
mandrill+70%	19/-4.4/-8.6	56/2.0/3.3	65/2.7/4.2	48/1.8/3.1	65/2.7/4.2	66/2.8/4.5	75/4.0/6.7
mandrill+90%	7/-5.2/-9.6	39/0.5/0.8	65/2.7/4.2	35/0.5/1.0	65/2.7/4.2	60/2.4/3.9	70/3.3/5.3
lake+10%	69/3.9/2.4	83/7.4/11.4	90/8.7/13.8	82/6.6/11.2	90/8.7/13.8	85/7.7/12.1	98/10.3/18.5
lake+30%	54/1.0/-1.8	81/7.1/10.6	89/8.5/13.4	80/6.3/10.6	89/8.5/13.4	84/7.6/11.8	97/10.1/17.9
lake+50%	38/-0.7/-3.9	76/6.4/9.6	87/8.2/12.9	77/6.0/9.8	87/8.2/12.8	82/7.3/11.3	96/9.8/17.0
lake+70%	23/-1.9/-5.3	49/3.9/6.3	86/7.9/12.2	56/4.4/7.3	86/7.9/12.2	81/7.0/10.7	94/9.3/15.9
lake+90%	8/-2.8/-6.4	24/0.9/1.4	83/7.4/11.2	21/1.0/1.8	84/7.5/11.1	63/5.0/8.1	88/8.2/13.3
jetplane+10%	75/2.3/-0.4	88/6.5/9.7	93/8.0/12.6	67/2.8/6.5	93/8.0/12.6	89/6.9/10.6	99/9.5/17.8
jetplane+30%	58/-0.9/-4.8	86/6.2/9.0	93/7.7/11.9	64/2.7/6.1	92/7.6/11.8	88/6.7/10.2	99/9.4/17.2
jetplane+50%	42/-2.7/-7.0	82/5.4/7.8	91/7.5/11.4	54/2.5/5.7	91/7.5/11.5	87/6.5/9.7	98/9.0/16.2
jetplane+70%	25/-3.9/-8.4	48/1.9/3.8	90/7.1/10.7	39/1.2/2.9	90/7.1/10.6	86/6.1/9.0	96/8.7/14.9
jetplane+90%	8/-4.9/-9.5	24/-1.3/-1.8	89/6.7/9.7	21/-1.9/-2.8	89/6.7/9.9	72/3.6/6.1	92/7.2/11.8
Mixed Impulse Noise (Half Random-Value Noise and Half Salt-and-Pepper Noise)
walkbridge+10%	62/2.4/1.9	74/4.8/8.5	72/4.6/8.2	77/5.1/9.2	81/5.6/10.1	76/5.0/9.0	93/7.4/14.0
walkbridge+30%	50/0.2/-1.9	71/4.5/7.9	65/3.9/7.2	74/4.8/8.6	79/5.4/9.6	74/4.8/8.7	87/6.5/12.0
walkbridge+50%	38/-1.3/-3.9	64/3.8/6.9	52/2.9/5.5	71/4.5/8.0	78/5.2/8.8	73/4.7/8.3	84/6.1/11.0
walkbridge+70%	27/-2.3/-5.3	48/2.4/4.4	38/1.6/3.2	59/3.3/6.0	74/4.5/7.1	70/4.4/7.8	81/5.6/10.1
walkbridge+90%	15/-3.2/-6.3	29/0.5/1.0	27/0.4/0.9	33/0.8/1.6	43/1.1/1.4	50/2.3/3.7	71/4.3/7.2
pepper+10%	80/4.2/2.6	94/9.1/14.5	93/8.5/13.7	69/5.1/10.0	96/9.7/15.9	94/9.0/14.7	99/11.1/19.8
pepper+30%	64/1.0/-1.8	91/8.4/13.0	87/6.4/10.9	68/4.9/9.6	96/9.7/15.8	93/8.8/14.1	99/10.9/19.3
pepper+50%	49/-0.8/-3.9	84/6.7/10.2	68/4.2/7.5	66/4.7/9.1	96/9.4/15.0	92/8.5/13.5	98/10.5/18.2
pepper+70%	33/-2.1/-5.4	61/3.5/5.2	43/2.3/4.0	54/3.4/6.3	94/8.4/11.4	90/7.9/12.2	97/10.0/16.8
pepper+90%	17/-3.1/-6.4	31/0.9/1.3	27/0.7/1.2	32/1.1/1.8	55/2.0/1.5	60/3.2/4.9	92/8.2/11.7
mandrill+10%	58/1.1/-0.2	67/2.9/4.7	65/2.7/4.3	53/2.1/3.8	67/2.9/4.6	68/3.1/5.0	85/5.0/9.4
mandrill+30%	47/-0.9/-3.7	66/2.8/4.5	62/2.4/4.0	52/2.1/3.7	68/3.0/4.6	67/3.0/4.8	76/4.0/6.8
mandrill+50%	36/-2.3/-5.7	64/2.6/4.2	54/1.9/3.3	51/2.0/3.4	68/2.9/4.6	66/2.9/4.6	74/3.7/6.3
mandrill+70%	25/-3.3/-7.0	54/1.9/3.2	43/1.1/2.2	47/1.7/3.1	67/2.8/4.2	65/2.7/4.4	71/3.4/5.4
mandrill+90%	14/-4.2/-8.1	38/0.4/0.7	36/0.4/0.8	35/0.5/0.9	50/1.3/1.4	48/1.0/1.1	66/2.8/4.3
lake+10%	70/4.3/3.5	83/7.5/11.5	83/7.4/11.4	82/6.6/11.3	89/8.6/13.8	84/7.7/12.1	97/10.0/17.9
lake+30%	56/1.7/-0.5	80/7.0/10.6	74/5.8/9.6	80/6.3/10.6	88/8.4/13.3	83/7.5/11.6	94/9.5/16.2
lake+50%	42/0.1/-2.6	73/6.0/9.3	45/4.0/7.0	77/6.0/9.8	88/8.1/11.8	82/7.2/11.1	92/9.1/15.1
lake+70%	29/-1.0/-4.0	40/2.9/5.1	27/2.3/4.0	51/4.1/6.8	84/7.4/10.8	79/6.8/10.3	89/8.5/13.5
lake+90%	15/-2.0/-5.0	18/0.7/1.2	17/0.7/1.3	18/0.9/1.6	32/1.4/1.5	55/3.8/5.4	81/6.8/9.8
jetplane+10%	76/2.8/0.6	88/6.7/9.9	89/6.4/9.8	66/2.8/6.5	93/7.9/12.5	89/6.8/10.5	98/9.1/16.6
jetplane+30%	60/-0.2/-3.6	86/6.2/8.9	79/4.1/7.5	66/2.7/6.1	93/7.8/11.8	88/6.6/9.9	97/8.8/15.6
jetplane+50%	45/-1.9/-5.7	81/5.0/7.5	44/1.9/4.2	51/2.5/5.6	91/7.1/10.3	87/6.4/9.5	95/8.4/14.1
jetplane+70%	30/-3.1/-7.1	39/0.7/2.2	25/0.0/1.0	32/0.8/2.2	89/6.4/8.4	85/5.9/8.7	93/7.7/12.2
jetplane+90%	15/-4.1/-8.2	16/-1.6/-2.1	16/-1.6/-2.0	22/-2.0/-3.0	30/-1.1/-1.8	56/2.0/3.1	86/5.8/7.8

Table 5. Table V: CPU time (in seconds) comparisons. First row: image denoising; second row: image deblurring.

${\ell }_{1}TV$-$SBM$	$TSM$	${\ell }_{p}TV$-ADMM	${\ell }_{02}TV$-AOP	${\ell }_{0}TV$-PDA	${\ell }_{0}TV$-$PADMM$
5$\pm$4	6$\pm$4	15$\pm$4	30$\pm$5	17$\pm$3	14$\pm$4
15$\pm 8$	16$\pm 7$	38$\pm$8	62$\pm$4	39$\pm$7	35$\pm$8

Equations

$$\mathbf{b}=\left((\mathbf{Ku})\odot {\epsilon }_m\right)+{\epsilon }_a,$$

$${\min }_{\mathbf{u}}\text{leavevmode}\ell (\mathbf{Ku},\mathbf{b})+\lambda \text{leavevmode}\Omega ({\nabla }_x\mathbf{u},{\nabla }_y\mathbf{u}),$$

$$\mathcal{T}({\mathbf{u}}_i)=\{\begin{array}{ll}{d}_i,& \text{with probability}\text{leavevmode}{r}_{rv};\\ (\mathbf{Ku}{)}_i,& \text{with probability}\text{leavevmode}1-r_{rv}.\end{array}$$

$$\mathcal{T}({\mathbf{u}}_i)=\{\begin{array}{ll}{u}_{\min },& \text{with probability}\text{leavevmode}{r}_{sp}/2;\\ u_{\max },& \text{with probability}\text{leavevmode}{r}_{sp}/2;\\ (\mathbf{Ku}{)}_i,& \text{with probability}\text{leavevmode}1-r_{sp}.\end{array}$$

$$p(\mathbf{u}|\mathcal{T}(\mathbf{u}))=p(\mathbf{u})\cdot p(\mathcal{T}(\mathbf{u})|\mathbf{u})\text{leavevmode}/\text{leavevmode}p(\mathcal{T}(\mathbf{u})).$$

$${\max }_{\mathbf{u}}\text{leavevmode}\log p(\mathcal{T}(\mathbf{u})|\mathbf{u})+\log p(\mathbf{u}).$$

$$p(\mathcal{T}(\mathbf{u})|\mathbf{u})$$

$$\displaystyle\textstyle\min_{\bm{\mathbf{u}}}\leavevmode\nobreak\ \|\mathbf{Ku}-\mathbf{b}\|_{0}+\lambda\sum_{i=1}^{n}\Big{[}|(\bm{\mathbf{\nabla}}_{x}\mathbf{u})_{i}|^{p}+|(\bm{\mathbf{\nabla}}_{y}\mathbf{u})_{i}|^{p}\Big{]}^{1/p},$$

$$\Vert \mathbf{x}{\Vert }_{p,1}\triangleq {\sum }_{i=1}^n(|{\mathbf{x}}_i{|}^p+|{\mathbf{x}}_{n+i}{|}^p{)}^{\frac{1}{p}};\text{leavevmode}\nabla \triangleq \left[\genfrac{}{}{0ex}{}{{\nabla }_x}{{\nabla }_y}\right]\in {\mathbb{R}}^{2n\times n}.$$

$${\min }_{0\le \mathbf{u}\le 1}\text{leavevmode}\Vert \mathbf{o}\odot \left(\mathbf{Ku}-\mathbf{b}\right){\Vert }_0+\lambda \Vert \nabla \mathbf{u}{\Vert }_{p,1},$$

$$\Vert \mathbf{w}{\Vert }_0=\underset{0\le \mathbf{v}\le 1}{\min }\text{leavevmode}⟨1,1-\mathbf{v}⟩,\text{leavevmode}s.t.\text{leavevmode}\mathbf{v}\odot |\mathbf{w}|=0,$$

$$\begin{array}{rl}& {\min }_{0\le \mathbf{u},\mathbf{v}\le 1}\text{leavevmode}⟨1,1-\mathbf{v}⟩+\lambda \Vert \nabla \mathbf{u}{\Vert }_{p,1}\\ & s.t.\mathbf{v}\odot |\mathbf{o}\odot (\mathbf{Ku}-\mathbf{b})|=0.\end{array}$$

$$\underset{0\le \mathbf{u},\mathbf{v}\le 1,\text{leavevmode}\mathbf{x},\text{leavevmode}\mathbf{y}}{\min }\text{leavevmode}⟨1,1-\mathbf{v}⟩+\lambda \Vert \mathbf{x}{\Vert }_{p,1}$$

$$s.t.\nabla \mathbf{u}=\mathbf{x},\text{leavevmode}\mathbf{Ku}-\mathbf{b}=\mathbf{y},\text{leavevmode}\mathbf{v}\odot \mathbf{o}\odot |\mathbf{y}|=0.$$

$$\displaystyle\textstyle\mathcal{L}(\bm{\mathbf{u}},\bm{\mathbf{v}},\bm{\mathbf{x}},\bm{\mathbf{y}},\bm{\mathbf{\xi}},\bm{\mathbf{\zeta}},\bm{\mathbf{\pi}}):=\langle\bm{\mathbf{1}},\bm{\mathbf{1}}-\bm{\mathbf{v}}\rangle+\lambda\|\bm{\mathbf{x}}\|_{p,1}+\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\$$

$$\displaystyle\textstyle\langle\bm{\mathbf{\nabla}}\bm{\mathbf{u}}-\bm{\mathbf{x}},\bm{\mathbf{\xi}}\rangle+\frac{\beta}{2}\|\bm{\mathbf{\nabla}}\bm{\mathbf{u}}-\bm{\mathbf{x}}\|^{2}+\langle\bm{\mathbf{Ku}}-\bm{\mathbf{b}}-\bm{\mathbf{y}},\bm{\mathbf{\zeta}}\rangle+\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\$$

$$\displaystyle\textstyle\frac{\beta}{2}\|\bm{\mathbf{Ku}}-\bm{\mathbf{b}}-\bm{\mathbf{y}}\|^{2}+\langle\bm{\mathbf{v}}\odot\bm{\mathbf{o}}\odot|\bm{\mathbf{y}}|,\bm{\mathbf{\pi}}\rangle+\frac{\beta}{2}\|\bm{\mathbf{v}}\odot\bm{\mathbf{o}}\odot|\bm{\mathbf{y}}|\|^{2},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\$$

$$\displaystyle\begin{split}\left[\begin{array}[]{c}\bm{\mathbf{u}}^{k+1}\\ \bm{\mathbf{v}}^{k+1}\\ \end{array}\right]=\mathop{\arg\min}_{\bm{\mathbf{0}}\leq\bm{\mathbf{u}},\bm{\mathbf{v}}\leq\bm{\mathbf{1}}}\mathcal{L}(\bm{\mathbf{u}},\bm{\mathbf{v}},\bm{\mathbf{x}}^{k},\bm{\mathbf{y}}^{k},\bm{\mathbf{\xi}}^{k},\bm{\mathbf{\zeta}}^{k},\bm{\mathbf{\pi}}^{k})\\ +\tfrac{1}{2}\|\bm{\mathbf{u}}-\bm{\mathbf{u}}^{k}\|^{2}_{\bm{\mathbf{D}}}+\tfrac{1}{2}\|\bm{\mathbf{v}}-\bm{\mathbf{v}}^{k}\|^{2}_{\bm{\mathbf{E}}}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \end{split}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\$$

$$\displaystyle\left[\begin{array}[]{c}\bm{\mathbf{x}}^{k+1}\\ \bm{\mathbf{y}}^{k+1}\\ \end{array}\right]=\mathop{\arg\min}_{\bm{\mathbf{x}},\bm{\mathbf{y}}}\mathcal{L}(\bm{\mathbf{u}}^{k+1},\bm{\mathbf{v}}^{k+1},\bm{\mathbf{x}},\bm{\mathbf{y}},\bm{\mathbf{\xi}}^{k},\bm{\mathbf{\zeta}}^{k},\bm{\mathbf{\pi}}^{k})$$

$${\xi }^{k+1}={\xi }^k+\gamma \beta (\nabla {\mathbf{u}}^k-{\mathbf{x}}^k),$$

$${\zeta }^{k+1}={\zeta }^k+\gamma \beta (\mathbf{K}{\mathbf{u}}^k-\mathbf{b}-{\mathbf{y}}^k),$$

$${\pi }^{k+1}={\pi }^k+\gamma \beta (\mathbf{o}\odot {\mathbf{v}}^k\odot |{\mathbf{y}}^k|).$$

$$\begin{array}{r}{\mathbf{u}}^{k+1}={\mathrm{argmin}}_{0\le \mathbf{u}\le 1}\frac{\beta }{2}\Vert \nabla \mathbf{u}-{\mathbf{x}}^k+{\xi }^k/\beta {\Vert }^2+\\ \frac{\beta }{2}\Vert \mathbf{Ku}-\mathbf{b}-{\mathbf{y}}^k+{\zeta }^k/\beta {\Vert }^2+\frac{1}{2}\Vert \mathbf{u}-{\mathbf{u}}^k{\Vert }_{\mathbf{D}}^2.\end{array}$$

$${\mathbf{u}}^{k+1}={\mathrm{argmin}}_{0\le \mathbf{u}\le 1}\frac{1}{2}\Vert \mathbf{u}-({\mathbf{u}}^k-{\mathbf{g}}^k/L){\Vert }^2$$

$${\mathbf{u}}^{k+1}=\min (1,\max (0,{\mathbf{u}}^k-{\mathbf{g}}^k/L)).$$

$${\mathbf{v}}^{k+1}={\mathrm{argmin}}_{0\le \mathbf{v}\le 1}\text{leavevmode}\frac{1}{2}{\sum }_{i=1}^n{\mathbf{s}}_i^k{\mathbf{v}}_i^2+⟨\mathbf{v},{\mathbf{c}}^k⟩,$$

$${\mathbf{v}}^{k+1}=\min (1,\max (0,-\frac{{\mathbf{c}}^k}{{\mathbf{s}}^k})).$$

$${\mathbf{x}}^{k+1}={\mathrm{argmin}}_{\mathbf{x}\in {\mathbb{R}}^{2n}}\frac{\beta }{2}\Vert \mathbf{x}-{\mathbf{h}}^k{\Vert }^2+\lambda \Vert \mathbf{x}{\Vert }_{p,1},$$

$$\displaystyle\bm{\mathbf{x}}^{k+1}={\rm sign}\big{(}\bm{\mathbf{h}}^{k}\big{)}\odot\max\big{(}|\bm{\mathbf{h}}^{k}|-\lambda/\beta,0\big{)},$$

$$\displaystyle\textstyle\left[\begin{array}[]{c}\bm{\mathbf{x}}_{i}^{k+1}\\ \bm{\mathbf{x}}_{i+n}^{k+1}\\ \end{array}\right]=\big{(}\max(0,1-\frac{\lambda/\beta}{\|(\bm{\mathbf{h}}_{i}^{k};\bm{\mathbf{h}}_{i+n}^{k})\|})\big{)}\left[\begin{array}[]{c}\bm{\mathbf{h}}_{i}^{k}\\ \bm{\mathbf{h}}_{i+n}^{k}\\ \end{array}\right]$$

$${\mathbf{y}}^{k+1}=\arg {\min }_{\mathbf{y}}\text{leavevmode}\frac{\beta }{2}\Vert \mathbf{y}-{\mathbf{q}}^k{\Vert }^2+\frac{\beta }{2}\Vert {\mathbf{w}}^k\odot |\mathbf{y}|+{\pi }^k/\beta {\Vert }^2,$$

$$\displaystyle\textstyle\bm{\mathbf{y}}^{k+1}=\text{sign}(\bm{\mathbf{q}}^{k})\odot\max\big{(}0,\frac{|\bm{\mathbf{q}}^{k}|-\bm{\mathbf{\pi}}^{k}\odot\bm{\mathbf{w}}^{k}/\beta}{1+\bm{\mathbf{v}}^{k}\odot\bm{\mathbf{w}}^{k}}\big{)},$$

$${\min }_{\mathbf{L}}\text{leavevmode}\Vert \mathbf{B}-\mathbf{L}{\Vert }_0+\lambda \text{leavevmode}rank(\mathbf{L}),$$

Full text

1 Introduction

Image restoration is an inverse problem, which aims at estimating the original clean image $\bm{\mathbf{u}}$ from a blurry and/or noisy observation $\bm{\mathbf{b}}$. Mathematically, this problem is formulated as:

[TABLE]

where $\bm{\mathbf{K}}$ is a linear operator, $\bm{\mathbf{\varepsilon}}_{m}$ and $\bm{\mathbf{\varepsilon}}_{a}$ are the noise vectors, and $\odot$ denotes an elementwise product. Let $\bm{\mathbf{1}}$ and $\bm{\mathbf{0}}$ be column vectors of all entries equal to one and zero, respectively. When $\bm{\mathbf{\varepsilon}}_{m}=\bm{\mathbf{1}}$ and $\bm{\mathbf{\varepsilon}}_{a}\neq\bm{\mathbf{0}}$ (or $\bm{\mathbf{\varepsilon}}_{m}\neq\bm{\mathbf{0}}$ and $\bm{\mathbf{\varepsilon}}_{a}=\bm{\mathbf{0}}$), (1) corresponds to the additive (or multiplicative) noise model. For convenience, we adopt the vector representation for images, where a 2D $M\times N$ image is column-wise stacked into a vector $\bm{\mathbf{u}}\in\mathbb{R}^{n\times 1}$ with $n=M\times N$. So, for completeness, we have $\bm{\mathbf{1}},\bm{\mathbf{0}},\bm{\mathbf{b}},\bm{\mathbf{u}},\bm{\mathbf{\varepsilon}}_{a},\bm{\mathbf{\varepsilon}}_{m}\in\mathbb{R}^{n}$, and $\bm{\mathbf{K}}\in\mathbb{R}^{n\times n}$. Before proceeding, we present an image restoration example on the well-known ‘barbara’ image using our proposed method for solving impulse noise removal in Figure 1.

In general image restoration problems, $\bm{\mathbf{K}}$ represents a certain linear operator, e.g. convolution, wavelet transform, etc., and recovering $\bm{\mathbf{u}}$ from $\bm{\mathbf{b}}$ is known as image deconvolution or image deblurring. When $\bm{\mathbf{K}}$ is the identity operator, estimating $\bm{\mathbf{u}}$ from $\bm{\mathbf{b}}$ is referred to as image denoising [50]. The problem of estimating $\bm{\mathbf{u}}$ from $\bm{\mathbf{b}}$ is called a linear inverse problem which, for most scenarios of practical interest, is ill-posed due to the singularity and/or the ill-conditioning of $\bm{\mathbf{K}}$. Therefore, in order to stabilize the recovery of $\bm{\mathbf{u}}$, it is necessary to incorporate prior-enforcing regularization on the solution. Therefore, image restoration can be modelled globally as the following optimization problem:

[TABLE]

where $\ell(\bm{\mathbf{Ku}},\bm{\mathbf{b}})$ measures the data fidelity between $\bm{\mathbf{Ku}}$ and the observation $\bm{\mathbf{b}}$, $\bm{\mathbf{\nabla}}_{x}\in\mathbb{R}^{n\times n}$ and $\bm{\mathbf{\nabla}}_{y}\in\mathbb{R}^{n\times n}$ are two suitable linear transformation matrices such that $\bm{\mathbf{\nabla}}_{x}\bm{\mathbf{u}}\in\mathbb{R}^{n}$ and $\bm{\mathbf{\nabla}}_{y}\bm{\mathbf{u}}\in\mathbb{R}^{n}$ compute the discrete gradients of the image $\bm{\mathbf{u}}$ along the $x$-axis and $y$-axis, respectively111In practice, one does not need to compute and store the matrices $\bm{\mathbf{\nabla}}_{x}$ and $\bm{\mathbf{\nabla}}_{y}$ explicitly. Since the adjoint of the gradient operator ${\bm{\mathbf{\nabla}}}$ is the negative divergence operator $-\bm{\mathbf{\text{div}}}$, i.e., $\langle\bm{\mathbf{r}},\bm{\mathbf{\nabla}}_{x}\bm{\mathbf{u}}\rangle=\langle-\bm{\mathbf{\text{div}}}_{x}\bm{\mathbf{r}},\bm{\mathbf{u}}\rangle,\langle\bm{\mathbf{s}},\bm{\mathbf{\nabla}}_{y}\bm{\mathbf{u}}\rangle=\langle-\bm{\mathbf{\text{div}}}_{y}\bm{\mathbf{s}},\bm{\mathbf{u}}\rangle$ for any $\bm{\mathbf{r}},\bm{\mathbf{s}}\in\mathbb{R}^{n}$, the inner product between vectors can be evaluated efficiently. Fore more details on the computation of $\bm{\mathbf{\nabla}}$ and div operators, please refer to [14, 51, 4]., $\Omega(\bm{\mathbf{\nabla}}_{x}\bm{\mathbf{u}},\bm{\mathbf{\nabla}}_{y}\bm{\mathbf{u}})$ is the regularizer on $\bm{\mathbf{\nabla}}_{x}\bm{\mathbf{u}}$ and $\bm{\mathbf{\nabla}}_{y}\bm{\mathbf{u}}$, and $\lambda$ is a positive parameter used to balance the two terms for minimization. Apart from regularization, other prior information such as bound constraints [5, 70] or hard constraints can be incorporated into the general optimization framework in (2).

1.1 Related Work

This subsection presents a brief review of existing TV methods, from the viewpoint of data fidelity models, regularization models and optimization algorithms.

Data Fidelity Models: The fidelity function $\ell(\cdot,\cdot)$ in (2) usually penalizes the difference between $\bm{\mathbf{Ku}}$ and $\bm{\mathbf{b}}$ by using different norms/divergences. Its form depends on the assumed distribution of the noise model. Some typical noise models and their corresponding fidelity terms are listed in Table I. The classical TV model [47] only considers TV minimization involving the squared $\ell_{2}$-norm fidelity term for recovering images corrupted by additive Gaussian noise. However, this model is far from optimal when the noise is not Gaussian. Other works [60, 23] extend classical TV to use the $\ell_{1}$-norm in the fidelity term. Since the $\ell_{1}$-norm fidelity term coincides with the probability density function of Laplace distribution, it is suitable for image restoration in the presence of Laplace noise. Moreover, additive uniform noise [22, 51], multiplicative Poisson noise [36], and multiplicative Gamma noise [53] have been considered in the literature. Some extensions have been made to deal with mixed Rayleigh impulse noise and mixed Poisson impulse noise in [2]. Recently, a sparse noise model using an $\ell_{02}$-norm for data fidelity has been investigated in [59] to remove impulse and mixed Gaussian impulse noise. In this paper, we consider $\ell_{0}$-norm data fidelity and show that it is particularly suitable for reconstructing images corrupted with additive/multiplicative 222The impulse noise has a discrete nature (corrupted or uncorrupted), thus it can be viewed as additive noise or multiplicative noise. impulse noise.

Regularization Models: Several regularization models have been studied in the literature (see Table II). The Tikhonov-like regularization [1] function $\Omega_{\text{tik}}$ is quadratic and smooth, therefore it is relatively inexpensive to minimize with first-order smooth optimization methods. However, since this method tends to overly smooth images, it often erodes strong edges and texture details. To address this issue, the total variation (TV) regularizer was proposed by Rudin, Osher and Fatemi in [47] for image denoising. Several other variants of TV have been extensively studied. The original TV norm $\Omega_{\text{tv}_{2}}$ in [47] is isotropic, while an anisotropic variation $\Omega_{\text{tv}_{1}}$ is also used. From a numerical point of view, $\Omega_{\text{tv}_{2}}$ and $\Omega_{\text{tv}_{1}}$ cannot be directly minimized since they are not differentiable. A popular method is to use their smooth approximation $\Omega_{\text{stv}}$ and $\Omega_{\text{hub}}$ (see [46] for details). Very recently, the Potts model $\Omega_{\text{pot}}$ [29, 42, 9], which is based on the $\ell_{0}$-norm, has received much attention. It has been shown to be particularly effective for image smoothing [56] and motion deblurring [57].

Optimization Algorithms: The optimization problems involved in TV-based image restoration are usually difficult due to the non-differentiability of the TV norm and the high dimensionality of the image data. In the past several decades, a plethora of approaches have been proposed, which include PDE methods based on the Euler-Lagrange equation [47], the interior-point method [18], the semi-smooth Newton method [45], the second-order cone optimization method [31], the splitting Bregman method [32, 69], the fixed-point iterative method [21], Nesterov’s first-order optimal method [44, 5], and alternating direction methods [50, 20, 53]. Among these methods, some solve the TV problem in its primal form [50], while others consider its dual or primal-dual forms [18, 23]. In this paper, we handle the TV problem with $\ell_{0}$-norm data fidelity using a primal-dual formulation, where the resulting equality constrained optimization is solved using proximal Alternating Direction Method of Multipliers (PADMM). It is worthwhile to note that the Penalty Decomposition Algorithm (PDA) in [39] can also solve our problem, however, it lacks numerical stability. This motivates us to design a new $\ell_{0}$-norm optimization algorithm in this paper.

1.2 Contributions and Organization

The main contributions of this paper are two-fold. (1) $\ell_{0}$-norm data fidelity is proposed to address the TV-based image restoration problem333We are also aware of Ref. [19] where $\ell_{0}$-norm data fidelity is considered. However, their interpretation from the MAP viewpoint is not correct. . Compared with existing models, our model is particularly suitable for image restoration in the presence of impulse noise. (2) To deal with the resulting NP-hard 444The $\ell_{0}$ norm problem is known to be NP-hard [43], since it is equivalent to NP-complete subset selection problems. $\ell_{0}$ norm optimization, we propose a proximal ADMM to solve an equivalent MPEC form of the problem. A preliminary version of this paper appeared in [63].

The rest of the paper is organized as follows. Section 2 presents the motivation and formulation of the problem for impulse noise removal. Section 3 presents the equivalent MPEC problem and our proximal ADMM solution. Section 4 discusses the connection between our method and prior work. Section 5 provides extensive and comparative results in favor of our $\ell_{0}$TV method. Finally, Section 6 concludes the paper.

2 Motivation and Formulations

2.1 Motivation

This work focuses on image restoration in the presence of impulse noise, which is very common in data acquisition and transmission due to faulty sensors or analog-to-digital converter errors, etc. Moreover, scratches in photos and video sequences can be also viewed as a special type of impulse noise. However, removing this kind of noise is not easy, since corrupted pixels are randomly distributed in the image and the intensities at corrupted pixels are usually indistinguishable from those of their neighbors. There are two main types of impulse noise in the literature [23, 35]: random-valued and salt-and-pepper impulse noise. Let $[u_{\min},u_{\max}]$ be the dynamic range of an image, where $u_{\min}=0$ and $u_{\max}=1$ in this paper. We also denote the original and corrupted intensity values at position $i$ as $\bm{\mathbf{u}}_{i}$ and $\mathcal{T}(\bm{\mathbf{u}}_{i})$, respectively.

Random-valued impulse noise: A certain percentage of pixels are altered to take on a uniform random number $d_{i}\in[u_{\min},u_{\max}]$:

[TABLE]

Salt-and-pepper impulse noise: A certain percentage of pixels are altered to be either $u_{\min}$ or $u_{\max}$:

[TABLE]

The above definition means that impulse noise corrupts a portion of pixels in the image while keeping other pixels unaffected. Expectation maximization could be used to find the MAP estimate of $\bm{\mathbf{u}}$ by maximizing the conditional posterior probability $p(\bm{\mathbf{u}}|\mathcal{T}(\bm{\mathbf{u}}))$, the probability that $\mathbf{u}$ occurs when $\mathcal{T}(\mathbf{u})$ is observed. By the Bayes’ theorem, we have that

[TABLE]

Taking the negative logarithm of the above equation, the estimate is a solution of the following minimization problem:

[TABLE]

We now focus on the two terms in (5). (i) The expression $p(\mathcal{T}(\bm{\mathbf{u}})|\bm{\mathbf{u}})$ can be viewed as a fidelity term measuring the discrepancy between the estimate $\bm{\mathbf{u}}$ and the noisy image $\mathcal{T}(\mathbf{u})$. The choice of the likelihood $p(\mathcal{T}(\mathbf{u})|\mathbf{u})$ depends upon the property of noise. From the definition of impulse noise given above, we have that

[TABLE]

where $r$ is the noise density level as defined in (3) and (4) and $\|\cdot\|_{0}$ counts the number of non-zero elements in a vector. (ii) The term $p(\bm{\mathbf{u}})$ in (5) is used to regularize a solution that has a low probability. We use a prior which has the Gibbs form: $p(\bm{\mathbf{u}})=\frac{1}{\vartheta}\exp(-E(\bm{\mathbf{u}}))$ with $E(\bm{\mathbf{u}})=\sigma\cdot\Omega_{\text{tv}}(\bm{\mathbf{\nabla}}_{x}\mathbf{u},\bm{\mathbf{\nabla}}_{y}\mathbf{u})$. Here, $E(\bm{\mathbf{u}})$ is the TV prior energy functional, $\vartheta$ is a normalization factor such that the TV prior is a probability, and $\sigma$ is the free parameter of the Gibbs measure. Replacing $p(\mathcal{T}(\mathbf{u})|\mathbf{u})$ and $p(\mathbf{u})$ into (5) and ignoring a constant, we obtain the following $\ell_{0}TV$ model:

[TABLE]

where $\lambda$ is a positive number related to $n$, $\sigma$ and $r$. The parameter $p$ can be 1 (anisotropic TV) or $2$ (isotropic TV), and $(\bm{\mathbf{\nabla}}_{x}\mathbf{u})_{i}$ and $(\bm{\mathbf{\nabla}}_{y}\mathbf{u})_{i}$ denote the $i$th component of the vectors $\bm{\mathbf{\nabla}}_{x}\bm{\mathbf{u}}$ and $\bm{\mathbf{\nabla}}_{y}\bm{\mathbf{u}}$, respectively. For convenience, we define $\forall\bm{\mathbf{x}}\in\mathbb{R}^{2n}$:

[TABLE]

In order to make use of more prior information, we consider the following box-constrained model:

[TABLE]

where $\bm{\mathbf{o}}\in\{0,1\}^{n}$ is specified by the user. When $\bm{\mathbf{o}}_{i}$ is 0, it indicates the pixel in position $i$ is an outlier, while when $\bm{\mathbf{o}}_{i}$ is 1, it indicates the pixel in position $i$ is a potential outlier. For example, in our experiments, we set $\bm{\mathbf{o}}=\mathbf{1}$ for the random-valued impulse noise and $\mathbf{o}_{i}\scriptsize=\begin{cases}0,&{\bm{\mathbf{b}}_{i}=u_{\min}\leavevmode\nobreak\ \text{or}\leavevmode\nobreak\ u_{\max}}\\ 1,&{\text{otherwise}}\\ \end{cases}$ for the salt-and-pepper impulse noise. In what follows, we focus on optimizing the general formulation in (6).

2.2 Equivalent MPEC Reformulations

In this section, we reformulate the problem in (6) as an equivalent MPEC from a primal-dual viewpoint. First, we provide the variational characterization of the $\ell_{0}$-norm using the following lemma.

Lemma 1**.**

For any given $\bm{\mathbf{w}}\in\mathbb{R}^{n}$, it holds that

[TABLE]

and $\bm{\mathbf{v}}^{*}=\bm{\mathbf{1}}-{\rm sign}(|\bm{\mathbf{w}}|)$ is the unique optimal solution of the problem in (7). Here, the standard signum function sign is applied componentwise, and ${\rm sign}(0)=0$.

Proof.

The total number of zero elements in $\bm{\mathbf{w}}$ can be computed as $n-\|\bm{\mathbf{w}}\|_{0}=\max_{\bm{\mathbf{v}}\in\{0,1\}}\leavevmode\nobreak\ \sum_{i=1}^{n}\bm{\mathbf{v}}_{i},\leavevmode\nobreak\ s.t.\leavevmode\nobreak\ \bm{\mathbf{v}}\in\Phi$, where $\Phi\triangleq\{\bm{\mathbf{v}}\leavevmode\nobreak\ |\leavevmode\nobreak\ \bm{\mathbf{v}}_{i}\cdot|\bm{\mathbf{w}}_{i}|=0,\leavevmode\nobreak\ \forall i\in[n]\}$. Note that when $\bm{\mathbf{w}}_{i}=0$, $\bm{\mathbf{v}}_{i}=1$ will be achieved by maximization, when $\bm{\mathbf{w}}_{i}\neq 0$, $\bm{\mathbf{v}}_{i}=0$ will be enforced by the constraint. Thus, $\bm{\mathbf{v}}^{*}_{i}=1-{\rm sign}(|\bm{\mathbf{w}}_{i}|)$. Since the objective function is linear, maximization is always achieved at the boundaries of the feasible solution space. Thus, the constraint of $\bm{\mathbf{v}}_{i}\in\{0,1\}$ can be relaxed to $0\leq\bm{\mathbf{v}}_{i}\leq 1$, we have: $\|\bm{\mathbf{w}}\|_{0}=n-\max_{\mathbf{0}\leq\bm{\mathbf{v}}\leq\mathbf{1},\leavevmode\nobreak\ \bm{\mathbf{v}}\in\Phi}\leavevmode\nobreak\ \sum_{i=1}^{n}\bm{\mathbf{v}}_{i}=\min_{\mathbf{0}\leq\bm{\mathbf{v}}\leq\mathbf{1},\leavevmode\nobreak\ \bm{\mathbf{v}}\in\Phi}\leavevmode\nobreak\ {\langle\mathbf{1},\mathbf{1}-\bm{\mathbf{v}}\rangle}$.

∎

The result of Lemma 1 implies that the $\ell_{0}$-norm minimization problem in (6) is equivalent to

[TABLE]

If $\bm{\mathbf{u}}^{*}$ is a global optimal solution of (6), then $(\bm{\mathbf{u}}^{*},\bm{\mathbf{1}}-{\rm sign}(|\bm{\mathbf{Ku}}^{*}-\bm{\mathbf{b}}|))$ is globally optimal to (LABEL:eq:main:mpec). Conversely, if $(\bm{\mathbf{u}}^{*},\bm{\mathbf{1}}-{\rm sign}(|\bm{\mathbf{Ku}}^{*}-\bm{\mathbf{b}}|))$ is a global optimal solution of (LABEL:eq:main:mpec), then $\bm{\mathbf{u}}^{*}$ is globally optimal to (6).

Although the MPEC problem in (LABEL:eq:main:mpec) is obtained by increasing the dimension of the original $\ell_{0}$-norm problem in (6), this does not lead to additional local optimal solutions. Moreover, compared with (6), (LABEL:eq:main:mpec) is a non-smooth non-convex minimization problem and its non-convexity is only caused by the complementarity constraint $\bm{\mathbf{v}}\odot|\bm{\mathbf{o}}\odot(\bm{\mathbf{Ku}}-\bm{\mathbf{b}})|=\bm{\mathbf{0}}$.

Such a variational characterization of the $\ell_{0}$-norm is proposed in [25, 34, 27, 7, 6], but it is not used to develop any optimization algorithms for $\ell_{0}$-norm problems. We argue that, from a practical perspective, improved solutions to (6) can be obtained by reformulating the $\ell_{0}$-norm in terms of complementarity constraints [40, 63, 65, 64, 66, 67]. In the following section, we will develop an algorithm to solve (LABEL:eq:main:mpec) based on proximal ADMM and show that such a “lifting” technique can achieve a desirable solution of the original $\ell_{0}$-norm optimization problem.

3 Proposed Optimization Algorithm

This section is devoted to the solution of (LABEL:eq:main:mpec). This problem is rather difficult to solve, because it is neither convex nor smooth. Our solution is based on the proximal ADM method, which iteratively updates the primal and dual variables of the augmented Lagrangian function of (LABEL:eq:main:mpec).

First, we introduce two auxiliary vectors $\bm{\mathbf{x}}\in\mathbb{R}^{2n}$ and $\bm{\mathbf{y}}\in\mathbb{R}^{n}$ to reformulate (LABEL:eq:main:mpec) as:

[TABLE]

Let $\mathcal{L}:\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{2n}\times\mathbb{R}^{n}\times\mathbb{R}^{2n}\times\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}$ be the augmented Lagrangian function of (9).

[TABLE]

where $\bm{\mathbf{\xi}}$, $\bm{\mathbf{\zeta}}$ and $\bm{\mathbf{\pi}}$ are the Lagrange multipliers associated with the constraints ${\bm{\mathbf{\nabla u}}}=\bm{\mathbf{x}}$, $\bm{\mathbf{Ku}}-\bm{\mathbf{b}}=\bm{\mathbf{y}}$ and $\bm{\mathbf{v}}\odot\bm{\mathbf{o}}\odot|\bm{\mathbf{y}}|=0$, respectively, and $\beta>0$ is the penalty parameter. The detailed iteration steps of the proximal ADM for (9) are described in Algorithm 1. In simple terms, ADM updates are performed by optimizing for a set of primal variables at a time, while keeping all other primal and dual variables fixed. The dual variables are updated by gradient ascent on the resulting dual problem.

Next, we focus our attention on the solutions of the subproblems in (10) and (13) arising in Algorithm 1. We will show that the computation required in each iteration of Algorithm 1 is insignificant.

(i) $(\bm{\mathbf{u}},\bm{\mathbf{v}})$-subproblem. Proximal ADM introduces a convex proximal term to the objective. The specific form of $\bm{\mathbf{D}}$ is chosen to expedite the computation of the closed form solution. The introduction of $\mu$ is to guarantee strongly convexity of the subproblems.

$\bm{\mathbf{u}}$-subproblem in (10) reduces to the following minimization problem:

[TABLE]

After an elementary calculation, subproblem (17) can be simplified as

[TABLE]

with $\bm{\mathbf{g}}^{k}=\bm{\mathbf{\nabla}}^{T}\xi^{k}+\bm{\mathbf{K}}^{T}\bm{\mathbf{\zeta}}^{k}+\beta\bm{\mathbf{\nabla}}^{T}(\bm{\mathbf{x}}^{k}-\bm{\mathbf{\nabla}}\bm{\mathbf{u}}^{k})+\beta\bm{\mathbf{K}}^{T}(\bm{\mathbf{b}}+\bm{\mathbf{y}}^{k}-\bm{\mathbf{K}}\bm{\mathbf{u}}^{k})$. Then, the solution $\bm{\mathbf{u}}^{k}$ of (10) has the following closed form expression:

[TABLE]

Here the parameter $L$ depends on the spectral norm of the linear matrices $\bm{\mathbf{\nabla}}$ and $\bm{\mathbf{K}}$. Using the definition of $\bm{\mathbf{\nabla}}$ and the classical finite differences that $\|\bm{\mathbf{\nabla}}_{y}\|\leq 2$ and $\|\bm{\mathbf{\nabla}}_{y}\|\leq 2$ (see [4, 14, 70]), the spectral norm of $\bm{\mathbf{\nabla}}$ can be computed by: $\|\bm{\mathbf{\nabla}}\|=\|\left(\begin{smallmatrix}\bm{\mathbf{\nabla}}_{x}\\ \bm{\mathbf{0}}\\ \end{smallmatrix}\right)+\left(\begin{smallmatrix}\bm{\mathbf{0}}\\ \bm{\mathbf{\nabla}}_{y}\\ \end{smallmatrix}\right)\|\leq\|\left(\begin{smallmatrix}\bm{\mathbf{\nabla}}_{x}\\ \bm{\mathbf{0}}\\ \end{smallmatrix}\right)\|+\|\left(\begin{smallmatrix}\bm{\mathbf{0}}\\ \bm{\mathbf{\nabla}}_{y}\\ \end{smallmatrix}\right)\|=\|\bm{\mathbf{\nabla}}_{x}\|+\|\bm{\mathbf{\nabla}}_{y}\|\leq 4$.

$\bm{\mathbf{v}}$-subproblem in (10) reduces to the following minimization problem:

[TABLE]

where $\bm{\mathbf{c}}^{k}=o\odot\bm{\mathbf{\pi}}^{k}\odot|\bm{\mathbf{y}}^{k}|-\bm{\mathbf{1}}-\mu\bm{\mathbf{v}}^{k}$, $\bm{\mathbf{s}}^{k}=\beta\bm{\mathbf{o}}\odot\bm{\mathbf{y}}^{k}\odot\bm{\mathbf{y}}^{k}+\mu$. Therefore, the solution $\bm{\mathbf{v}}^{k}$ can be computed as:

[TABLE]

(iii) $(\bm{\mathbf{x}},\bm{\mathbf{y}})$-subproblem. Variable $\bm{\mathbf{x}}$ in (13) is updated by solving the following problem:

[TABLE]

where $\bm{\mathbf{h}}^{k}:=\bm{\mathbf{\nabla}}\bm{\mathbf{u}}^{k+1}+\bm{\mathbf{\xi}}^{k}/\beta$. It is not difficult to check that for $p=1$,

[TABLE]

and when $p=2$,

[TABLE]

Variable $\bm{\mathbf{y}}$ in (13) is updated by solving the following problem:

[TABLE]

where $\bm{\mathbf{q}}^{k}=\bm{\mathbf{Ku}}^{k+1}-\bm{\mathbf{b}}+\bm{\mathbf{\zeta}}^{k}/\beta$ and $\bm{\mathbf{w}}^{k}=\bm{\mathbf{o}}\odot\bm{\mathbf{v}}^{k+1}$. A simple computation yields that the solution $\bm{\mathbf{y}}^{k}$ can be computed in closed form as:

[TABLE]

Proximal ADM has excellent convergence in practice. The global convergence of ADM for convex problems was given by He and Yuan in [33, 20] under the variation inequality framework. However, since our optimization problem in (LABEL:eq:main:mpec) is non-convex, the convergence analysis for ADM needs additional conditions. By imposing some mild conditions, Wen et al. [52] managed to show that the sequence generated by ADM converges to a KKT point. Along a similar line, we establish the convergence property of proximal ADM. Specifically, we have the following convergence result.

Theorem 1**.**

Convergence of Algorithm 1. Let $X\triangleq(\bm{\mathbf{u,v,x,y}})$, $Y\triangleq(\bm{\mathbf{\xi,\zeta,\pi}})$ and $\{X^{k},Y^{k}\}_{k=1}^{\infty}$ be the sequence generated by Algorithm 1. Assume that $\{Y^{k}\}_{k=1}^{\infty}$ is bounded and satisfies $\sum_{k=0}^{\infty}\|Y^{k+1}-Y^{k}\|_{F}^{2}<\infty$. Then any accumulation point of sequence satisfies the KKT conditions of (9).

Proof.

Please refer to Appendix A. ∎

Remark 1. The condition $\sum_{k=0}^{\infty}\|Y^{k+1}-Y^{k}\|_{F}^{2}<\infty$ holds when the multiplier does not change in two consecutive iterations. By the boundedness of the penalty parameter $\beta$ and Eqs (14-16), this condition also indicates that the equality constraints in (9) are satisfied. This assumption can be checked by measuring the violation of the equality constraints. Theorem 1 indicates that when the equality constraint holds, PADMM converges to a KKT point. Though not satisfactory, it provides some assurance on the convergence of Algorithm 1.

Remark 2. Two reasons explain the good performance of our method. (i) It targets a solution to the original problem in (6). (ii) It has monotone and self-penalized properties owing to the complimentarity constraints brought on by the MPEC. Our method directly handles the complimentary constraints in (9)$:\bm{\mathbf{v}}\odot\bm{\mathbf{o}}\odot|\bm{\mathbf{y}}|=\bm{\mathbf{0}}$ with $\bm{\mathbf{v}}\geq\bm{\mathbf{0}}$. These constraints are the only sources of non-convexity for the optimization problem and they characterize the optimality of the KKT solution of (6). These special properties of MPEC distinguish it from general nonlinear optimization [65, 66, 64, 67]. We penalize the complimentary error of $\bm{\mathbf{v}}\odot\bm{\mathbf{o}}\odot|\bm{\mathbf{y}}|$ (which is always non-negative) and ensure that the error is decreasing in every iteration.

4 Connection with Existing Work

In this section, we discuss the connection between the proposed method $\ell_{0}TV$-PADM and prior work.

4.1 Sparse Plus Low-Rank Matrix Decomposition

Sparse plus low-rank matrix decomposition [54, 35] is becoming a powerful tool that effectively corrects large errors in structured data in the last decade. It aims at decomposing a given corrupted image $\bm{\mathbf{B}}$ (which is of matrix form) into its sparse component ($\bm{\mathbf{S}}$) and low-rank component ($\bm{\mathbf{L}}$) by solving: $\min_{\bm{\mathbf{B}},\bm{\mathbf{L}}}\leavevmode\nobreak\ \|\bm{\mathbf{S}}\|_{0}+\lambda\leavevmode\nobreak\ rank(\bm{\mathbf{L}}),\leavevmode\nobreak\ s.t.\leavevmode\nobreak\ \bm{\mathbf{B}}=\bm{\mathbf{L}}+\bm{\mathbf{S}}$. Here the sparse component represents the foreground of an image which can be treated as outliers or impulse noise, while the low-rank component corresponds to the background, which is highly correlated. This is equivalent to the following optimization problem:

[TABLE]

which is also based on $\ell_{0}$-norm data fidelity. While they consider the low-rank prior in their objective function, we consider the Total Variation (TV) prior in ours.

4.2 Convex Optimization Method $\ell_{1}TV$

The goal of image restoration in the presence of impulse noise has been pursued by a number of authors (see, e.g., [60, 23]) using $\ell_{1}TV$, which can be formulated as follows:

[TABLE]

It is generally believed that $\ell_{1}TV$ is able to remove the impulse noise properly. This is because $\ell_{1}$-norm provides the tightest convex relaxation for the $\ell_{0}$-norm over the unit ball in the sense of $\ell_{\infty}$-norm. It is shown in [12] that the problem of minimizing $\|\bm{\mathbf{Ku}}-\bm{\mathbf{b}}\|_{1}$ is equivalent to $\|\bm{\mathbf{Ku}}-\bm{\mathbf{b}}\|_{0}$ with high probability under the assumptions that (i) $\bm{\mathbf{Ku}}-\bm{\mathbf{b}}$ is sparse at the optimal solution $\bm{\mathbf{u}}^{*}$ and (ii) $\bm{\mathbf{K}}$ is a random Gaussian matrix and sufficiently “incoherent” (i.e., number of rows in $\bm{\mathbf{K}}$ is greater than its number of columns). However, these two assumptions required in [12] do not necessarily hold true for our $\ell_{0}TV$ optimization problem. Specifically, when the noise level of the impulse noise is high, $\bm{\mathbf{Ku}}-\bm{\mathbf{b}}$ may not be sparse at the optimal solution $\bm{\mathbf{u}}^{*}$. Moreover, the matrix $\bm{\mathbf{K}}$ is a square identity or ill-conditioned matrix. Generally, $\ell_{1}TV$ will only lead to a sub-optimal solution.

4.3 Adaptive Outlier Pursuit Algorithm

Very recently, Yan [59] proposed the following new model for image restoration in the presence of impulse noise and mixed Gaussian impulse noise:

[TABLE]

where $\chi>0$ is the regularization parameter. They further reformulate the problem above into $\min_{\bm{\mathbf{u}},\leavevmode\nobreak\ \bm{\mathbf{v}}}\leavevmode\nobreak\ \|\bm{\mathbf{v}}\odot(\bm{\mathbf{K}}\bm{\mathbf{u}}-\bm{\mathbf{b}})\|_{2}^{2}+\lambda\leavevmode\nobreak\ \|\bm{\mathbf{\nabla}}\bm{\mathbf{u}}\|_{p,1},\leavevmode\nobreak\ s.t.\leavevmode\nobreak\ \bm{\mathbf{0}}\leq\bm{\mathbf{v}}\leq\bm{\mathbf{1}},\leavevmode\nobreak\ \langle\bm{\mathbf{v}},\bm{\mathbf{1}}\rangle\leq n-k$ and then solve this problem using an Adaptive Outlier Pursuit(AOP) algorithm. The AOP algorithm is actually an alternating minimization method, which separates the minimization problem over $\bm{\mathbf{u}}$ and $\bm{\mathbf{v}}$ into two steps. By iteratively restoring the images and updating the set of damaged pixels, it is shown that AOP algorithm outperforms existing state-of-the-art methods for impulse noise denoising, by a large margin.

Despite the merits of the AOP algorithm, we must point out that it incurs three drawbacks, which are unappealing in practice. First, the formulation in (24) is only suitable for mixed Gaussian impulse noise, i.e. it produces a sub-optimal solution when the observed image is corrupted by pure impulse noise. (ii) Secondly, AOP is a multiple-stage algorithm. Since the minimization sub-problem over $\bm{\mathbf{u}}$555It actually reduces to the $\ell_{2}TV$ optimization problem. needs to be solved exactly in each stage, the algorithm may suffer from slow convergence. (iii) As a by-product of (i), AOP inevitably introduces an additional parameter (that specifies the Gaussian noise level), which is not necessarily readily available in practical impulse denoising problems.

In contrast, our proposed $\ell_{0}$TV method is free from these problems. Specifically, (i) as have been analyzed in Section 2, i.e. our $\ell_{0}$-norm model is optimal for impulse noise removal. Thus, our method is expected to produce higher quality image restorations, as seen in our results. (ii) Secondly, we have integrated $\ell_{0}$-norm minimization into a unified proximal ADM optimization framework, it is thus expected to be faster than the multiple stage approach of AOP. (iii) Lastly, while the optimization problem in (24) contains two parameters, our model only contains one single parameter.

4.4 Other $\ell_{0}$-Norm Optimization Techniques

Actually, the optimization technique for the $\ell_{0}$-norm regularization problem is the key to removing impulse noise. However, existing solutions are not appealing. The $\ell_{0}$-norm problem can be reformulated as a $0\text{-}1$ mixed integer programming [8]problem which can be solved by a tailored branch-and-bound algorithm but it involves high computational complexity. The simple projection methods are inapplicable to our model since they assume the objective function is smooth. Similar to the $\ell_{1}$ relaxation, the convex methods such as $k$-support norm relaxation [41], $k$-largest norm relaxation [62], QCQP and SDP relaxations [15] only provide loose approximation of the original problem. The non-convex methods such as Schatten $\ell_{p}$ norm [28, 37], re-weighted $\ell_{1}$ norm [13], $\ell_{\text{1-2}}$ norm DC (difference of convex) approximation [61], the Smoothly Clipped Absolute Deviation (SCAD) penalty method[68], the Minimax Concave Plus (MCP) penalty method [26] only produce sub-optimal results since they give approximate solutions for the $\ell_{0}TV$ problem or incur high computational overhead.

We take $\ell_{p}$ norm approximation method for example and it may suffer two issues. First, it involves an additional hyper-parameter $p$ which may not be appealing in practice. Second, the $\ell_{p}$ regularized norm problem for general $p$ could be difficult to solve. This includes the iterative re-weighted least square method [38] and proximal point method. The former approximates $\|\bm{\mathbf{x}}\|_{p}^{p}$ by $\sum_{i=1}^{n}(\bm{\mathbf{x}}_{i}^{2}+\epsilon)^{p/2}$ with a small parameter $\epsilon$ and solves the resulting re-weighted least squares sub-problem which reduces to a weighted $\ell_{2}TV$ problem. The latter needs to evaluate a relatively expensive proximal operator $\Pi(\bm{\mathbf{a}})=\min_{\bm{\mathbf{x}}}\leavevmode\nobreak\ \frac{1}{2}\|\bm{\mathbf{x}}-\bm{\mathbf{a}}\|_{2}^{2}+\lambda\|\bm{\mathbf{x}}\|_{p}^{p}$ in general, except that it has a closed form solution for some special values such as $p=\frac{1}{2}$ and $p=\frac{2}{3}$ [58].

Recently, Lu et al. propose a Penalty Decomposition Algorithm (PDA) for solving the $\ell_{0}$-norm optimization algorithm [39]. As has been remarked in [39], direct ADM on the $\ell_{0}$ norm problem can also be used for solving $\ell_{0}TV$ minimization simply by replacing the quadratic penalty functions in the PDA by augmented Lagrangian functions. Nevertheless, as observed in our preliminary experiments and theirs, the practical performance of direct ADM is worse than that of PDA.

Actually, in our experiments, we found PDA is unstable. The penalty function can reach very large values $(\geq 10^{8})$, and the solution can be degenerate when the minimization problem of the augmented Lagrangian function in each iteration is not exactly solved. This motivates us to design a new $\ell_{0}$-norm optimization algorithm in this paper. We consider a proximal ADM algorithm to the MPEC formulation of $\ell_{0}$-norm since it has a primal-dual interpretation. Extensive experiments have demonstrated that proximal ADM for solving the “lifting” MPEC formulation for $\ell_{0}TV$ produces better image restoration qualities.

5 Experimental Validation

In this section, we provide empirical validation for our proposed $\ell_{0}TV$-PADMM method by conducting extensive image denoising experiments and performing a thorough comparative analysis with the state-of-the-art.

In our experiments, we use 5 well-known test images of size $512\times 512$. All code is implemented in MATLAB using a 3.20GHz CPU and 8GB RAM. Since past studies [11, 21] have shown that the isotropic TV model performs better than the anisotropic one, we choose $p=2$ as the order of the TV norm here. In our experiments, we apply the following algorithms:

(i) BM3D is an image denoising strategy based on an enhanced sparse representation in transform-domain. The enhancement of the sparsity is achieved by grouping similar 2D image blocks into 3D data arrays [24].

(ii) MFM, Median Filter Methods. We utilize adaptive median filtering to remove salt-and-pepper impulse noise and adaptive center-weighted median filtering to remove random-valued impulse noise.

(iii) $\ell_{1}TV$-SBM, the Split Bregman Method (SBM) of [32], which has been implemented in [30]. We use this convex optimization method as our baseline implementation.

(iv) TSM, the Two Stage Method[16, 17, 10]. The method first detects the damaged pixels by MFM and then solves the TV image inpainting problem.

(v) $\ell_{p}TV$-ADMM (direct). We directly use ADMM (Alternating Direction Method of Multipliers) to solve the non-smooth non-convex $\ell_{p}$ problem with proximal operator being computed analytically. We only consider $p=\frac{1}{2}$ in our experiments [58].

(vi) $\ell_{02}TV$-AOP, the Adaptive Outlier Pursuit (AOP) method described in [59]. We use the implementation provided by the author. Here, we note that AOP iteratively calls the $\ell_{1}TV$-$SBM$ procedure, mentioned above.

(vii) $\ell_{0}TV$-PDA, the Penalty Decomposition Algorithm (PDA) [39] for solving the $\ell_{0}TV$ optimization problem in (6).

(viii) $\ell_{0}TV$-PADMM, the proximal ADMM described in Algorithm 1 for solving the $\ell_{0}TV$ optimization problem in (6). We set the relaxation parameter to 1.618 and the strongly convex parameter $\mu$ to $0.01$. All MATLAB codes to reproduce the experiments of this paper are available online at the authors’ research webpages.

5.1 Experiment Setup

For the denoising and deblurring test, we use the following strategies to generate artificial noisy images.

(a) Denoising problem. We corrupt the original image by injecting random-value, salt-and-pepper noise, and mixed noise (half random-value and half salt-and-pepper) with different densities (10% to 90%) to the images.

(b) Deblurring problem. Although blurring kernel estimation has been pursued by many studies (e.g. [55]), here we assume that the blurring kernel is known beforehand. We blur the original images with a $9\times 9$ Gaussian blurring kernel and add impulse noise with different densities (10% to 90%). We use the following MATLAB scripts to generate a blurring kernel of radius $r$ ($r$ is set to 7 in the experiments):

[TABLE]

We run all the previously mentioned algorithms on the generated noisy and blurry images. For $\ell_{02}TV$-AOP, we adapt the author’s image denoising implementation to the image deblurring setting. Since both BM3D and Median Filter Methods (MFM) are not convenient to solve the deblurring problems, we do not test them in the deblurring problem. We terminate $\ell_{0}TV$-PADMM whenever $\|\bm{\mathbf{\nabla}}\bm{\mathbf{u}}^{k}-\bm{\mathbf{x}}^{k}\|_{2}\leq\frac{1}{255}$ and $\|\bm{\mathbf{K}}\bm{\mathbf{u}}^{k}-\bm{\mathbf{b}}-\bm{\mathbf{y}}^{k}\|_{2}\leq\frac{1}{255}$ and $\|\bm{\mathbf{o}}\odot\bm{\mathbf{v}}^{k}\odot|\bm{\mathbf{y}}^{k}|\|_{2}\leq\frac{1}{255}$. For $\ell_{p}TV$-PADMM, $\ell_{0}TV$-PDA, and $\ell_{0}TV$-PADMM, we use the same stopping criterion to terminate the optimization. For $\ell_{1}TV$-SBM and $\ell_{02}TV$-AOP, we adopt the default stopping conditions provided by the authors. For the regularization parameter $\lambda$, we swept over $\{0.1,0.6,1.1,...,9.6\}$. For the regularization parameter $\chi$ in $\ell_{02}TV$-AOP, we swept over $\{10,50,100,500,1000,5000,10000,50000\}$ and set $k$ to the number of corrupted pixels.

To evaluate these methods, we compute their Signal-to-Noise Ratios (SNRs). Since the corrupted pixels follow a Bernoulli-like distribution, it is generally hard to measure the data fidelity between the original images and the recovered images. Therefore, we consider three ways to measure SNR.

[TABLE]

where $\bm{u}^{0}$ is the original clean image and $\bar{\bm{u}}$ is the mean intensity value of $\bm{u}^{0}$, and $\|\cdot\|_{0\text{-}\epsilon}$ is the soft $\ell_{0}$-norm which counts the number of elements whose magnitude is greater than a threshold $\epsilon$. We adopt $\epsilon=\frac{20}{255}$ in our experiments.

5.2 Convergence of $\ell_{0}TV$-PADMM

Here, we verify the convergence property of our $\ell_{0}TV$-PADMM method on denoising and deblurring problems by considering the ‘cameraman’ image subject to 30% random-valued impulse noise. We set $\lambda=8$ for this problem. We record the objective and SNR values for $\ell_{0}TV$-PADMM at every iteration $k$ and plot these results in Figure 2.

We make two important observations from these results. (i)) The objective value (or the SNR value) does not necessarily decrease (or increase) monotonically, and we attribute this to the non-convexity of the optimization problem and the dynamic updates of the penalty factor in Algorithm 1. (ii) The objective and SNR values stabilize after the $120$th iteration, which means that our algorithm has converged, and the increase of the SNR value is negligible after the $80$th iteration. This implies that one may use a looser stopping criterion without sacrificing much restoration quality.

5.3 General Image Denoising Problems

In this subsection, we compare the performance of all 6 methods on general denoising problems. Table III shows image recovery results when random-value or salt-and-pepper or mixed impulse noise is added. Figure 3 shows image recovery results with varying the regularization parameter $\lambda$. For $\ell_{02}TV$ model in (24), the parameter $\chi$ is scaled to the range $[0,10]$ for better visualization. We make the following interesting observations. (i) The $\ell_{02}TV$-AOP method greatly improves upon $\ell_{1}TV$-SBM, MFM and TSM, by a large margin. These results are consistent with the reported results in [59]. (ii) The $\ell_{0}TV$-PDA method outperforms $\ell_{02}TV$-AOP in most test cases because it adopts the $\ell_{0}$-norm in the data fidelity term. (iii) In the case of random-value impulse noise, our $\ell_{0}TV$-PADMM method is better than $\ell_{0}TV$-PDA in $SNR_{0}$ value while it is comparable to $\ell_{0}TV$-PDA in $SNR_{1}$ and $SNR_{2}$. On the other hand, when salt-and-pepper impulse noise is added, we find that $\ell_{0}TV$-PADMM outperforms $\ell_{0}TV$-PDA in most test cases. Interestingly, the performance gap between $\ell_{0}TV$-PADMM and $\ell_{0}TV$-PDA grows larger, as the noise level increases. (iv) For the same noise level, $\ell_{0}TV$-PADMM achieves better recovery performance in the presence of salt-and-pepper impulse noise than random-valued impulse noise. This is primarily due to the fact that random-valued noise can take any value between 0 and 1, thus, making it more difficult to detect which pixels are corrupted.

5.4 General Image Deblurring Problems

In this subsection, we demonstrate the performance of all methods with their optimal regularization parameters on general deblurring problems. Table IV shows the recovery results for random-valued impulse noise, salt-and-pepper impulse noise, and mixed impulse noise, respectively. Figure 5 shows image recovery results with varying the regularization parameter. We have the following interesting observations. (i) $\ell_{02}TV$-AOP significantly outperforms $\ell_{1}TV$-SBM, and the performance gap becomes larger as the noise level increases. This is because the key assumption in the $\ell_{1}$ model is that $Ku-b$ is sparse at the optimal solution $u^{*}$. This does not hold when the noise level is high. (ii) $\ell_{0}TV$-PDA outperforms $\ell_{02}TV$-AOP for high level ($\geq 30\%$) random-valued impulse noise. However, for salt-and-pepper impulse noise, $\ell_{0}TV$-PDA gives worse performance than $\ell_{02}TV$-AOP in most cases. This phenomenon indicates that the Penalty Decomposition Algorithm is not stable for deblurring problems. (iii) By contrast, our $\ell_{0}TV$-PADMM consistently outperforms all methods, especially when the noise level is large. We attribute this result to the “lifting” technique that is used in our optimization algorithm.

Finally, we also report the performance of all methods with sweeping the radius parameter $r$ as in (25) over $\{1,4,7,...,20\}$ in Figure 5. We notice that the restoration quality degenerates as the radius of the kernel increases for all methods. However, our method consistently gives the best performance.

5.5 Scratched Image Denoising Problems

In this subsection, we demonstrate the superiority of the proposed $\ell_{0}TV$-PADMM in real-world image restoration problems. Specifically, we corrupt the images with scratches which can be viewed as impulse noise666Note that this is different from the classical image inpainting problem that assumes the mask is known. In our scratched image denoising problem, we assume the mask is unknown., see Figure 7. We only consider recovering images using $\ell_{02}TV$-AOP, $\ell_{0}TV$-PDA and $\ell_{0}TV$-PADMM. We show the recovered results in Figure 7. For better visualization of the images recovered by all methods, we also show auxiliary images $\bm{\mathbf{c}}$ in Figure 9, which show the complement of the absolute residual between the recovered image $\bm{\mathbf{u}}$ and the corrupted image $\bm{\mathbf{b}}$ (i.e., $\bm{\mathbf{c}}=\{\bm{\mathbf{1}}-|\bm{\mathbf{b}}-\bm{\mathbf{u}}|\}$). Note that when $\bm{\mathbf{c}}_{i}$ is approximately equal to $1$, the color of the corresponding pixel at position $i$ in the image is white. A conclusion can be drawn that our method $\ell_{0}TV$-PADMM generates more ‘white’ images $\bm{\mathbf{c}}$ than the other two methods, since it can identify the ‘right’ outliers in the corrupted image and make the correction using their neighborhood information.

5.6 Colored Image Denoising Problems

Our proposed method can be directly extended to its color version. Since color total variation is not the main theme of this paper, we only provide a basic implementation of it. Specifically, we compute the color total variation channel-by-channel, and take a $\ell_{1}$-norm of the resulting vectors. Suppose we have RGB channels, then we have the following optimization problem:

[TABLE]

where $\bm{\mathbf{o}}^{k}$ and $\bm{\mathbf{u}}^{k}$ are the prior and the solution of the $k$th channel. The grayscale proximal ADM algorithm in Algorithm 1 can be directly extended to solve the optimization above. We demonstrate its applicability in colored image denoising problems in Figure 9. The regularization parameter $\lambda$ is set to $8$ for the three images in our experiments.

5.7 Running Time Comparisons

We provide some running time comparisons for the methods $\ell_{1}TV$-SBM, TSM, $\ell_{p}TV$-ADMM, $\ell_{02}TV$-AOP, $\ell_{0}TV$-PDA, and $\ell_{0}TV$-PADMM on grayscale image ‘cameraman’ corrupted by 50% random-value impulse noise. For RGB color images, the running time is three times the amount of grayscale images since the colored image recovery problem can be decomposed into dependent subproblems. Table V shows the average CPU time for five runs. Generally, our method is efficient and comparable with existing solutions. This is expected since our method is an alternating optimization algorithm.

6 Conclusions

In this paper, we propose a new method for image restoration based on total variation (TV) with $\ell_{0}$-norm data fidelity, which is particularly suitable for removing impulse noise. Although the resulting optimization model is non-convex, we design an efficient and effective proximal ADM method for solving the equivalent MPEC problem of the original $\ell_{0}$-norm minimization problem. Extensive numerical experiments indicate that the proposed $\ell_{0}$TV model significantly outperforms the state-of-the-art in the presence of impulse noise. In particular, our proposed proximal ADM solver is more effective than the penalty decomposition algorithm used for solving the $\ell_{0}$TV problem [39].

Acknowledgments. We would like to thank Prof. Shaohua Pan for her helpful discussions on this paper. We also thank Prof. Ming Yan for sharing his code with us. This work was supported by the King Abdullah University of Science and Technology (KAUST) Office of Sponsored Research and, in part, by the NSF-China (61772570, 61402182).

Appendix A Proof of Theorem 1

Proof.

We define $Z\triangleq(X,Y)$ and denote $I(\cdot)$ as the indicator function on the constrained set $\Delta\triangleq\{z\leavevmode\nobreak\ |\leavevmode\nobreak\ \bm{\mathbf{0}}\leq\bm{\mathbf{z}}\leq\bm{\mathbf{1}}\}$. First of all, we present the first-order KKT conditions of the MPEC reformulation. Based on the augmented Lagrangian function $\mathcal{L}$, we naturally derive the following KKT conditions for $\{\bm{\mathbf{u}}^{*},\bm{\mathbf{v}}^{*},\bm{\mathbf{x}}^{*},\bm{\mathbf{y}}^{*},\bm{\mathbf{\xi}}^{*},\bm{\mathbf{\zeta}}^{*},\bm{\mathbf{\pi}}^{*}\}$:

[TABLE]

Secondly, we prove that the solution is convergent: $Z^{k+1}-Z^{k}\rightarrow 0$. We observe that $\mathcal{L}$ can be rewritten as:

[TABLE]

Since $Y\triangleq(\bm{\mathbf{\xi,\zeta,\pi}})$ is bounded by assumption, $\mathcal{L}(Z)$ is bounded below for all $Z$. We now define $\mathcal{J}(Z)$ as:

[TABLE]

where $\bm{\mathbf{u}}^{\prime}$ and $\bm{\mathbf{v}}^{\prime}$ denote the values of $\bm{\mathbf{u}}$ and $\bm{\mathbf{v}}$ in the previous iteration. We define $Z^{-1}=Z^{0}$, and the variable $Z$ in $\mathcal{J}(Z)$ is in the range of $\{Z^{0},\leavevmode\nobreak\ Z^{1},\leavevmode\nobreak\ Z^{2},...\}$. Since $\mathcal{J}(Z)$ is strongly and jointly convex with respect to $\{\bm{\mathbf{u}},\bm{\mathbf{v}}\}$ and $\{\bm{\mathbf{u}}^{k+1},\leavevmode\nobreak\ \bm{\mathbf{v}}^{k+1}\}$ is the minimizer of $\min_{\bm{\mathbf{u}},\bm{\mathbf{v}}}\leavevmode\nobreak\ \mathcal{J}(\bm{\mathbf{u}},\leavevmode\nobreak\ \bm{\mathbf{v}},\leavevmode\nobreak\ \bm{\mathbf{x}}^{k},\leavevmode\nobreak\ \bm{\mathbf{y}}^{k},\leavevmode\nobreak\ Y^{k})$ which is based on $\{\bm{\mathbf{u}}^{k},\leavevmode\nobreak\ \bm{\mathbf{v}}^{k}\}$, using the second order growth condition, we have:

[TABLE]

Using the same methodology for the variable $\bm{\mathbf{x}}$ and $\bm{\mathbf{y}}$, we have the following inequalities:

[TABLE]

Denoting $\rho=\tfrac{1}{2}\min(\mu,\beta)$ and combining (27) and (28), we obtain:

[TABLE]

Using the definition of $\mathcal{J}$ and the update rule of the multipliers, we have:

[TABLE]

Combining (29) and (30), we have:

[TABLE]

Taking summation of the above inequality and using the boundedness of $\mathcal{J}(Z)$, we have that:

[TABLE]

Since the second term in the inequality above is bounded, i.e. $\sum_{k=0}^{\infty}\lim_{k\rightarrow\infty}\|Y^{k}-Y^{k+1}\|_{F}^{2}=0$, we obtain that $\sum_{k=0}^{\infty}\lim_{k\rightarrow\infty}\|X^{k}-X^{k+1}\|_{F}^{2}=0$ and $X^{k}-X^{k+1}\rightarrow 0$.

Finally, we are ready to prove the result of the theorem. By the update rule of $Y^{k}$, we have:

[TABLE]

Using the convergence of $Y$ that $Y^{k}-Y^{k+1}\rightarrow 0$ and the optimality of $X^{k+1}$ with respect to $\mathcal{J}(\cdot)$, we have:

[TABLE]

Combining the convergence of $X$ that: $X^{k}-X^{k+1}\rightarrow 0$, we have

[TABLE]

which coincides with the KKT condition in (A). Therefore, $Z^{k+1}$ asymptotically converges to the KKT point.

∎