A Convex Similarity Index for Sparse Recovery of Missing Image Samples

TL;DR

This paper introduces a new convex perceptual similarity index called CSIM for sparse recovery of missing image samples, enabling globally optimal solutions with improved image reconstruction quality.

Contribution

It proposes a convex similarity measure based on SSIM, along with an ADMM-based algorithm for sparse recovery that guarantees convergence to the global optimum.

Findings

CSIM outperforms traditional metrics like MSE in image recovery tasks.

The proposed method guarantees convergence to the global optimum.

Simulation results validate the effectiveness of CSIM and the recovery algorithm.

Abstract

This paper investigates the problem of recovering missing samples using methods based on sparse representation adapted especially for image signals. Instead of $l_2$-norm or Mean Square Error (MSE), a new perceptual quality measure is used as the similarity criterion between the original and the reconstructed images. The proposed criterion called Convex SIMilarity (CSIM) index is a modified version of the Structural SIMilarity (SSIM) index, which despite its predecessor, is convex and uni-modal. We derive mathematical properties for the proposed index and show how to optimally choose the parameters of the proposed criterion, investigating the Restricted Isometry (RIP) and error-sensitivity properties. We also propose an iterative sparse recovery method based on a constrained $l_1$-norm minimization problem, incorporating CSIM as the fidelity criterion. The resulting convex optimization problem is solved via an algorithm based on Alternating Direction Method of Multipliers (ADMM). Taking advantage of the convexity of the CSIM index, we also prove the convergence of the algorithm to the globally optimal solution of the proposed optimization problem, starting from any arbitrary point. Simulation results confirm the performance of the new similarity index as well as the proposed algorithm for missing sample recovery of image patch signals.

Figures

Tables

Table 1. TABLE I: Performance Comparison of the FIR denoising filters of order m 𝑚 m for simulations in part V-A (Input SNR= 1 ​ d ​ B 1 d B 1\mathrm{dB} )

		$m=6$			$m=12$
		MSE Filter	CSIM Filter	SSIM Filter	MSE Filter	CSIM Filter	SSIM Filter
Peppers	PSNR (dB)	21.64793	21.7143	21.5236	21.49739	21.5663	21.547
	SSIM	0.483457	0.487075	0.48695	0.541758	0.547965	0.545625
	FSIM	0.737972	0.743853	0.74152	0.760147	0.762156	0.759996
	Runtime (s)	8.84375	8.671875	123.4844	12.5	12.90625	126.6406
Lena	PSNR (dB)	23.4301	23.51699	23.5393	23.27777	23.29914	23.31125
	SSIM	0.509773	0.51945	0.516033	0.554572	0.559539	0.55748
	FSIM	0.76562	0.770744	0.769893	0.779354	0.780314	0.780307
	Runtime (s)	8.546875	8.875	121.7813	12.57813	12.73438	125.0469
Barbara	PSNR (dB)	22.07703	22.09299	22.06136	21.86944	21.86258	21.80966
	SSIM	0.471681	0.474548	0.47164	0.478401	0.479339	0.473757
	FSIM	0.757897	0.759037	0.755325	0.760586	0.763256	0.756748
	Runtime (s)	8.609375	9.09375	122.1875	12.51563	13.15625	125.1719
House	PSNR (dB)	23.43906	23.5731	23.52321	24.12948	24.16865	24.10161
	SSIM	0.386421	0.394148	0.390097	0.476552	0.48155	0.479555
	FSIM	0.693577	0.695039	0.694965	0.731814	0.733746	0.732695
	Runtime (s)	8.6875	9.09375	124.75	12.64063	13.20313	129.3281
Cameraman	PSNR (dB)	22.06733	22.07957	22.05864	22.0714	22.06688	22.10787
	SSIM	0.382708	0.396937	0.39702	0.453836	0.457089	0.454604
	FSIM	0.674499	0.683098	0.681842	0.687971	0.693188	0.689301
	Runtime (s)	8.96875	8.921875	125.5	12.75	13.14063	123.0781
Couple	PSNR (dB)	22.70483	22.67239	22.69083	22.37438	22.37977	22.38362
	SSIM	0.46712	0.467596	0.470518	0.458069	0.457024	0.454858
	FSIM	0.758213	0.75606	0.758453	0.730827	0.732256	0.728856
	Runtime (s)	8.625	9.25	124.0313	12.29688	13.42188	125.4844

Equations

$$\mathrm{SSIM}({\bf{x}},{\bf{y}})=\Big{(}\dfrac{2\mu_{x}\mu_{y}+C_{1}}{\mu_{x}^{2}+\mu_{y}^{2}+C_{1}}\Big{)}\Big{(}\dfrac{2\sigma_{x,y}+C_{2}}{\sigma_{x}^{2}+\sigma_{y}^{2}+C_{2}}\Big{)}$$

$$\mathrm{CSIM}(\mathbf{x},\mathbf{y})=k_1({\mu }_x^2+{\mu }_y^2-2{\mu }_x{\mu }_y)+k_2({\sigma }_x^2+{\sigma }_y^2-2{\sigma }_{x,y})$$

$$\sigma_{x}^{2}=\frac{1}{n-1}||{\bf{x}}-{\mbox{\boldmath$\mu$}}_{x}||_{2}^{2},\quad\sigma_{x,y}=\frac{1}{n-1}({\bf{x}}-{\mbox{\boldmath$\mu$}}_{x})^{T}({\bf{y}}-{\mbox{\boldmath$\mu$}}_{y})$$

$$\mathrm{CSIM}(\mathbf{x},\mathbf{y})$$

$$\mathrm{CSIM}(\mathbf{e})$$

$$={\mathbf{e}}^T\left(\frac{k_1}{n^2}{1}_n{1}_n^T+\frac{k_2}{n-1}\mathbf{M}\right)\mathbf{e}={\mathbf{e}}^T\mathbf{W}\mathbf{e}$$

$$\mathbf{W}=\frac{k_2}{n-1}{\mathbf{I}}_n+\left(\frac{k_1}{n^2}-\frac{k_2}{n(n-1)}\right)1_n{1}_n^T={\theta }_1{\mathbf{I}}_n+{\theta }_2{1}_n{1}_n^T$$

$$\mathbf{W}^{1/2}=\sqrt{\frac{k_{2}}{n-1}}\mathbf{I}+\frac{1}{n}\Bigg{(}\sqrt{\frac{k_{1}}{n}}-\sqrt{\frac{k_{2}}{n-1}}\Bigg{)}{\bf{1}}_{n}{\bf{1}}_{n}^{T}$$

$$\mathrm{CSIM}(\mathbf{x},\mathbf{y})=(\mathbf{x}-\mathbf{y}{)}^T{{\mathbf{W}}^{1/2}}^T{\mathbf{W}}^{1/2}(\mathbf{x}-\mathbf{y})=||\mathbf{T}(\mathbf{x}-\mathbf{y})|{|}_2$$

$$r_S=\frac{\mathbb{E}[\mathrm{CSIM}(\mathbf{x},{\mathbf{y}}_1)]}{\mathbb{E}[\mathrm{CSIM}(\mathbf{x},{\mathbf{y}}_2)]}=\frac{\mathbb{E}[\mathrm{CSIM}({\mathbf{e}}_1)]}{\mathbb{E}[\mathrm{CSIM}({\mathbf{e}}_2)]}=\frac{\mathbb{E}[{\mathbf{e}}_1^T\mathbf{W}{\mathbf{e}}_1]}{{\mathbf{e}}_2^T\mathbf{W}{\mathbf{e}}_2}$$

$$\mathbb{E}[{\mathbf{e}}_1^T\mathbf{W}{\mathbf{e}}_1]$$

$$=\mathrm{Trace}(\mathbb{E}[{\mathbf{e}}_1{\mathbf{e}}_1^T]\mathbf{W})=a^2\mathrm{Trace}(\mathbf{W})$$

$$\frac{\mathbb{E}[{\mathbf{e}}_1^T\mathbf{W}{\mathbf{e}}_1]}{{\mathbf{e}}_2^T\mathbf{W}{\mathbf{e}}_2}$$

$$\underset{\mathbf{s}}{\min }\mathrm{CSIM}(\mathbf{D}\mathbf{s},\mathbf{y})$$

$$=||{\mathbf{W}}^{1/2}(\mathbf{D}\mathbf{s}-\mathbf{y})|{|}_2^2=||{\mathbf{D}}^{\prime }\mathbf{s}-{\mathbf{y}}^{\prime }|{|}_2^2$$

$$\kappa (\tilde{\mathbf{D}})=\frac{{\sigma }_{\max }(\tilde{\mathbf{D}})}{{\sigma }_{\min }(\tilde{\mathbf{D}})}$$

$$1<\frac{k_2}{k_1}\le \frac{{\kappa }_{\max }-\nu }{\xi }$$

$$\underset{\mathbf{s}}{\min }||\mathbf{s}|{|}_1,\text{s.t.}||{\mathbf{D}}^{\prime }\mathbf{s}-{\mathbf{y}}^{\prime }|{|}_2^2\le ϵ$$

$$1$$

$$n$$

$$2<2k$$

$$C_1$$

$$C_2$$

$$n_{\max }$$

$$\frac{k_2}{k_1}=\min \left\{\frac{C_1}{C_2-{\delta }_{2k}^{\prime }},\frac{{\kappa }_{\max }-\nu }{\xi }\right\}$$

$$\min_{{\bf{x}},{\bf{s}}}||{\bf{s}}||_{1}\quad\text{s.t.}\quad\left\{\begin{array}[]{l}||\mathbf{M}{\bf{x}}-{\bf{y}}||_{2}\leq\epsilon_{n}\\ {\bf{x}}=\mathbf{D}{\bf{s}}\end{array}\right.$$

$$\underset{\mathbf{x},\mathbf{s}}{\min }||\mathbf{M}\mathbf{x}-\mathbf{y}|{|}_2^2+\alpha ||\mathbf{s}|{|}_1\text{s.t.}\mathbf{x}=\mathbf{D}\mathbf{s}$$

$$\max_{{\bf{x}},{\bf{s}}}\mathrm{SSIM}({\bf{x}},\mathbf{D}{\bf{s}})\quad\text{s.t.}\quad\left\{\!\!\begin{array}[]{ll}\mathbf{M}{\bf{x}}={\bf{y}}\\ ||{\bf{s}}||_{0}\leq T\end{array}\right.$$

$$\min_{{\bf{x}},{\bf{s}}}{\mathrm{CSIM}}(\mathbf{M}\mathbf{D}{\bf{s}},{\bf{y}})\quad\text{s.t.}\quad\left\{\!\!\begin{array}[]{ll}||\mathbf{M}\mathbf{D}{\bf{s}}-{\bf{y}}||<\epsilon_{n}\\ ||{\bf{s}}||_{1}\leq T\end{array}\right.$$

$$\underset{\mathbf{x},\mathbf{s}}{\min }\mathrm{CSIM}(\mathbf{MD}\mathbf{s},\mathbf{y})+\alpha ||\mathbf{s}|{|}_1+\gamma ||\mathbf{MD}\mathbf{s}-\mathbf{y}|{|}_2^2$$

$$\min_{{\bf{x}},{\bf{s}},{\bf{z}}}{\mathrm{CSIM}}({\bf{z}})+\alpha||{\bf{s}}||_{1}+\gamma||{\bf{z}}||_{2}^{2}\quad\text{s.t.}\quad\left\{\!\!\begin{array}[]{ll}{\bf{x}}=\mathbf{D}{\bf{s}}\\ {\bf{z}}=\mathbf{M}{\bf{x}}-{\bf{y}}\end{array}\right.$$

$$\underset{\mathbf{x},\mathbf{s},\mathbf{z}}{\min }\mathcal{L}(\mathbf{x},\mathbf{s},\mathbf{z})$$

$$\displaystyle+{\mbox{\boldmath$\mu$}}_{1}^{T}({\bf{x}}-\mathbf{D}{\bf{s}})+\frac{\sigma_{1}}{2}||{\bf{x}}-\mathbf{D}{\bf{s}}||_{2}^{2}$$

$$\displaystyle+{\mbox{\boldmath$\mu$}}_{2}^{T}({\bf{z}}-\mathbf{M}{\bf{x}}+{\bf{y}})+\frac{\sigma_{2}}{2}||{\bf{z}}-\mathbf{M}{\bf{x}}+{\bf{y}}||_{2}^{2}$$

$${\mathbf{x}}^{(t+1)}$$

Full text

I Introduction

The missing sample recovery problem arises in many applications in the literature of signal processing [1, 2, 3]. It is also known as inpainting in the context of audio and image processing. Audio inpainting is investigated in [4], while image inpainting is discussed in [5] and [6].

Among numerous algorithms for missing sample recovery and inpainting, some of them exploit the sparsity of the signals [2, 7]. In this paper, we restrict ourselves to these sparse representation based algorithms. In this regard, it is assumed that the signal is sparse within a domain such as Discrete Cosine Transform (DCT), Discrete Wavelet Transform (DWT) [8] or any other predefined or learned complete or overcomplete dictionary. The sparsity of the signal on a dictionary-based representation, means that the vector of coefficients of the signal in the transform domain has many zeros (or nearly zeros) and only a few of its elements are nonzero. Neglecting the insignificant (zero) coefficients, it is possible to reconstruct the signal with considerably low error. The sparsity of the signal gives us the ability to reconstruct it from very few random measurements far below the Nyquist rate. This is well known as Compressed Sensing (CS) [9, 10], which has had many applications in the past decade [11]. The problem of reconstruction of the sparse signal from a few random measurements, is also known as sparse recovery. Many algorithms are proposed for sparse recovery of signals in different applications in audio and image processing.

In a missing sample recovery problem, some samples of the signal are missed due to physical impairment, unavailability of measurements, or distortion and disturbances. In such cases, it is shown that the corrupted samples would better be omitted throughout the reconstruction process [2]. Thus, the discarded samples may be considered as missed. Even with these missing samples, the signal can still be reconstructed, given the sparsifying basis or dictionary and the corresponding sparse coefficients. Many algorithms and optimization problems are suggested to recover the sparse samples in this regard. The fundamental problem in a sparse recovery method is to maximize the sparsity which is principally stated in terms of the $\ell_{0}$-norm. There are a class of greedy algorithms for strictly sparse signal recovery based on $\ell_{0}$-norm minimization. These include Matching Pursuit (MP) [12], Orthogonal Matching Pursuit (OMP) [13], Regularized OMP (ROMP) [14], Compressive Sampling Matching Pursuit (CoSaMP) [15] and Generalized Matching Pursuit (GOMP) [16]. There are also iterative methods based on majorization minimization technique proposed for approximate $\ell_{0}$-norm minimization using surrogate functions. The Iterative Hard Thresholding (IHT) algorithm [17] is the first member of this class. There is also a modified version which uses adaptive thresholding named as Iterative Method with Adaptive Thresholding (IMAT) [18] with different variants including IMATI and IMATCS [19]. A recent improved version of this method called INPMAT is also proposed in [20]. Furthermore there is an approach for sparse approximation based on Smoothed-$\ell_{0}$ (SL0) norm minimization presented in [21]. The $\ell_{0}$-minimization algorithms are mostly used in cases where the signal has exactly sparse support and the sparsity is known. But in many practical situations the sparsity is unknown or the signal is not strictly sparse but instead compressible, meaning that most coefficients are negligible (despite being precisely zero) compared to the significant elements. A good and a common alternative is to use the $\ell_{1}$-norm as the nearest convex approximation of $\ell_{0}$-norm. This approach is called $\ell_{1}$-minimization or the basis-pursuit method [22]. There are many algorithms presented for $\ell_{1}$-norm minimization including Iterative Soft Thresholding Algorithm (ISTA) [23] and the fast version FISTA [24], $\ell_{1}$ Least Squares (L1-LS) [25], Primal and Dual Augmented Lagrangian Methods (PALM and DALM)[26], Iterative Bayesian Algorithm (IBA) [27], Sparse Bayesian Learning (SBL) [28] and Bayesian Compressed Sensing (BCS) [29]. There are also more general $p$-norm minimization based algorithms available for solving the sparse recovery problem in the literature [30]. For detailed survey on sparse recovery methods, one can refer to [31].

In this paper, we propose an alternative $\ell_{1}$-minimization method for sparse recovery of signals. In particular, we consider the sparse recovery of image patches with missing samples which has application in image inpainting and restoration. We use a criterion for measuring the similarity between two image signals, called Convex SIMilarity (CSIM) Index [32]. Although it is derived from the Structural SIMilarity (SSIM) index, the well-known perceptual quality assessment criterion [33], it has feasible mathematical features unlike its predecessor. These properties including convexity and error-sensitivity, which are investigated in this paper, result in simplified methods for solving the optimization problem incorporating CSIM index as the measure of similarity. We also present analyses on how to choose optimal parameters of the similarity index proposed in [32]. In this paper, we use CSIM criterion in our proposed optimization problem for missing sample recovery. Similar to [26], an iterative algorithm is presented for $\ell_{1}$-minimization based on Alternating Direction Method of Multipliers (ADMM) method. Simulation results show the efficiency of the proposed method called CSIM minimization via Augmented Lagrangian Method (CSIM-ALM) compared to some popular existing algorithms.

II Metrics for Image Quality Assessment

There are different criteria for Image Quality Assessment (IQA) which can be used as measures of similarity between two signals. The most popular fidelity metric for measuring the similarity between two signals ${\bf{x}},{\bf{y}}\in\mathbb{R}^{n}$ is MSE which is defined as $\mathrm{MSE}({\bf{x}},{\bf{y}})=\frac{1}{n}||{\bf{x}}-{\bf{y}}||^{2}_{2}$. This criterion is widely used, because MSE or equivalently $\ell_{2}$-norm is mathematically a well-defined function of the difference between the reference and the test signal. This function has desirable mathematical features such as convexity and differentiability. Hence, the optimization problem is simple and tractable using convex programming methods. Nevertheless, there are cases in which the MSE criterion seems to be inefficient to accurately recover the original image signal, especially in the presence of noise. One reason is that this fidelity metric is indifferent toward the error distribution, i.e., the statistics of the error signal. For instance consider two scenarios in which a signal is corrupted. In the first case the reference signal is perturbed by noise and in the second, a constant amplitude is added to the original image. Both corrupted images have the same MSE distance with respect to the original signal, whilst the noisy image is definitely more visually deteriorated. Thus, there are a class of perceptual criteria introduced for measuring visual quality of images. The most popular perceptual metric is SSIM which is defined as [33]:

[TABLE]

where $C_{1},C_{2}>0$ are constant and $\mu$ and $\sigma$ denote the mean and the variance (cross-covariance), respectively. This function whose mathematical properties are discussed in [34], is non-convex and multi-modal implying that the problem of optimizing this criterion is hard to solve.

II-A The CSIM index

As mentioned earlier, SSIM function is non-convex and multi-modal which results in hard optimization problems. In [32], we have proposed a simplified criterion which inherits error-sensitivity property of SSIM. The proposed index called CSIM, is defined as follows:

[TABLE]

where $k_{1}$ and $k_{2}$ are positive constants and $\mu$ and $\sigma$ denote statistical estimates for mean and variance/cross-covariance. It is assumed that $k_{1}<k_{2}$, to ensure biased sensitivity toward random disturbance or noise compared to uniform level change. In fact unlike MSE, the proposed criterion has noise-sensitive variation. In other words, a constant change in the brightness level of the image does not alter this criterion as much as noise. This is logically because constant change in the amplitude only affects the signal mean value; it does not change variance/cross-covariance. Thus, as far as $k_{1}<k_{2}$, the CSIM index is only slightly influenced by mere brightness level change. This index also benefits some feasible mathematical features, including convexity and differentiability.

Theorem 1**.**

For positive $k_{1}$ and $k_{2}$, the fidelity criterion defined by equation (2) is strictly convex with respect to ${\bf{x}}$ or ${\bf{y}}$. Furthermore, there exists a transformation matrix $\mathbf{T}$, such that the proposed criterion can be restated in terms of $\ell_{2}$-norm, after the transformation.

Proof:

Suppose ${\bf{x}},{\bf{y}}\in\mathbb{R}^{n}$. If we use unbiased estimate for variance and covariance, we will have:

[TABLE]

where {\mbox{\boldmath\mu}}_{x}=\mu_{x}{\bf{1}}_{n} and ${\bf{1}}_{n}=(1,\ldots,1)^{T}\in\mathbb{R}^{n}$. Therefore, \sigma_{x}^{2}+\sigma_{y}^{2}-2\sigma_{x,y}=\frac{1}{n-1}||({\bf{x}}-{\mbox{\boldmath\mu}}_{x})-({\bf{y}}-{\mbox{\boldmath\mu}}_{y})||_{2}^{2} and $\mu_{x}^{2}+\mu_{y}^{2}-2\mu_{x}\mu_{y}$ equals $(\mu_{x}-\mu_{y})^{2}$. Hence, if we define the error signal as the difference between test signals, i.e., ${\bf{e}}={\bf{x}}-{\bf{y}}$, equation (2) is simplified to:

[TABLE]

which implies ${\mathrm{CSIM}}$ is always non-negative. Likewise, if we consider ${\bf{x}}$ and ${\bf{y}}$ as samples of random variables $X$ and $Y$, we have ${\mathrm{CSIM}}(X,Y)=k_{1}\mathbb{E}^{2}[X-Y]+k_{2}\mathrm{\mathbf{var}}[X-Y]\geq 0$ which equals zero if and only if $P(X=Y)=1$. Now, using simplifications in (3), one can expand equation (4) as follows:

[TABLE]

where $\mathbf{M}=\mathbf{M}^{T}\mathbf{M}=\mathbf{I}_{n}-\frac{1}{n}{\bf{1}}_{n}{\bf{1}}_{n}^{T}$ and $\mathbf{I}_{n}$ is the identity matrix of size $n$. Therefore, the matrix $\mathbf{W}\in\mathbb{R}^{n\times n}$ is obtained by:

[TABLE]

Since ${\mathrm{CSIM}}$ is a quadratic function of ${\bf{e}}$, the positive-definiteness of the Hessian matrix $\mathbf{W}$, implies its strict convexity as well. On the other hand, since ${\bf{e}}={\bf{x}}-{\bf{y}}$ is affine, it will also be concluded that ${\mathrm{CSIM}}$ is strictly convex with respect to both ${\bf{x}}$ and ${\bf{y}}$. Now, to prove the second part of the theorem, we proceed as follows. Since $\mathbf{W}$ is Hermitian and positive-definite, it is diagonalizable and thus, $\mathbf{W}^{1/2}$ exists:

[TABLE]

Hence, we can write:

[TABLE]

which means there exists a transformation denoted by $\mathbf{T}=\mathbf{W}^{1/2}$, after which, we can restate the proposed index in terms of $\ell_{2}$-norm distance. ∎

II-B The choice of the parameters $k_{1}$ and $k_{2}$

Now, the problem is how to find $k_{1}$ and $k_{2}$ so that the proposed criterion has utmost sensitivity toward random perturbation compared to uniform change.

Definition 1**.**

Assume ${\bf{e}}_{1}\in\mathbb{R}^{n}$ is a random binary signal with i.i.d. elements taking values $\{a,-a\}$ with equal probability, i.e., $C_{{\bf{e}}_{1}}=\mathbb{E}[{\bf{e}}_{1}{\bf{e}}_{1}^{T}]=a^{2}\mathbf{I}_{n}$, and ${\bf{e}}_{2}=a{\bf{1}}_{n}$ is a deterministic constant-amplitude signal. Let ${\bf{y}}_{1}={\bf{x}}+{\bf{e}}_{1}$ and ${\bf{y}}_{2}={\bf{x}}+{\bf{e}}_{2}$ where ${\bf{x}}\in\mathbb{R}^{n}$ is the reference signal. We define the ratio of sensitivity as:

[TABLE]

Proposition 1**.**

The sensitivity ratio with respect to CSIM criterion with kernel $\mathbf{W}$, is obtained by $r_{S}=\dfrac{\sum_{i}w_{i,i}}{\sum_{i}\sum_{j}w_{i,j}}$.

Proof:

Since ${\bf{e}}_{1}^{T}\mathbf{W}{\bf{e}}_{1}$ is scalar and $\mathbb{E}$ and $\mathrm{Trace}$ are linear operators, we may write:

[TABLE]

Thus:

[TABLE]

∎

Now for $\mathbf{W}$ defined in (85), the sensitivity ratio is $r_{S}=\frac{k_{2}}{k_{1}}+\frac{1}{n}$. This states that the greater the ratio $k_{2}/k_{1}$ is, the more sensitive to noise, the proposed CSIM index will be. But there are other conditions which impose constraints on the eigenvalues of $\mathbf{W}$ and accordingly the values of $k_{1}$ and $k_{2}$.

II-B1 Condition number

Consider the optimization problem below with $\mathbf{D}\in\mathbb{R}^{n\times p}$:

[TABLE]

where $\mathbf{D}^{\prime}=\mathbf{W}^{1/2}\mathbf{D}$ and ${\bf{y}}^{\prime}=\mathbf{W}^{1/2}{\bf{y}}$. If $\mathbf{D}^{\prime}$ is full-column rank ($n>p$), the solution to (12) will be obtained by $\mathbf{D}^{\prime\dagger}{\bf{y}}^{\prime}=(\mathbf{D}^{T}\mathbf{W}\mathbf{D})^{-1}\mathbf{D}^{T}\mathbf{W}{\bf{y}}=\tilde{\mathbf{D}}{\bf{y}}$111In case where $\mathbf{D}^{\prime}$ is full-row rank ($n<p$), we have $\mathbf{D}^{\prime\dagger}{\bf{y}}^{\prime}=\mathbf{D}^{\dagger}{\bf{y}}$. Thus, the matrix $\mathbf{W}$ does not change the condition number. To have a robust (reliable) solution, the matrix $\tilde{\mathbf{D}}$ must not be ill-conditioned. Now assuming $\mathbf{D}$ or equivalently $\tilde{\mathbf{D}}$ is full-rank, using the extended definition of condition number for non-square matrices, we can write:

[TABLE]

where $\sigma_{\min}(\tilde{\mathbf{D}})$ denotes the minimum (non-zero) singular value of $\tilde{\mathbf{D}}$. The following theorem states how the constraint on the maximum value of $\kappa(\tilde{\mathbf{D}})$ corresponds to the ratio of the parameters $k_{1}$ and $k_{2}$.

Theorem 2**.**

Define $\xi=\kappa(\mathbf{D})\left(\frac{n}{n-1}\right)\left(\frac{1}{\kappa^{2}(\mathbf{D})}-\frac{1}{n}\frac{{\bf{1}}^{T}\mathbf{D}\mathbf{D}^{T}{\bf{1}}}{\sigma^{2}_{\max}(\mathbf{D})}\right)$ and $\nu=\kappa(\mathbf{D})\left(\frac{1}{n}\frac{{\bf{1}}^{T}\mathbf{D}\mathbf{D}^{T}{\bf{1}}}{\sigma^{2}_{\max}(\mathbf{D})}\right)$. The matrix $\tilde{\mathbf{D}}$ satisfies $\kappa(\tilde{\mathbf{D}})\leq\kappa_{\max}$ with $\kappa_{\max}>\xi+\nu$, if:

[TABLE]

Proof:

Refer to Appendix A. ∎

II-B2 RIP condition

Another constraint is imposed when using the proposed criterion in a $\ell_{1}$-minimization problem. In particular, consider the Basis Pursuit Denoising (BPDN) problem [35] with CSIM as the fidelity criterion:

[TABLE]

In order for (15) to uniquely recover $k$-sparse signals, the matrix $\mathbf{D}^{\prime}$ should satisfy RIP of order $2k$ with constant $\delta_{2k}<\sqrt{2}-1$. Assume the columns of $\mathbf{D}$ have unit norm and let $\mu(\mathbf{D})$ denote the mutual coherence between these column vectors. Investigating the RIP condition for $\mathbf{D}^{\prime}$, yields the following theorem:

Theorem 3**.**

Assume the columns of $\mathbf{D}$ are normalized. The matrix $\mathbf{D}^{\prime}=\mathbf{W}^{1/2}\mathbf{D}\in\mathbb{R}^{n\times p}$ satisfies RIP condition of order $2k$ with constant $\delta^{\prime}_{2k}<\sqrt{2}-1$, if the following conditions are satisfied:

[TABLE]

where

[TABLE]

.

Proof:

Refer to Appendix B for the proof. ∎

Corollary 1**.**

To gain the maximum ratio of sensitivity, while satisfying the condition number and the RIP restrictions, we choose:

[TABLE]

where in practice, we choose $\kappa_{\max}=4$ and $\delta^{\prime}_{2k}=0.4$. Furthermore, in case where the sparsity of the signal is unknown, we assume 10% sparsity, i.e., $k=\lfloor 0.1n\rfloor$. For more details go to section V.

III Statement of The Problem

The problem of recovering an image with samples missed at random, is equivalent to random sampling reconstruction of a signal. This problem is also addressed in the literature as block loss restoration due to error in the transmission channel [36, 37]. Or alternatively, known as image inpainting in applications where the sampling mask is known and the objective is to fill in the gaps or remove occlusion or specific objects from the image [38, 6, 39]. Suppose ${\bf{x}}\in\mathbb{R}^{N}$ is the vectorized image signal and $\mathbf{M}$ is the sampling matrix by which the pattern of sampling of the image signal is determined. In other words $\mathbf{M}\in\mathbb{R}^{N\times N}$ is the identity matrix with $N-m$ of its diagonals set to 0, where $m$ denotes the number of available samples. The observed image signal with missing samples, is also denoted by ${\bf{y}}\in\mathbb{R}^{N}$. Among many approaches for missing sample recovery of images, there are also a class of inpainting algorithms which use sparse representation for image restoration. In [5, 19, 40, 41], there are iterative algorithms proposed for recovery of missing samples exploiting the sparsity of representation based on redundant (over-complete) dictionaries. If we assume that ${\bf{x}}$ has approximately a sparse representation based on the atoms of a dictionary specified by the matrix $\mathbf{D}$, the regular optimization problem for sparse recovery of the missing samples is formulated as follows:

[TABLE]

where ${\bf{s}}$ denotes the sparse vector of representation coefficients and $\epsilon_{n}$ denotes the variance of the additive observation noise. This problem is indeed the extension of the BPDN problem for missing sample recovery. Using a proper value of $\alpha$, this problem is also shown to be equivalent to [35]:

[TABLE]

Also, there exist methods optimized in terms of different image quality measures. In [42], an exemplar-based method for image completion is proposed which uses dictionary learning based on SSIM for spare recovery of local image patches. In this algorithm, for each candidate patch, the optimization problem below is solved within the sparse coding step:

[TABLE]

This problem is iteratively solved using a matching pursuit approach, i.e., in each step, the support of the sparse vector is retrieved and the coefficients are subsequently obtained solving unconstrained (20). Although this problem addresses perceptual quality of reconstruction, it is non-convex and thus solved using time-consuming linear search methods. In this paper, we propose an alternative method for sparse recovery of image signals which can be applied in the sparse coding step of an adaptive dictionary learning method for image inpainting. We also use a perceptual metric, instead of $\ell_{2}$-norm, for visually enhanced reconstruction of the missing samples.

IV The Proposed Algorithm

As discussed earlier in section III, most of the algorithms use $\ell_{2}$-norm as fidelity criterion for image reconstruction. But, there are also inpainting methods based on local sparse representation which use perceptual IQA metrics namely SSIM, for recovery of the missing samples. Here, we propose to use CSIM instead of SSIM. Hence, to solve the missing sample recovery problem defined in (18), we incorporate our proposed perceptual metric for reconstruction of the image samples. In fact, among all the solutions with equal $\ell_{2}$-norm distance from ${\bf{y}}$, we seek for those with maximum similarity in terms of ${\mathrm{CSIM}}$. Hence, the proposed optimization problem is as follows:

[TABLE]

which is equivalent to

[TABLE]

where $\alpha_{1}$ and $\gamma$ are chosen such that Karush-Kuhn-Tucker (KKT) conditions are satisfied [43]. Note that since ${\mathrm{CSIM}}$ is convex and uni-modal, finding the local minima of (21) is sufficient to attain its global optimum. Now, introducing auxiliary variables ${\bf{x}}=\mathbf{D}{\bf{s}}$ and ${\bf{z}}=\mathbf{M}{\bf{x}}-{\bf{y}}$ the optimization problem (22) would change to:

[TABLE]

The auxiliary variables are used to separate the main problem into simpler sub-problems. Since the CSIM function is convex, it is guaranteed to use the Alternating Direction Method of Multipliers (ADMM) [44] to solve (II-B). Hence, the final cost function to be optimized is the augmented Lagrangian function:

[TABLE]

The ADMM alternatively minimizes (IV) with respect to each variable while assuming the other variables fixed. Hence, at each iteration of the ADMM, the problem (IV) is split into three sub-problems as follows:

IV-A ${\bf{x}}$* sub-problem:*

The augmented Lagrangian cost function with respect to ${\bf{x}}$ is a quadratic function. Hence, the optimization sub-problem associating with ${\bf{x}}$ at $t$-th iteration of the ADMM is:

[TABLE]

where $\mathbf{Q}=\sigma_{2}\mathbf{M}^{T}\mathbf{M}+\sigma_{1}\mathbf{I}$ and {{\bf{b}}^{(t)}}=\sigma_{1}\mathbf{D}{\bf{s}}^{(t)}-{\mbox{\boldmath\mu}}_{1}^{(t)}+\mathbf{M}^{T}(\sigma_{2}({\bf{z}}^{(t)}+{\bf{y}})+{\mbox{\boldmath\mu}}_{2}^{(t)}). The solution to this problem is simply obtained by ${\bf{x}}^{(t+1)}=\left(\sigma_{1}\mathbf{I}+\sigma_{2}\mathbf{M}^{T}\mathbf{M}\right)^{-1}{\bf{b}}^{(t)}$. Now, since $\mathbf{M}\mathbf{M}^{T}=\mathbf{I}$, using Sherman-Morrison-Woodbury lemma [45], the solution to (25) is simplified to:

[TABLE]

IV-B ${\bf{s}}$* sub-problem:*

The optimization sub-problem associating with ${\bf{s}}$ is:

[TABLE]

Using the Majorization Minimization (MM) technique [46], as proposed in [23], this problem is shown [32] to be equivalent to minimization of a surrogate cost function denoted by ${\cal L}^{S}({\bf{s}},{\bf{s}}_{0})$:

[TABLE]

where {\bf{a}}({\bf{s}}_{0})={\bf{s}}_{0}+\frac{1}{\lambda}\mathbf{D}^{T}\big{(}{\bf{x}}^{(t+1)}-\mathbf{D}{\bf{s}}_{0}+\frac{1}{\sigma_{1}}{\mbox{\boldmath\mu}}_{1}^{(t)}\big{)} and $\lambda>\lambda_{\max}(\mathbf{D})$. Let us set ${\bf{s}}_{0}={\bf{s}}^{(t)}$, where $t$ denotes the iteration number. Now the solution to (28) is obtained using the soft-thresholding operator and ${\bf{s}}^{(t)}$ is updated according to:

[TABLE]

where ${\cal S}_{\tau}({\bf{x}})$ is a vector with component ${\cal S}_{\tau}({\bf{x}})_{i}$ obtained by:

[TABLE]

IV-C ${\bf{z}}$* sub-problem*

The sub-problem associating with ${\bf{z}}$ is as follows:

[TABLE]

Now substituting ${\mathrm{CSIM}}$ from (2), the resulting cost function will be in the quadratic form below:

[TABLE]

where $\mathbf{K}=\sigma_{2}\mathbf{I}+2(\mathbf{W}+\gamma\mathbf{I})=\zeta_{1}\mathbf{I}_{n}+\zeta_{2}{\bf{1}}_{n}{\bf{1}}_{n}^{T}$ with $\zeta_{1}=\sigma_{2}+2\frac{k_{2}}{n-1}+2\gamma$ and $\zeta_{2}=2\left(\frac{k_{1}}{n^{2}}-\frac{k_{2}}{n(n-1)}\right)$ and {{\bf{c}}^{(t)}}=\sigma_{2}(\mathbf{M}{\bf{x}}^{(t+1)}-{\bf{y}})-{\mbox{\boldmath\mu}}_{2}^{(t)}. The solution to (32) is obtained by:

[TABLE]

where in the latter, we have applied the matrix inverse lemma.

IV-D Multipliers update

The final step of the ADMM is to update the Lagrangian multipliers associated with the equality constraints. Hence, we have:

[TABLE]

Theorem 4**.**

If $\lambda>\lambda_{\max}(\mathbf{D})$ is satisfied and the parameters of CSIM, i.e., $k_{1},k_{2}$ are positive, then, the proposed algorithm is guaranteed to converge to the optimal solution of problem (22) for fixed positive $\alpha$, $\sigma_{1}$, $\sigma_{2}$ and $\gamma$

Proof:

The proof of convergence is given in Appendix C. ∎

IV-E Remark

We can enhance the performance of the proposed method using two techniques, namely backtracking and continuation. First of all, since ${\bf{s}}^{(t+1)}$ is the minimizer of ${\cal L}^{S}({\bf{s}},{\bf{s}}^{(t)})$, we have:

[TABLE]

where the first inequality implies from the fact that the surrogate function is the majorization of the original cost function. To ensure this condition is fulfilled, we may use a backtracking procedure to choose the appropriate value of $\lambda$. This method as proposed in [24], solves the optimization problem (28) and checks whether the solution ${\bf{s}}^{\ast}$ satisfies ${\cal L}({\bf{s}}^{\ast})\leq{\cal L}^{S}({\bf{s}}^{\ast},{\bf{s}}^{(t)})$. If true, the value of ${\bf{s}}^{(t+1)}$ is set to ${\bf{s}}^{\ast}$ and if not, it multiplies the value of $\lambda$ by a constant $\beta>1$. We may also improve the rate of convergence, using continuation on the regularizing parameter $\alpha$, i.e., we use a decreasing sequence for $\alpha^{(t)}$ instead of a fixed regularizing parameter. To this purpose, we choose exponentially decreasing $\alpha^{(t)}$ which alike [11], results in an adaptive approach to specify the threshold $\frac{\alpha}{\lambda\sigma_{1}}$. We also substitute corresponding samples of the reconstructed image ${\bf{x}}^{(t+1)}$ with known samples of the observed signal at each iteration. There is of course a minimum threshold for $\alpha^{(t)}$ specified by the precision tolerance of the final solution. Thus, the enhanced version of the proposed method with $\alpha$ continuation is shown in Alg. 1.

V Simulation Results

V-A Denoising

In this part, we conduct an experiment to show the performance of the proposed quality assessment criterion compared to some popular criteria, namely MSE, SSIM and FSIM [47]. Consider ${\bf{x}}\in\mathbb{R}^{N}$ is the reference image signal and ${\bf{y}}\in\mathbb{R}^{N}$ denotes the noisy observed image, i.e., ${\bf{y}}={\bf{x}}+{\bf{n}}$ where ${\bf{n}}\in\mathbb{R}^{N}$ is the noise signal with distribution ${\bf{n}}\sim{\cal{N}}(0,\sigma_{n}^{2})$. Suppose that the image is divided into small patches of size $\sqrt{n}\times\sqrt{n}$. The problem is to find a linear denoising filter ${\bf{h}}_{j}\in\mathbb{R}^{m}(m\leq n)$ whose convolution with the $j$th patch of the image denoted by ${\bf{y}}_{j}\in\mathbb{R}^{n}$ gives an estimate of the original signal denoted by $\hat{{\bf{x}}}_{j}$. i.e.

[TABLE]

After the small patches are denoised, the entire image is then reconstructed by superposition of the recovered patches.

In this scenario, the original image signal is unknown. In fact, ${\bf{x}}$ is the vector of spatial samples of a random process $x[i]$, which is assumed to be ergodic and WSS stationary. The clean image is then perturbed by a white Gaussian noise process $n[i]$ with variance $\sigma_{n}^{2}$, whose samples are denoted by ${\bf{n}}$. The observed noisy signal ${\bf{y}}$ is thus, modelled by the sum of theses two processes. Hence, in this case, the problem of finding the equalizer filter, is indeed a linear estimation problem as discussed in [48]. But, here we assume stationariness within the spatial domain of each patch. Consequently, we confine ourselves to local patch denoising for estimation of ${\bf{h}}_{j}$. The usual fidelity criterion for estimation is MSE, which leads to the so-called Wiener-Hopf [49] equations:

[TABLE]

where $\mathbf{R}_{y_{j},y_{j}}$ and ${\bf{r}}_{x_{j},y_{j}}$ denote the auto-correlation matrix and the the vector of cross-correlation components respectively, i.e., $\mathbf{R}_{y_{j},y_{j}}[k,l]=r_{y_{j},y_{j}}[k-l],\quad{\bf{r}}_{x_{j},y_{j}}[k]=r_{x_{j},y_{j}}[k]$. Now, assume the noise process is independent from ${\bf{x}}$ and is distributed homogeneously across the whole image. Using ${\bf{y}}_{j}={\bf{x}}_{j}+{\bf{n}}_{j},\,{\bf{n}}_{j}\bot{\bf{x}}_{j}$ with $\mu_{n_{j}}=0$ and $c_{n_{j},n_{j}}[k]=\sigma_{n}^{2}\delta[k]$, it can be shown that [48]

[TABLE]

Instead of the MSE criterion, we may use CSIM or SSIM in our estimation problem. If we use the statistical definition of CSIM, the optimization problem for finding the denoising filter ${\bf{h}}_{j}$ would be:

[TABLE]

Substituting ${\hat{x}_{j}}$ from (36) we have the following relations:

[TABLE]

If we denote the covariance matrix of ${\bf{y}}_{j}$ by $\mathbf{C}_{y_{j},y_{j}}$, the optimization problem (39) is simplified to:

[TABLE]

This problem is quadratic in terms of ${\bf{h}}_{j}$. To solve this, we use (38) and we differentiate with respect to ${\bf{h}}_{j}$ which gives:

[TABLE]

To reduce the complexity of calculating the inverse above, having $\mathbf{C}_{y_{j},y_{j}}^{-1}$, we can use matrix inverse lemma.

Now, let us turn to SSIM fidelity criterion. The SSIM optimization problem for estimating the optimal equalization filter, tries to maximize the cost function below with respect to ${\bf{h}}_{j}$:

[TABLE]

This optimization problem is solved via conversion to a quasi-convex cost, using bi-sectional search method [48].

In practical simulations for denoising image patches, the covariance matrix $\mathbf{C}_{y_{j},y_{j}}$ is empirically estimated using unbiased estimation and we use equation (38) for the values of $\sigma_{x_{j}}$ and $\mu_{x_{j}}$. The simulation results for denoising with FIR filters are given in table I. It is clear that the proposed CSIM index is outperforming other criteria in terms of estimation quality at low SNR. It is even faster in several cases and generally more efficient than MSE and SSIM in the sense of FSIM performance as a tertiary index.

V-B Sparse recovery

In this experiment, we compare the quality performance of the proposed CSIM-ALM method for recovery of missing samples of image patches with some popular sparse recovery algorithms. We use IMATCS222http://ee.sharif.edu/$\sim$imat/, L1-LS, DALM333https://people.eecs.berkeley.edu/$\sim$yang/software/l1benchmark/, TV444http://www.caam.rice.edu/$\sim$optimization/L1/TVAL3/ [50], FISTA, SL0555http://ee.sharif.edu/$\sim$SLzero/, GOMP666http://islab.snu.ac.kr/paper/gOMP.zip, and the method in [42], which we call SSIM-based Matching Pursuit (SSIM-MP).

For simulations of this part, we use $8\times 8$ image patches vectorized using raster scanning, and we select 100 patch vectors of size $n=64$ denoted by ${\bf{x}}_{i}$ at random. For each patch, a binary random sampling mask $\mathbf{M}_{j}$ with $m_{j}$ ones, is generated and the observed image signal for each experiment is acquired by ${\bf{y}}_{i,j}=\mathbf{M}_{j}\odot{\bf{x}}_{i}$ where $\odot$ denotes pointwise (Hadamard) product. The locations of $1$s are chosen uniformly at random and the sampling ratio of the signal defined as $SR_{j}=\frac{m_{j}}{n}$ varies between $(0,1)$. We use complete ($64\times 64$) and over-complete ($64\times 128$) DCT and Haar Wavelet Packet (WP)777MATLAB commands wmpdictionary and wphaar dictionaries for sparse representation.

Natural image signals

In this scenario, we extract $8\times 8$ patches from natural grayscale images. Thus, the exact sparsity is unknown. Hence, to use matching pursuit algorithms, we consider $10\%$ sparsity. After recovery of missed samples, we then average over random experiments (random ${\bf{x}}_{i}$s and random masks with same $m_{j}$) and plot the PSNR and SSIM, versus the sampling rate. The parameters for SL0 and TV are set to their defaults. The values of the parameters for IMATCS are set to $\alpha=0.2$, $\beta=0.5||\mathbf{D}^{T}\mathbf{M}^{T}{\bf{y}}||_{\infty}$, $\lambda_{\min}=1e-3$. The stopping criterion for DALM and L1-LS are set to their defaults, meaning that the algorithms stop when the duality gap falls below a certain tolerance. The stopping criterion for the remaining algorithms including IMATCS, FISTA, SL0 and CSIM-ALM is set to the maximum iteration count which is 50. The parameters for CSIM-ALM are chosen as $\sigma_{1}=0.4\frac{m}{n}$, $\sigma_{2}=2\frac{m}{n}$, $\gamma=1$, $\xi=0.1$, $\eta=0.95$, $\beta=1.1$, $\alpha_{\min}1e-4$, $k_{2}=n-1=63$ and $k_{1}=0.25k_{2}$. The choices for $k_{1}$ and $k_{2}$ are empirically obtained considering constraints in 1. Fig. 1 shows the results of the sparse recovery from random samples. The top row figures depict reconstruction via DCT atoms and the bottom row show the results via Haar Wavelet Packet dictionary. As shown in these figures, the proposed CSIM-ALM algorithm, mostly outperforms the state of the art methods for sparse recovery via $\ell_{1}$-norm minimization specifically at $SR>0.3$. It mainly provides a better reconstruction quality compared to DALM which commonly uses the ADMM technique to solve the $\ell_{1}$ optimization problem. This superiority is particularly more apparent in terms of SSIM performance. In fact, the proposed algorithm also outperforms SSIM-MP which is based on non-convex SSIM maximization.

Artificially generated signals

In this scenario, the patch signals are artificially produced by the product of generated sparse vectors $\{{\bf{s}}_{i}\}$ with $(64\times 64)$ or $(64\times 128)$ DCT atoms as sparsifying basis. Similar to the previous experiment, we assume $10\%$ sparsity, i.e., $k=\lceil 0.1p\rceil$ where $p=64$ or $128$ denotes the number of atoms. The $k$ nonzero coefficients of sparse vectors, are selected uniformly at random and the values are chosen according to Normal distribution. For performance comparison, we consider relative error defined as:

[TABLE]

where $\hat{{\bf{s}}}_{i}$ denotes the recovered sparse vector, given the observation samples $\mathbf{M}\mathbf{D}{\bf{s}}_{i}$ where $\mathbf{M}$ denotes the random sampling mask with sampling ratios $0.4,0.6$ and $0.8$.. Fig. 2 shows these results for random sampling reconstruction of artificially generated signals versus the number of iterations. As shown in this figure, the rate of convergence of CSIM-ALM is comparably faster than IMATCS and FISTA while the final error is significantly lower. This is specifically, more obvious for recovery via $(64\times 64)$ DCT atoms. The running time of these methods for sparse reconstruction from complete DCT dictionary, is also given in Fig. 3. The vertical axis in this figure represents the time (in seconds) consumed upon each iteration, and the horizontal axis shows the number of iterations. Of course, the number of iterations for L1-LS which is determined by the stopping criterion of the algorithm, may be lass than the maximum limit. Hence, in Fig. 2 and Fig. 3 only some iterations of this method are given and beyond this limit the algorithm usually starts to diverge. According to Fig. 3, the least complex algorithms are SL0, DALM and IMATCS but they need more iterations to reach acceptable recovery performance compared to CSIM-ALM. Indeed according to Fig. 2, CSIM-ALM converges much faster than SL0 at high $SR$s and DALM at low sampling rates. This is in particular, more apparent for recovery via over-complete DCT dictionary where SL0 and L1-LS usually fail to converge.

VI Conclusion

In this paper, a fidelity metric called CSIM is introduced, which is convex and error-sensitive and can be applied for image quality assessment. The proposed index like MSE, is well suited for mathematical manipulations and like SSIM, has perceptual meanings. We investigate mathematical features of CSIM index in this paper and incorporate it as fidelity criterion for solving missing sample recovery problem based on sparse representation for the image signal. This recovery algorithm can also be applied in any sparse coding application such as dictionary learning and sparse approximation for image signals. Furthermore, an iterative ADMM-based algorithm is proposed to solve the $\ell_{1}$-minimization inverse problem for recovery of missing samples. The convexity of the optimization function leads to proof of convergence for the algorithm. Simulation results show the efficiency of CSIM and the the proposed iterative algorithm over counterpart methods for missing sample recovery of image signals.

Appendix A Proof of Condition Number (Theorem 2)

A-1 Condition number of $\mathbf{W}$

Let us assume $\mathbf{P}=\left(\frac{k_{1}}{n^{2}}-\frac{k_{2}}{n(n-1)}\right){\bf{1}}_{n}{\bf{1}}_{n}^{T}$, $\mathbf{Q}=\mathbf{W}$ and $\gamma=\frac{k_{2}}{n-1}$. Now, using lemma 1, we can obtain the eigenvalues of $\mathbf{W}$. But first, we need to compute the eigenvalues of $\mathbf{P}$. Since $\mathbf{P}$ is symmetric and rank-one, it is implied that there is only a single (real) non-zero eigenvalue denoted by $\lambda_{\mathbf{P}_{n}}$ (with the corresponding eigenvector ${\bf{v}}_{\mathbf{P}_{n}}$) which satisfies:

[TABLE]

It also implies that all the remaining (sorted) eigenvalues $\lambda_{\mathbf{P}_{i}},\,1\leq i\leq n-1$ are zero. Now, if we set ${\bf{v}}_{\mathbf{P}_{n}}={\bf{1}}_{n}$, then $\lambda_{\mathbf{P}_{n}}$ is given by:

[TABLE]

Consequently, the eigenvalues of $\mathbf{Q}$ are obtained by:

[TABLE]

Now assuming $k_{2}>k_{1}$, it yields $\kappa(\mathbf{W})=\frac{\lambda_{\max}(\mathbf{W})}{\lambda_{\min}(\mathbf{W})}=\frac{k_{2}}{k_{1}}\left(\frac{n}{n-1}\right)=\frac{\theta_{1}}{n\theta_{2}+\theta_{1}}$.

A-2 Condition number of $\tilde{\mathbf{D}}$

Assume $\mathbf{A}=\mathbf{D}^{T}\mathbf{W}\mathbf{D}\in\mathbb{R}^{p\times p}$. Let $\mathbf{A}=\mathbf{B}+\mathbf{C}$ where $\mathbf{B}=\theta_{1}\mathbf{D}^{T}\mathbf{D}$ and $\mathbf{C}=\theta_{2}\mathbf{D}^{T}{\bf{1}}{\bf{1}}^{T}\mathbf{D}$. Also assume $\mathbf{D}$ is full column rank ($p<n$) and the eigenvalues are sorted in descending order. Using Weyl’s inequality [51] for the eigenvalues of sum of Hermitian matrices $\mathbf{B}$ and $\mathbf{C}$, we may write:

[TABLE]

where $\lambda_{\mathbf{C}}=\theta_{2}{\bf{1}}^{T}\mathbf{D}\mathbf{D}^{T}{\bf{1}}$ is the only nonzero eigenvalue of $\mathbf{C}$. Now using lemma 2, we have:

[TABLE]

Taking advantage of (A-2), it is concluded that:

[TABLE]

Appendix B Proof of RIP Condition (Theorem 3)

Assume $\mathbf{D}^{\prime}=\mathbf{W}^{1/2}\mathbf{D}\in\mathbb{R}^{n\times p}$ satisfies the RIP condition with constant $\delta^{\prime}_{2k}$. If we denote the set of all $2k$-sparse vectors by $\Lambda_{2k}$ ($2k\leq p$), we have:

[TABLE]

where $\mathbf{B}=\mathbf{D}^{\prime T}\mathbf{D}^{\prime}=\mathbf{D}^{T}\mathbf{W}\mathbf{D}=\theta_{1}\mathbf{D}^{T}\mathbf{D}+\theta_{2}\mathbf{D}^{T}{\bf{1}}{\bf{1}}^{T}\mathbf{D}$. Also let $I_{2k}$ denote the set of indices corresponding with the non-zero elements of ${\bf{s}}$ and $\mathbf{B}_{2k}=\mathbf{D}_{2k}^{\prime T}\mathbf{D}^{\prime}_{2k}$ where $\mathbf{D}^{\prime}_{2k}=[{\bf{d}}^{\prime}_{i}],i\in I_{2k}$ and ${\bf{d}}^{\prime}_{i}$ denotes the $i$-th column of $\mathbf{D}^{\prime}$. Hence, incorporating the Min-Max theorem, it is sufficient to find $\delta^{\prime}_{2k}$ which satisfies:

[TABLE]

Now, similar to [52], using Gershgorin Circle Theorem (GCT) while assuming $||{\bf{d}}_{i}||=1,\,\forall i$, we may write:

[TABLE]

To avoid $\lambda_{\min}(\mathbf{B}_{2k})=0$, we assume $\mathbf{B}_{2k}$ is full rank, i.e., $n\geq 2k$. Thus, using (52), the upper bound for the maximum eigenvalue of $\mathbf{B}_{2k}$ is obtained as follows:

[TABLE]

where we have applied the Cauchy-Schwarz inequality for \big{|}\sum_{p}\sum_{q\neq p}d_{i,p}d_{j,q}\big{|}\leq(n-1)||{\bf{d}}_{i}||||{\bf{d}}_{j}||. Furthermore, for the minimum eigenvalue, we can write:

[TABLE]

Thus, it is sufficient to choose $\delta^{\prime}_{2k}$ such that:

[TABLE]

Assume $k_{2}=n-1$ and $\rho=k_{2}/k_{1}>1$. Substituting $\theta_{1}$ and $\theta_{2}$ from (6), we yield $\theta_{1}=1$ and $\theta_{2}<0$. Hence, we can simplify the condition above as follows:

[TABLE]

which leads to:

[TABLE]

where $C_{1}$ and $C_{2}$ are obtained by:

[TABLE]

If $n>1$ and $k>\frac{1}{2}(1+\frac{n}{n-1})>1$, we have $C_{1}>0$ and $C_{2}>1>\delta^{\prime}_{2k}$. Now to fulfil the assumption $\rho>1$, we must have $C_{1}>C_{2}-\delta^{\prime}_{2k}$ which implies:

[TABLE]

Since $h(n)$ is quadratic with respect to $n$, it has two roots specified by $n_{\min}$ and $n_{\max}$ where:

[TABLE]

The function \Delta(\mu(\mathbf{D}))=1-4\left(\frac{\delta^{\prime}_{2k}}{2k-1}-\mu(\mathbf{D})\right)\big{(}1-\mu(\mathbf{D})\big{)} itself, has two roots. The greater root satisfies $\mu_{\max}>1$ and the smaller one obtained by:

[TABLE]

satisfies $\mu_{\min}<\frac{\delta^{\prime}_{2k}}{2k-1}$. Now, since $\frac{\delta^{\prime}_{2k}}{2k-1}<1$ for $k>1$, we will have the following cases:

If $0<\mu(\mathbf{D})<\max\{0,\mu_{\min}\}$: Then $\Delta(\mu(\mathbf{D}))<0$ and for any $n\geq 2k$, $h(n)$ is positive. But for $\mu_{\min}$ to be positive, we need to have $k<\frac{4\delta^{\prime}_{2k}+1}{2}$ which only occurs for $k=2$ and $\delta^{\prime}_{2k}>\frac{3}{4}$. 2. 2.

$\max\{0,\mu_{\min}\}<\mu(\mathbf{D})<\frac{\delta^{\prime}_{2k}}{2k-1}$: In this case $0<\Delta(\mu(\mathbf{D}))<1$ and the function $h(n)$ is positive if $n>n_{\max}$ or $n<n_{\min}$. 3. 3.

$\frac{\delta^{\prime}_{2k}}{2k-1}<\mu(\mathbf{D})<1$: In this case $\Delta(\mu(\mathbf{D}))>1$ and hence $h(n)$ has two distinct roots. But only $n_{\min}$ is positive and thus inequality (59) holds if $n<n_{\min}$.

With our initial assumptions about the values of $k$ and $n$, the case $n<n_{\min}$ is not feasible. Hence, considering $\delta_{2k}<\sqrt{2}-1$, the only feasible case is when $0<\mu(\mathbf{D})<\frac{\delta^{\prime}_{2k}}{2k-1}$ and $n>n_{\max}$ where $k$ satisfies:

[TABLE]

Appendix C Proof of Convergence (Theorem 4)

We use similar procedure as in [26] for the proof of convergence of the proposed algorithm (Alg. LABEL:Algorithm_0) with fixed regularizing parameter $\alpha$. Before we proceed, let us define two notations used in this proof. We define $\delta{\bf{u}}^{(t+1)}={\bf{u}}^{(t+1)}-{\bf{u}}^{(t)}$ and $\Delta{\bf{u}}^{(t+1)}={\bf{u}}^{(t+1)}-{\bf{u}}^{\ast}$ where {\bf{u}}^{\ast}=\left({\bf{s}}^{\ast},{\bf{z}}^{\ast},{\mbox{\boldmath\mu}}^{\ast}_{1},{\mbox{\boldmath\mu}}^{\ast}_{2}\right)^{T} is the optimal solution to (23). Since CSIM is convex, the cost function is also convex. Hence, the KKT optimality conditions for problem (23) yields its global optimum:

[TABLE]

Consider equation (25) for the update of ${\bf{x}}$. Substituting {\mbox{\boldmath\mu}}^{(t)}_{1} and {\mbox{\boldmath\mu}}^{(t)}_{2} from (IV-D) and incorporating {\mbox{\boldmath\mu}}_{1}^{\ast}-\mathbf{M}^{T}{\mbox{\boldmath\mu}}_{2}^{\ast}={\bf{0}}, we get:

[TABLE]

Furthermore, optimizing problem (28) with respect to ${\bf{s}}$, while ${\bf{s}}_{0}={\bf{s}}^{(t)}$, we obtain the following:

[TABLE]

Now, according to lemma 3, we have:

[TABLE]

Using $\partial||{\bf{s}}^{\ast}||_{1}$ from (C) and $\partial||{\bf{s}}^{(t+1)}||_{1}$ from (65) and discarding $\alpha>0$, we can conclude that:

[TABLE]

Next, substituting \Delta{\mbox{\boldmath\mu}}_{1}^{(t+1)} from (64) we have:

[TABLE]

or equivalently:

[TABLE]

Now, using equation (32) for the update of ${\bf{z}}$ as well as the optimality condition 2(\mathbf{W}+\gamma\mathbf{I}){\bf{z}}^{\ast}=-{\mbox{\boldmath\mu}}^{\ast}_{2}, we yield:

[TABLE]

On the other side, incorporating the Lagrange multipliers update equations from (IV-D), we obtain:

[TABLE]

Hence, expanding (C), we have:

[TABLE]

Now, we define:

[TABLE]

where $\zeta>1$. Substituting the brackets in (C) from (64), while incorporating (70), we finally yield:

[TABLE]

Multiplying (C) by $2$, adding the term $-(\zeta-1)\sigma_{2}\Delta{{\bf{z}}^{(t+1)}}^{T}\delta{{\bf{z}}^{(t+1)}}^{T}$ to both sides of the inequality and finally using the Cauchy-Schwartz inequality $2{\bf{a}}^{T}{\bf{b}}\geq-\frac{1}{\beta}||{\bf{a}}||^{2}-{\beta}||{\bf{b}}||^{2}$ ($\beta>0$), we get the following:

[TABLE]

If $\beta_{3}$ and $\zeta$ are chosen such that $4(\mathbf{W}+\gamma\mathbf{I})-\beta_{3}(\zeta-1)\sigma_{2}\mathbf{I}$ is positive semi-definite, then we can discard the terms involving $\Delta{{\bf{z}}^{(t+1)}}$, in right hand side of (C). To fulfil this, we must have $4\lambda_{\min}(\mathbf{W})+4\gamma-\beta_{3}(\zeta-1)\sigma_{2}\geq 0$. Now, using (46) with $k_{2}=n-1$ and $\rho=k_{2}/k_{1}>1$, we conclude that:

[TABLE]

Next, using $\Delta{{\bf{u}}^{(t)}}=\Delta{{\bf{u}}^{(t+1)}}-\delta{{\bf{u}}^{(t+1)}}$ we proceed as follows:

[TABLE]

where the constant $\eta>0$ determines the rate of convergence of the algorithm. Combining (C) and (C) and discarding the terms with $\Delta{{\bf{z}}^{(t+1)}}$, we have:

[TABLE]

From $\eta>0$, we conclude that $\sigma_{1}<\beta_{1}<\frac{\lambda}{||\mathbf{D}||^{2}}\sigma_{1}$ and $\sigma_{2}<\beta_{2}<\sigma_{2}(\zeta-\frac{1}{\beta_{3}}(\zeta-1))$; which in turn implies $\lambda>\lambda_{\max}(\mathbf{D})=||\mathbf{D}||^{2}$ and $1<\beta_{3}$. Now, define $\delta_{1}=\frac{\lambda}{||\mathbf{D}||^{2}}-1$, $\beta_{1}=\sigma_{1}(1+r_{1}\delta_{1})$, $\delta_{2}=(\zeta-1)(1-1/\beta_{3})$, $\beta_{2}=\sigma_{2}(1+r_{2}\delta_{2})$, $\delta_{3}=\frac{\sigma_{2}}{4(\gamma+\frac{n-1}{\rho n})}$ and $\zeta=1+\frac{1}{r_{3}\delta_{3}}$. Subsisting these into (C), taking advantage of (76), it implies that we if choose $0<r_{1},r_{2}<1$ and $r_{3}>1$, then there will exist $\eta>0$ which satisfies (C) and is obtained by:

[TABLE]

Thus, using (C) to prove the convergence, we may write:

[TABLE]

Since $\sum_{t=0}^{\infty}||\delta{{\bf{u}}^{(t+1)}}||_{\mathbf{G}}^{2}>0$ is bounded from above, we can conclude that $||\delta{{\bf{u}}^{(t)}}||_{\mathbf{G}}\rightarrow 0$ or equivalently $\delta{\bf{u}}^{(t)}\rightarrow{\bf{0}}$ (since $\mathbf{G}$ is positive definite). This implies that the algorithm converges to a solution denoted by ${\tilde{{\bf{u}}}}=\lim_{t\rightarrow\infty}{\bf{u}}^{(t)}$. Now, since \delta{\mbox{\boldmath\mu}}^{(t)}_{1}\rightarrow{\bf{0}} and \delta{\mbox{\boldmath\mu}}^{(t)}_{2}\rightarrow{\bf{0}}, from (IV-D) we obtain:

[TABLE]

which implies $\tilde{{\bf{x}}}=\mathbf{D}\tilde{{\bf{s}}}$ and $\tilde{{\bf{z}}}=\mathbf{M}\tilde{{\bf{x}}}-{\bf{y}}$. Also using (32) and (64) we get:

[TABLE]

Hence, we can conclude that \tilde{{\mbox{\boldmath\mu}}}_{1}=-\mathbf{M}^{T}\tilde{{\mbox{\boldmath\mu}}}_{2} and 2(\mathbf{W}+\gamma\mathbf{I})\tilde{{\bf{z}}}=-\tilde{{\mbox{\boldmath\mu}}}_{2}. Together with (C) this implies that every limit point of $\{{\bf{u}}^{(t)}\}$ is an optimal solution to (23), i.e., ${\bf{u}}^{(t)}\rightarrow{\bf{u}}^{\ast}$.

It can also be easily verified that any solution to (23) is also a solution to (22).

Appendix D Lemmas and Proofs

Lemma 1**.**

Let $\mathbf{Q}=\mathbf{P}+\gamma\mathbf{I}_{n}$ where $\mathbf{P},\,\mathbf{Q}\in\mathbb{R}^{n\times n}$. If all the eigenvalues of $\mathbf{P}$ are given as $\lambda_{\mathbf{P}_{1}},\lambda_{\mathbf{P}_{2}},\ldots\lambda_{\mathbf{P}_{n}}$, then the eigenvalues of $\mathbf{Q}$ will be obtained using the following equation:

[TABLE]

Proof:

Let ${\bf{v}}_{\mathbf{P}_{i}}$ be an eigenvector of $\mathbf{P}$. We thus have:

[TABLE]

which implies that ${\bf{v}}_{\mathbf{Q}_{i}}={\bf{v}}_{\mathbf{P}_{i}}$ is an eigenvector and $\lambda_{\mathbf{Q}_{i}}=\lambda_{\mathbf{P}_{i}}+\gamma$ is an eigenvalue of $\mathbf{Q}$. ∎

Lemma 2**.**

For any matrices $\mathbf{A}$ and $\mathbf{B}$ with $\kappa$ defined in (13), we have:

[TABLE]

Proof:

Using the the definitions of $\sigma_{\min}(\mathbf{A})=\min_{{\bf{x}}\neq{\bf{0}}}\frac{||\mathbf{A}{\bf{x}}||}{||{\bf{x}}||}$ and $\sigma_{\max}(\mathbf{A})=\max_{{\bf{x}}\neq{\bf{0}}}\frac{||\mathbf{A}{\bf{x}}||}{||{\bf{x}}||}$ it is implied that $\sigma_{\min}(\mathbf{A})||{\bf{x}}||\leq||\mathbf{A}{\bf{x}}||\leq\sigma_{\max}(\mathbf{A})||{\bf{x}}||$. Consequently if we substitute ${\bf{x}}$ by $\mathbf{B}{\bf{y}}$, it is easy to show that $\sigma_{\min}(\mathbf{A})\sigma_{\min}(\mathbf{B})||{\bf{x}}||\leq||\mathbf{A}\mathbf{B}{\bf{x}}||\leq\sigma_{\max}(\mathbf{A})\sigma_{\max}(\mathbf{B})||{\bf{x}}||$ which concludes the proof. ∎

Lemma 3**.**

For any two vectors ${\bf{s}}^{\ast}$ and ${\bf{s}}^{(t+1)}$, we have:

[TABLE]

Proof:

Since $\ell_{1}$ norm is convex, exploiting the definition of sub-gradient of $f=||{\bf{s}}||_{1}$ at point ${\bf{s}}^{\ast}$, we obtain:

[TABLE]

Next, choose ${\bf{s}}^{(t+1)}$ as the starting point and ${\bf{s}}^{\ast}$ as the test point. Using a similar inequality, yields:

[TABLE]

Now, reversing (87) and adding it to (88), we get:

[TABLE]

∎