ICR: Iterative Convex Refinement for Sparse Signal Recovery Using Spike and Slab Priors

Abstract

In this letter, we address sparse signal recovery using spike and slab priors. In particular, we focus on a Bayesian framework where sparsity is enforced on reconstruction coefficients via probabilistic priors. The optimization resulting from spike and slab prior maximization is known to be a hard non-convex problem, and existing solutions involve simplifying assumptions and/or relaxations. We propose an approach called Iterative Convex Refinement (ICR) that aims to solve the aforementioned optimization problem directly allowing for greater generality in the sparse structure. Essentially, ICR solves a sequence of convex optimization problems such that sequence of solutions converges to a sub-optimal solution of the original hard optimization problem. We propose two versions of our algorithm: a.) an unconstrained version, and b.) with a non-negativity constraint on sparse coefficients, which may be required in some real-world problems. Experimental validation is performed on both synthetic data and for a real-world image recovery problem, which illustrates merits of ICR over state of the art alternatives.

I Introduction

SParse signal approximation and compressive sensing (CS) have recently gained considerable interest both in signal and image processing as well as statistics. Sparsity is often a natural assumption in inverse problems and sparse reconstruction or representation has variety of applications in image/signal classification [1, 2, 3, 4, 5, 6], dictionary learning [7, 8, 9, 10, 11, 12], signal recovery [13, 14], image denoising and inpainting [15], super resolution [16] and MRI image reconstruction [17]. Typically, sparse models assume that a signal can be efficiently represented as sparse linear combination of atoms in a given or learned dictionary [1, 18]. In other words, from CS viewpoint, a sparse signal can be recovered from fewer number of observations [19, 20, 21].

A typical sparse reconstruction algorithm aims to recover a sparse signal $\boldsymbol{x}\in\mathbb{R}^{p}$ from a set of fewer measurements $\boldsymbol{y}\in\mathbb{R}^{q}$ ($q\ll p$) according to the following model:

$$\boldsymbol{y}=\boldsymbol{A}\boldsymbol{x}+\boldsymbol{n},$$ (1)

where $\boldsymbol{A}\in\mathbb{R}^{q\times p}$ is the measurement matrix (Dictionary) and $\boldsymbol{n}\in\mathbb{R}^{q}$ models the additive Gaussian noise with variance $\sigma^{2}$.

In recent years, many sparse recovery algorithms have been proposed including but not limited to the following: proposing sparsity promoting optimization problems involving different regularizers such as $\ell_{1}$ norm, $\ell_{0}$ pseudo norm, greedy algorithms [14, 22, 23, 24], Bayesian-based methods [21, 25, 26] or general sparse approximation algorithms such as SpaRSA, ADMM, etc. [13, 27, 28, 29].

In this letter, we focus on sparse recovery from a Bayesian perspective by using hierarchical priors. In Bayesian sparse recovery, the choice of priors plays a key role in promoting sparsity and improving performance. Examples of such priors are Laplacian [30], generalized Pareto [31], Spike and Slab [32], etc. Amongst these priors, a well-suited sparsity promoting prior is spike and slab prior which is widely used in sparse recovery and Bayesian inference for variable selection and regression [33, 20, 17, 34]. In fact, it is acknowledged that spike and slab prior is indeed the gold standard for inducing sparsity in Bayesian inference [35].

Using these priors for sparse recovery leads to non-convex, non-smooth, mixed integer programming optimization problems which are often solved by means of relaxation or simplifying assumptions on the model parameters [36, 3, 17]. However, in this work we aim to solve the spike and slab optimization problem directly in its general form. Motivated by this, the Main Contributions of our work are as follows: (1) We propose a novel Iterative Convex Refinement (ICR) method to solve the optimization problem resulting from exploiting spike and slab priors. Essentially, the sequence of solutions from these convex problems approaches a sub-optimal solution of the hard non-convex problem. (2) We propose two versions of ICR: a.) an unconstrained version, and b.) with a non-negativity constraint on sparse coefficients, which may be required in some real-world problems such as image recovery. (3) Finally, we perform experimental validation on both synthetic data and a realistic image recovery problem, which reveals the benefits of ICR over other state-of-the-art recovery methods using spike and slab priors. Further, we compare the solution of various sparse recovery methods against the global solution for a small-scale problem, and remarkably the proposed ICR finds the most agreement with the global solution. Finally, convergence analysis is provided in support of the proposed ICR algorithm.

II Spike and Slab Sparse Signal Recovery

Introducing priors for capturing sparsity is a particular example of Bayesian inference where the signal recovery can be enhanced by exploiting contextual and prior information. As suggested by [37, 38], sparsity can be induced via solving the following optimization problem:

$$\max_{\boldsymbol{x}}P_{\boldsymbol{x}}(\boldsymbol{x}) ||\boldsymbol{y}-\boldsymbol{A}\boldsymbol{x}||_{2}<\epsilon.$$ subject to (2)

where $P_{\boldsymbol{x}}$ is the probability distribution function of $\boldsymbol{x}$ that captures sparsity. The most common example is the i.i.d. Laplacian prior which is equivalent to $\ell_{1}$ norm minimization [31]. In this work, we focus on the spike and slab prior for inducing sparsity on $\boldsymbol{x}$. Using this prior, every coefficient $x_{i}$ is modeled as a mixture of two densities as follows:

$$x_{i}\sim(1-\gamma_{i})\delta_{0}+\gamma_{i}P_{i}(x_{i})$$ (3)

where $\delta_{0}$ is the Dirac function at zero (spike) and $P_{i}$ (slab) is an appropriate prior distribution for nonzero values of $x_{i}$ (e.g. Gaussian). $\gamma_{i}\in[0,1]$ controls the structural sparsity of the signal. If $\gamma_{i}$ is chosen to be close to zero $x_{i}$ tends to remain zero. On the contrary, by choosing $\gamma_{i}$ close to 1, $P_{i}$ will be the dominant distribution encouraging $x_{i}$ to take a non-zero value.

Optimization Problem (Hierarchical Bayesian Framework): Inspired by Bayesian compressive sensing (CS) [21, 34], we employ a hierarchical Bayesian framework for signal recovery. In this model, priors are employed on $\boldsymbol{y}$ and $\boldsymbol{x}$. We also define $\gamma_{i}$ to be the indicator variable for the coefficient $x_{i}$, i.e. $\gamma_{i}\triangleq\mathbb{I}(x_{i}\neq 0)$. It takes the value one only if the corresponding coefficient $x_{i}$ is not zero, and zero otherwise. More precisely, the Bayesian formulation is as follows:

$$\boldsymbol{y}|\boldsymbol{A},\boldsymbol{x},\boldsymbol{\gamma},\sigma^{2} \sim \mathcal{N}\left(\boldsymbol{A}\boldsymbol{x},\sigma^{2}\boldsymbol{I}\right)$$ (4)

$$\boldsymbol{x}|\boldsymbol{\gamma},\lambda,\sigma^{2} \sim \prod_{i=1}^{p}\penalty\ \gamma_{i}\mathcal{N}(0,\sigma^{2}\lambda^{-1})+(1-\gamma_{i})\delta_{0}$$ (5)

$$\boldsymbol{\gamma}|\boldsymbol{\kappa} \sim \prod_{i=1}^{p}\penalty\ \mbox{Bernoulli}(\kappa_{i})$$ (6)

where $\mathcal{N}(.)$ represents the Gaussian distribution. Also note that in (5) each coefficient of $\boldsymbol{x}$ is modeled as i.i.d spike and slab prior. In addition, a Bernoulli distribution is used to model the indicator variable $\gamma_{i}$ with parameter $\kappa_{i}$, which controls the sparseness of the signal. Motivated by a recent maximum a posteriori (MAP) estimation technique proposed in [36] the optimal $\boldsymbol{x},\boldsymbol{\gamma}$ are obtained by the following MAP estimate.

$$(\boldsymbol{x}^{\ast},\boldsymbol{\gamma}^{\ast}) = \arg\max_{\boldsymbol{x},\boldsymbol{\gamma}}\left\{f(\boldsymbol{x},\boldsymbol{\gamma}|\boldsymbol{A},\boldsymbol{y},\boldsymbol{\kappa},\lambda,\sigma^{2})\right\}.$$ (7)

Proposition 1.

The MAP estimation above is equivalent to the following minimization problem:

$$(\boldsymbol{x}^{\ast},\boldsymbol{\gamma}^{\ast}) = \arg\min_{\boldsymbol{x},\boldsymbol{\gamma}}\penalty\ \penalty\ ||\boldsymbol{y}-\boldsymbol{A}\boldsymbol{x}||_{2}^{2}+\lambda||\boldsymbol{x}||_{2}^{2}+\sum_{i=1}^{p}\rho_{i}\gamma_{i}$$ (8)

where $\rho_{i}\triangleq\sigma^{2}\log\left(\frac{2\pi\sigma^{2}(1-\kappa_{i})^{2}}{\lambda\kappa_{i}^{2}}\right)$.

Proof.

See supplementary material. ∎

Remark: The optimization problem in (8) is a non-convex mixed integer programming involving the binary indicator variable $\boldsymbol{\gamma}$ and is not easily solvable using conventional optimization algorithms. It is worth mentioning that this is a more general formulation than the framework proposed in [3] or [36] where authors simplified the optimization problem by assuming the same $\kappa$ for each coefficient $x_{i}$. This assumption changes the last term in (8) to $\rho||\boldsymbol{x}||_{0}$ and the resulting optimization is solved in [36] by using Majorization-Minimization Methods. Further, a relaxation of $\ell_{0}$ to $\ell_{1}$ norm reduces the problem to the well-known Elastic-Net [39]. The framework in (8) therefore offers greater generality in capturing the sparsity of $\boldsymbol{x}$. As an example, consider the scenario in a reconstruction or classification problem where some dictionary (training) columns are more important than others[40]. It is then possible to encourage their contribution to the linear model by assigning higher values to the corresponding $\kappa_{i}$’s, which in turn makes it more likely that the $i^{th}$ coefficient $x_{i}$ becomes activated.

III Iterative Convex Refinement (ICR)

We first develop a solution to (8) for the case when the entries of $\boldsymbol{x}$ are non-negative. Then, we propose our method in its general form with no constraints.

The central idea of the proposed Iterative Convex Refinement (ICR) algorithm – see Algorithm 1 – is to generate a sequence of optimization problems that refines the solution of previous iteration based on solving a modified convex problem. At iteration $n$ of ICR, the indicator variable $\gamma_{i}$ is replaced with the normalized ratio $\frac{x_{i}}{\mu_{i}^{(n-1)}}$ and the convex optimization problem in (9) is solved which is a simple quadratic programming with non-negativity constraint. Note that, $\mu_{i}^{(n-1)}$ is intuitively the average value of optimal $x_{i}^{\ast}$’s obtained from iteration $1$ up to $n-1$ and is rigorously defined as in (11). The motivation for this substitution is that, if the sequence of solutions $\boldsymbol{x}^{(n)}$ converges to a point in $\mathbb{R}^{p}$ we also expect $\frac{x_{i}}{\mu_{i}^{(n-1)}}$ to converge to $\gamma_{i}$. Essentially, ICR is solving a sequence of convex quadratic programming problem that their solution converges to a sub-optimal solution of (8).

To generalize ICR to the unconstrained case, a simple modification is needed at each iteration. In fact, at each iteration (10) is solved instead of (9). Note that (10) is still convex and we solve it by alternating direction method of multipliers [29]. Again we expect the ratio $\frac{|x_{i}|}{|\mu_{i}^{(n-1)}|}$ to converge to the value of optimal $\gamma_{i}$ and the result of ICR be a sub-optimal solution for (8). ICR in both its versions is summarized in Algorithm 1.

Algorithm 1 Iterative Convex Refinement (ICR)

0:  $\boldsymbol{A},\boldsymbol{\kappa},\boldsymbol{y}$.initialize: $\boldsymbol{\mu}^{(0)}=\boldsymbol{A}^{T}\boldsymbol{y}$, iteration index $n=1$.

 while Stopping criterion not met do

  (1) Solve the convex optimization problem at iteration $n$:

  (Non-negative) For non-negative ICR solve

$$\boldsymbol{x}^{(n)}=\arg\min_{\boldsymbol{x}\succcurlyeq\boldsymbol{0}}||\boldsymbol{y}-\boldsymbol{A}\boldsymbol{x}||_{2}^{2}+\lambda||\boldsymbol{x}||_{2}^{2}+\sum_{i=1}^{p}\rho_{i}\frac{x_{i}}{\mu_{i}^{(n-1)}}$$ (9)

  (Unconstrained) For unconstrained ICR solve

$$\boldsymbol{x}^{(n)}=\arg\min_{\boldsymbol{x}}||\boldsymbol{y}-\boldsymbol{A}\boldsymbol{x}||_{2}^{2}+\lambda||\boldsymbol{x}||_{2}^{2}+\sum_{i=1}^{p}\rho_{i}\frac{|x_{i}|}{\big|\mu_{i}^{(n-1)}\big|}$$ (10)

  (2)

$$\text{Update\penalty\ }\mu_{i}^{(n)}:\penalty\ \penalty\ \penalty\ \mu_{i}^{(n)}=\frac{1}{n}\sum_{k=1}^{n}x_{i}^{(k)}\penalty\ \penalty\ \penalty\ i=1,...,p$$ (11)

  (3) Increase iteration index $n$.

 end while if $||\boldsymbol{x}^{(n)}-\boldsymbol{x}^{(n-1)}||\leq tol$

 $\boldsymbol{x}^{\ast}=\penalty\ \boldsymbol{x}^{(n-1)},\penalty\ \gamma_{i}^{\ast}=\frac{x_{i}^{\ast}}{\mu_{i}^{(n-1)}}$ for all $i=1,...,p$.

To analyze the convergence properties of ICR, we first define the function $f_{n}:\mathbb{R}^{p}\rightarrow\mathbb{R}$ as follows:

$$f_{n}(\boldsymbol{x})=\boldsymbol{x}^{T}(\boldsymbol{A}^{T}\boldsymbol{A}+\lambda\boldsymbol{I})\boldsymbol{x}-2\boldsymbol{y}^{T}\boldsymbol{A}\boldsymbol{x}+\sum_{i=1}^{p}\frac{\rho_{i}}{\big|\mu_{i}^{(n-1)}\big|}|x_{i}|$$ (12)

which is another form of the functions to be minimized at each iteration of ICR. With this definition and assuming $\alpha$ is a constant that $\alpha<\frac{1}{2(q+p)}$, we propose the following two lemmas with proofs in the supplementary material¹:

Lemma 1.

If —$\mu_{j}^{(n_{0})}\big|<\alpha\rho_{j}$, then $x_{j}^{(n_{0}+1)}=0$. ($\gamma_{j}\thickapprox\frac{x_{j}}{\mu_{j}^{(n_{0})}}=0$)

This lemma also implies that if $\big|\mu_{j}^{(n)}\big|<\alpha\rho_{j}$ for some $n$, then $x_{j}$ will remain zero for all the following iterations.

Lemma 2.

If $\big|\mu_{j}^{(n)}\big|\geq\alpha\rho_{j}$ for all $n\geq n_{0}$, then there exists $N_{j}\geq n_{0}$ such that for all $n>N_{j}$ we have

$$\bigg|\frac{1}{\big|\mu_{j}^{(n+1)}\big|}-\frac{1}{\big|\mu_{j}^{(n)}\big|}\bigg|\leq\frac{c}{n+1}$$ (13)

where $c$ is some positive constant.

Another interpretation of this lemma is that as the number of iterations grows, the cost functions at each iteration of ICR get closer to each other. In view of these two lemmas, we can show that the sequence of optimal cost function values obtained from ICR algorithm forms a Quasi-Cauchy sequence [41]. In other words, this is a sequence of bounded values that their difference at two consecutive iterations gets smaller.

Theorem 1.

After a sufficiently large $n$, the sequence of optimal cost function values obtained from ICR forms a Quasi-Cauchy sequence. i.e. $a_{n}=f_{n}(\boldsymbol{x}^{(n)})$ is a Quasi-Cauchy sequence of numbers.

$$\big|f_{n+1}(\boldsymbol{x}^{(n+1)})-f_{n}(\boldsymbol{x}^{(n)})\big|\leq\frac{c^{\prime}}{n}.$$ (14)

Proof.

We provide a sketch of the proof here, for more details please refer to the supplementary material¹.

Before proving the theorem, note that we can assume for a sufficiently large $N_{0}$, if $n\geq N_{0}$, then $\big|\mu_{j}^{(n)}\big|$ is either always less that $\alpha\rho_{j}$ or always greater (details can be found in the supplementary material). We now proceed to prove the Theorem and show that for $n>N_{0}$, the sequence of $f_{n}(\boldsymbol{x}^{(n)})$ satisfies the following property (assuming $|x_{i}|<1$):

$$\big|f_{n+1}(\boldsymbol{x}^{(n)})-f_{n}(\boldsymbol{x}^{(n)})\big|=\Big|\sum_{i=1}^{p}\rho_{i}\Big(\frac{1}{\big|\mu_{i}^{(n)}\big|}-\frac{1}{\big|\mu_{i}^{(n-1)}\big|}\Big)|x_{i}|\Big|$$ (15)

$$\leq \sum_{|\mu_{i}^{(n-1)}|<\alpha\rho_{i}}\rho_{i}\Bigg|\frac{1}{\big|\mu_{i}^{(n)}\big|}-\frac{1}{\big|\mu_{i}^{(n-1)}\big|}\Bigg||x_{i}|$$

$$\penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ +\sum_{|\mu_{i}^{(n-1)}|\geq\alpha\rho_{i}}\rho_{i}\Bigg|\frac{1}{\big|\mu_{i}^{(n)}\big|}-\frac{1}{\big|\mu_{i}^{(n-1)}\big|}\Bigg||x_{i}|$$

$$\leq \sum_{|\mu_{i}^{(n-1)}|\geq\alpha\rho_{i}}\rho_{i}\frac{c}{n}|x_{i}|\penalty\ \leq\penalty\ p\max{\{\rho_{i}\}}\frac{c}{n}\penalty\ \leq\penalty\ \frac{c^{\prime}}{n}.$$

This property also holds for $\boldsymbol{x}^{(n+1)}$. Finally, We show that for $n>N_{0}$, $a_{n}=f_{n}(\boldsymbol{x}^{(n)})$ is Quasi-Cauchy. Since the minimum value $f_{n+1}(\boldsymbol{x}^{(n+1)})$ is smaller than $f_{n+1}(\boldsymbol{x}^{(n)})$, we can write:

$$f_{n+1}(\boldsymbol{x}^{(n+1)})-f_{n}(\boldsymbol{x}^{(n)})\leq f_{n+1}(\boldsymbol{x}^{(n)})-f_{n}(\boldsymbol{x}^{(n)})\leq\frac{c^{\prime}}{n},$$

where we used (15) for $n>N_{0}$. With the same reasoning for $n>N_{0}$ we have:

$$f_{n+1}(\boldsymbol{x}^{(n+1)})-f_{n}(\boldsymbol{x}^{(n)})\geq f_{n+1}(\boldsymbol{x}^{(n+1)})-f_{n}(\boldsymbol{x}^{(n+1)})\geq-\frac{c^{\prime}}{n}.$$

Combining these two inequalities results (14) for $n>N_{0}$. ∎

Combination of this theorem with a reasonable stopping criterion guarantees the termination of the ICR algorithm. The stopping criteria used in this case is the norm of difference in the solutions $\boldsymbol{x}^{(n)}$ in consecutive iterations. At termination where the solution converges, the ratio $\frac{x_{i}}{\mu_{i}^{(n)}}$ will be zero for zero coefficients and approaches 1 for nonzero coefficients, which matches the value of $\gamma_{i}$ in both cases.

IV Experimental Validation

We now apply the ICR method to sparse signal recovery problem using spike and slab priors. Two experimental scenarios are considered: 1.) synthetic data and 2.) a real-world image recovery problem. In each case, comparisons are made against state of the art alternatives.

Synthetic data: We set up a typical experiment for sparse recovery as in [24, 36] with a randomly generated Gaussian matrix $\boldsymbol{A}\in\mathbb{R}^{q\times p}$ and a sparse vector $\boldsymbol{x}_{0}\in\mathbb{R}^{p}$. Based on $\boldsymbol{A}$ and $\boldsymbol{x}_{0}$, we form the observation vector $\boldsymbol{y}\in\mathbb{R}^{q}$ according to the additive noise model: $\boldsymbol{y}=\boldsymbol{A}\boldsymbol{x}_{0}+\boldsymbol{n}$ with $\sigma=0.01$. The competitive state-of-the-art methods for spike and slab sparse recovery that we compare against are: (1) SpaRSA [13, 42] which is a powerful method to solve the problems of the form (8). (2) Majorization Minimization (MM) algorithm [36] which aims to solve the spike and slab signal recovery problem through a majorization minimization approach. (3) Adaptive Elastic Net [39, 43] based on a $\ell_{1}$ relaxation of cost function. Initialization for all methods is consistent as suggested in [42].

Table I reports the experimental results for a small scale problem. We chose to first report results on a small scale problem in order to be able to use the IBM ILOG CPLEX optimizer [44] which is a very powerful optimization toolbox for solving many different optimization problems. It can also find the global solution to non-convex and mixed-integer programming problems. We used this feature of CPLEX to compare ICR’s solution with the global minimizer. For obtaining the results in Table I, we choose $p=64$, $q=32$ and the sparsity level of $\boldsymbol{x}_{0}$ is $10$. We generated $1000$ realizations of $\boldsymbol{A},\boldsymbol{x}_{0}$ and $\boldsymbol{n}$ and recovered $\boldsymbol{x}$ using different methods. Two different methods are used for evaluation of different sparse recovery methods: First, we compare different methods in terms of cost function value averaged over realizations, which is a direct measure of the quality of the solution to (8). Second, we compare performance of different methods from the sparse recovery point of view, and used the following figures of merit: mean square error (MSE) with respect to the global solution ($\boldsymbol{x}_{g}$) obtained by CPLEX optimizer, “Support Match” (SM) measure indicating how much the support of each solution matches to that of $\boldsymbol{x}_{g}$.

As can be seen from Table I, ICR outperforms the competing methods in many different aspects. In particular from the first row, we infer that ICR is a better solution to (8) since it achieves a better minimum in average sense. Moreover, significantly higher support match (SM $=97.25\%$ ) measure for ICR shows that ICR’s solution shows much more agreement with the global solution. Finally, the ICR solution is also the closest to the global solution obtained from CPLEX optimizer in the sense of MSE (by more than one order of magnitude in comparison with competing solutions).

Table I: Comparison of methods for $p=64$ and $q=32$. On Average, stopping criterion for ICR is achieved at iteration $n=11$.

Method	SpaRSA	MM	Elastic Net	ICR
Avg $f(\boldsymbol{x}^{*})$	1.8244E-2	1.4721E-2	3.1938E-2	1.3379E-2
MSE vs. $\boldsymbol{x}_{g}$	9.2668E-4	2.2290E-3	2.9210E-4	5.7851E-5
SM vs. $\boldsymbol{x}_{g}$ ($\%$)	84.46	83.12	67.46	97.25

Table II: Comparison of methods for $p=512$ and $q=128$. On Average, stopping criterion for ICR is achieved at iteration $n=15$.

Method	SpaRSA	MM	Elastic Net	ICR
Avg $f(\boldsymbol{x}^{*})$	4.7510E-2	3.9953E-2	4.7885E-2	3.9127E-2
MSE vs. $\boldsymbol{x}_{0}$	1.1740E-3	2.5685E-3	9.3378E-4	3.9277E-4
Sparsity Level	55.91	64.18	85.30	28.82
SM vs. $\boldsymbol{x}_{0}$ ($\%$)	89.68	82.41	87.61	96.33

Figure 1: Comparison of average MSE each method versus sparsity level of $\boldsymbol{x}_{0}$.

Next, we present results for a typical larger scale problem. We chose $p=512$, $q=128$ and set the sparsity level of $\boldsymbol{x}_{0}$ to be $30$ and carry out the same experiment as before. Because of the scale of the problem, the global solution is now unavailable and therefore, we compare the results against $\boldsymbol{x}_{0}$ which is the “ground truth”. Results are reported in Table II. Table II also additionally reports the average sparsity level of the solution and it can be seen that the sparsity level of ICR is the closest to the true sparsity level of $\boldsymbol{x}_{0}$. In all other figures of merit, viz. the cost function value (averaged over realizations), MSE and support match vs. $\boldsymbol{x}_{0}$, ICR is again the best. Fig. 1 shows an alternate result as the MSE plotted against the sparsity level; once again the merits of ICR are readily apparent.

Image reconstruction: In this part we aim to apply our ICR algorithm to real data for reconstruction of handwritten digit images from the well-known MNIST dataset [45]. The MNIST dataset contains 60000 digit images ($0$ to $9$) of size $28\times 28$ pixels. Most of pixels in these images are inactive and zero and only a few take non-zero values. Thus, these images are naturally sparse and fit into the spike and slab model. We set up this experiments such that a sparse signal $\boldsymbol{x}$ (vectorized image) is to be reconstructed from a smaller set of random measurements $\boldsymbol{y}$. For any particular image, we assume the smaller set of random measurement (150 measurements) is obtained by a Gaussian measurement matrix $\boldsymbol{A}\in\mathbb{R}^{150\times 784}$ with added noise according to (1). We compare our result against the following state-of-the-art image recovery methods for sparse images: 1.) SALSA-TV which uses the variable splitting proposed by Figueiredo et al. [46] combined with Total Variation (TV) regularizers [47]. 2.) A Bayesian Image Reconstruction (BIR) [25], based on a more recent version of Bayesian image reconstruction method [26] proposed by Hero et al.. We also compare our results with Adaptive Elastic Net method [39] which is commonly used in sparse image recovery problems. Finally, we also show the result of the non-negative version of ICR (ICR-NN) which explicitly enforces a non-negativity constraint on $\boldsymbol{x}$ which in this case corresponds to the intensity of reconstructed image pixels. Recovered images are shown in Fig. 2 and the corresponding average reconstruction error (MSE) for the whole database for different methods appears next to each method. Clearly, ICR and ICR-NN outperform the other methods both visually and based on MSE value. It is also intuitively satisfying that ICR-NN which captures the non-negativity constraint natural to this problem, provides the best result overall.

Figure 2: Examples of reconstructed images from MNIST dataset using different methods. The Numbers appeared next to each method is the average MSE for that method. On average, stoping criteria for ICR and ICR-NN are achieved at iteration $n=15$ and $n=29$, respectively.

V Conclusion

We develop a novel algorithm (ICR) for sparse recovery under spike and slab priors. Unlike known existing approaches, ICR does not simplify the optimization by assumptions/relaxations and hence affords a more general sparse structure. Experiments on synthetic data as well as a real-world image recovery problem confirms practical merits of ICR. Future research may investigate further analysis of ICR properties and extensions to multi-task sparse recovery under collaborative spike and slab priors.

[Experimental Validation] In this Appendix, we show more experimental results from our framework to further support its significance in comparison with other state-of-the-art methods for spike and slab sparse recovery problem. Following the same experimental setup for synthetic data as in the letter, we illustrate the performance of the ICR in comparison with others as the sparsity level of $\boldsymbol{x}_{0}$ ($||\boldsymbol{x}_{0}||_{0}$) changes. We vary the true sparsity level from only $5$ non-zero elements in $\boldsymbol{x}_{0}$ up to $95$ and compared MSE, support match percentage and the resulting sparsity level of the solutions from each method. Again we choose the length of sparse signal to be $p=512$ and number of observation to be $q=128$. Matrix $\boldsymbol{A}$ and sparse vector $\boldsymbol{x}_{0}$ are randomly generated and observation vector $\boldsymbol{y}$ is obtained by (1) with $\sigma=0.01$. $1000$ realization of $\boldsymbol{A}$, $\boldsymbol{x}_{0}$ and $\boldsymbol{n}$ are generated for each sparsity level and the results are averaged over these $1000$ realizations. Figures 1, 3 and 4 illustrate these results.

Fig. 3 illustrates that the support of ICR’s solution is the closest to the support of $\boldsymbol{x}_{0}$. More than $90\%$ match between the support of ICR’s solution and that of $\boldsymbol{x}_{0}$ for a wide range of sparsity levels makes ICR very valuable to variable selection problems specially in Bayesian framework. Fig. 4 shows the actual sparsity level of solution for different methods. The dashed line corresponds to the true level of sparsity and ICR’s solutions is the closest to the dashed line implying that the level of sparsity of ICR’s solution matches the level of sparsity of $\boldsymbol{x}_{0}$ more than other methods. This also support the results obtained from Fig. 3.

Figure 3: Comparison of average support match of the solution in $\%$ for each method versus sparsity level of $\boldsymbol{x}_{0}$.

Figure 4: Comparison of average sparsity level obtained by each method versus sparsity level of $\boldsymbol{x}_{0}$. Dashed line shows the true level of sparsity

[Analytical Results] In this appendix, we present the proofs to the theoretical lemmas and theorems in the paper. For the rest of our analysis, without loss of generality we assume that $|y_{i}|\leq 1$, $i=1...q$, $|x_{i}|\leq 1$, $i=1...p$ and columns of $\boldsymbol{A}$ have unity norm. We first begin with the proof to Proposition 1.

Proposition 1.

The MAP estimation in (7) is equivalent to the following minimization problem:

$$(\boldsymbol{x}^{\ast},\boldsymbol{\gamma}^{\ast}) = \arg\min_{\boldsymbol{x},\boldsymbol{\gamma}}\penalty\ \penalty\ ||\boldsymbol{y}-\boldsymbol{A}\boldsymbol{x}||_{2}^{2}+\lambda||\boldsymbol{x}||_{2}^{2}+\sum_{i=1}^{p}\rho_{i}\gamma_{i}$$ (16)

where $\rho_{i}\triangleq\sigma^{2}\log\left(\frac{2\pi\sigma^{2}(1-\kappa_{i})^{2}}{\lambda\kappa_{i}^{2}}\right)$.

Proof.

To perform the MAP estimation, note that the posterior probability is given by:

$$f(\boldsymbol{x},\boldsymbol{\gamma},|\boldsymbol{A},\boldsymbol{y},\lambda,\boldsymbol{\kappa})\propto f(\boldsymbol{y}|\boldsymbol{A},\boldsymbol{x},\boldsymbol{\gamma},\sigma^{2})f(\boldsymbol{x}|\boldsymbol{\gamma},\sigma^{2},\lambda)f(\boldsymbol{\gamma}|\boldsymbol{\kappa}).$$ (17)

The optimal $\boldsymbol{x}^{\ast},\boldsymbol{\gamma}^{\ast}$ are obtained by MAP estimation as:

$$(\boldsymbol{x}^{\ast},\boldsymbol{\gamma}^{\ast})=\arg\min_{\boldsymbol{x},\boldsymbol{\gamma}}\left\{-2\log f(\boldsymbol{x},\boldsymbol{\gamma},|\boldsymbol{A},\boldsymbol{y},\lambda,\boldsymbol{\kappa})\right\}.$$ (18)

We now separately evaluate each term on the right hand side of (17). According to (4) we have:

$$f(\boldsymbol{y}|\boldsymbol{A},\boldsymbol{x},\boldsymbol{\gamma},\sigma^{2})=\frac{1}{(2\pi\sigma^{2})^{q/2}}\exp\left\{-\frac{1}{2\sigma^{2}}(\boldsymbol{y}-\boldsymbol{A}\boldsymbol{x})^{T}(\boldsymbol{y}-\boldsymbol{A}\boldsymbol{x})\right\}$$

$$\Rightarrow-2\log f(\boldsymbol{y}|\boldsymbol{A},\boldsymbol{x},\boldsymbol{\gamma},\sigma^{2})=q\log\sigma^{2}+q\log(2\pi)+\frac{1}{\sigma^{2}}||\boldsymbol{y}-\boldsymbol{A}\boldsymbol{x}||^{2}.$$

Since $\gamma_{i}$ is assumed to be the indicator variable and only takes values $1$ and $0$, we can rewrite (5) in the following form:

$$\boldsymbol{x}|\boldsymbol{\gamma},\lambda,\sigma^{2} \sim \prod_{i=1}^{p}\penalty\ \Big(\mathcal{N}(0,\sigma^{2}\lambda^{-1})\Big)^{\gamma_{i}}.\Big(\delta_{0}\Big)^{1-\gamma_{i}}$$

Therefore

$$f\big(\boldsymbol{x}|\boldsymbol{\gamma},\sigma^{2},\lambda\big)=$$

$$\prod_{i=1}^{p}\left(\frac{1}{(2\pi\sigma^{2}/\lambda)^{1/2}}\right)^{\gamma_{i}}\exp\left(-\frac{\gamma_{i}x_{i}^{2}}{2\sigma^{2}\lambda^{-1}}\right)\delta_{0}^{1-\gamma_{i}}$$

$$= \left(\frac{2\pi\sigma^{2}}{\lambda}\right)^{-\frac{1}{2}\sum_{i=1}^{p}\gamma_{i}}\exp\left\{-\frac{1}{2\sigma^{2}\lambda^{-1}}\sum_{i=1}^{p}\gamma_{i}x_{i}^{2}\right\}\prod_{i=1}^{p}\delta_{0}^{1-\gamma_{i}}$$

$$= \left(\frac{2\pi\sigma^{2}}{\lambda}\right)^{-\frac{1}{2}\sum_{i=1}^{p}\gamma_{i}}\exp\left(-\frac{||\boldsymbol{x}||_{2}^{2}}{2\sigma^{2}\lambda^{-1}}\right)\prod_{i=1}^{p}\delta_{0}^{1-\gamma_{i}}$$

$$\Rightarrow -2\log f(\boldsymbol{x}|\boldsymbol{\gamma},\sigma^{2},\lambda)=$$

$$\frac{||\boldsymbol{x}||_{2}^{2}}{\sigma^{2}\lambda^{-1}}+\log\left(\frac{2\pi\sigma^{2}}{\lambda}\right)\sum_{i=1}^{p}\gamma_{i}-2\sum_{i=1}^{p}(1-\gamma_{i})\log\delta_{0}.$$

In fact $\delta_{0}=\mathbb{I}(x_{i}=0)$ and the final term on the right hand side evaluates to zero, since $\mathbb{I}(x_{i}=0)=1\Rightarrow\log\mathbb{I}(x_{i}=0)=0$, and $\mathbb{I}(x_{i}=0)=0\Rightarrow x_{i}\neq 0\Rightarrow\gamma_{i}=1\Rightarrow(1-\gamma_{i})=0$.

Finally (6) implies that

$$f(\boldsymbol{\gamma}|\boldsymbol{\kappa}) = \prod_{i=1}^{p}\kappa_{i}^{\gamma_{i}}(1-\kappa_{i})^{1-\gamma_{i}}$$

$$\Rightarrow-2\log f(\boldsymbol{\gamma}|\boldsymbol{\kappa}) = -2\sum_{i=1}^{p}\log\kappa_{i}^{\gamma_{i}}+\log(1-\kappa_{i})^{1-\gamma_{i}}$$

$$= -2\sum_{i=1}^{p}\gamma_{i}\log\kappa_{i}+(1-\gamma_{i})\log(1-\kappa_{i})$$

$$= -2\sum_{i=1}^{p}\gamma_{i}\log\Big(\frac{\kappa_{i}}{1-\kappa_{i}}\Big)+\log(1-\kappa_{i})$$

$$= \sum_{i=1}^{p}\gamma_{i}\log\left(\frac{1-\kappa_{i}}{\kappa_{i}}\right)^{2}-2\sum_{i=1}^{p}\log(1-\kappa_{i}).$$

Plugging all these expressions back into (18) and neglecting constant terms, we obtain:

$$(\boldsymbol{x}^{\ast},\boldsymbol{\gamma}^{\ast}) = \arg\min_{\boldsymbol{x},\boldsymbol{\gamma}}q\log\sigma^{2}+\frac{1}{\sigma^{2}}||\boldsymbol{y}-\boldsymbol{A}\boldsymbol{x}||^{2}+\frac{||\boldsymbol{x}||_{2}^{2}}{\sigma^{2}\lambda^{-1}}$$ (19)

$$+\log\left(\frac{2\pi\sigma^{2}}{\lambda}\right)\sum_{i=1}^{p}\gamma_{i}+\sum_{i=1}^{p}\gamma_{i}\log\left(\frac{1-\kappa_{i}}{\kappa_{i}}\right)^{2}$$

Essentially, for fixed $\sigma^{2}$ The cost function will reduce to:

$$L(\boldsymbol{x},\boldsymbol{\gamma}) = ||\boldsymbol{y}-\boldsymbol{A}\boldsymbol{x}||_{2}^{2}+\lambda||\boldsymbol{x}||_{2}^{2}+\sum_{i=1}^{p}\rho_{i}\gamma_{i}$$ (20)

where $\rho_{i}\triangleq\sigma^{2}\log\left(\frac{2\pi\sigma^{2}(1-\kappa_{i})^{2}}{\lambda\kappa_{i}^{2}}\right)$. ∎

Lemma 1.

If —$\mu_{j}^{(n_{0})}\big|<\alpha\rho_{j}$, then $x_{j}^{(n_{0}+1)}=0$. ($\gamma_{j}\thickapprox\frac{x_{j}}{\mu_{j}^{(n_{0})}}=0$)

Proof.

Assume that for a specific $j$, $\big|\mu_{j}^{(n_{0})}\big|<\alpha\rho_{j}$. Then for the next iteration the cost function to be minimized is as follows:

$$f_{n_{0}+1}(\boldsymbol{x})=\boldsymbol{x}^{T}(\boldsymbol{A}^{T}\boldsymbol{A}+\lambda\boldsymbol{I})\boldsymbol{x}-2\boldsymbol{y}^{T}\boldsymbol{A}\boldsymbol{x}+\sum_{i=1}^{p}\frac{\rho_{i}}{\big|\mu_{i}^{(n_{0})}\big|}|x_{i}|$$ (21)

Assume that the argument that minimizes (21) is $\boldsymbol{x}^{(n_{0}+1)}$. we can rewrite it in the following form:

$$\boldsymbol{x}^{(n_{0}+1)}=\boldsymbol{x}_{b}+x_{j}\boldsymbol{e}_{j}$$ (22)

where $\boldsymbol{e}_{j}$ is the $j^{th}$ basis function with one at component $j$ and zeros elsewhere. $x_{j}$ is the $j^{th}$ element of $\boldsymbol{x}^{(n_{0}+1)}$ and $\boldsymbol{x}_{b}$ is equal to $\boldsymbol{x}^{(n_{0}+1)}$ except at $j^{th}$ element which is zero. We prove that if $\big|\mu_{j}^{(n_{0})}\big|<\alpha\rho_{j}$, then $x_{j}=0$.

$$f_{n_{0}+1}(\boldsymbol{x}_{b}) = \boldsymbol{x}_{b}^{T}(\boldsymbol{A}^{T}\boldsymbol{A}+\lambda\boldsymbol{I})\boldsymbol{x}_{b}-2\boldsymbol{y}^{T}\boldsymbol{A}\boldsymbol{x}_{b}+\sum_{i=1}^{p}\frac{\rho_{i}}{\big|\mu_{i}^{(n_{0})}\big|}|x_{b_{i}}|$$

$$f_{n_{0}+1}(\boldsymbol{x}^{(n_{0}+1)}) = (\boldsymbol{x}_{b}+x_{j}\boldsymbol{e}_{j})^{T}(\boldsymbol{A}^{T}\boldsymbol{A}+\lambda\boldsymbol{I})(\boldsymbol{x}_{b}+x_{j}\boldsymbol{e}_{j})$$

$$-2\boldsymbol{y}^{T}\boldsymbol{A}(\boldsymbol{x}_{b}+x_{j}\boldsymbol{e}_{j})+\sum_{i=1}^{p}\frac{\rho_{i}}{\big|\mu_{i}^{(n_{0})}\big|}|x_{b_{i}}|$$

$$+\frac{\rho_{j}}{\big|\mu_{j}^{(n_{0})}\big|}|x_{j}|$$

Therefore, their difference is:

$$f_{n_{0}+1}\big(\boldsymbol{x}^{(n_{0}+1)}\big)-f_{n_{0}+1}\big(\boldsymbol{x}_{b}\big)\penalty\ =\penalty\ x_{j}^{2}\boldsymbol{e}_{j}^{T}(\boldsymbol{A}^{T}\boldsymbol{A}+\lambda\boldsymbol{I})\boldsymbol{e}_{j}$$

$$\penalty\ \penalty\ \penalty\ \penalty\ \penalty\ +2x_{j}\boldsymbol{x}_{b}^{T}(\boldsymbol{A}^{T}\boldsymbol{A}+\lambda\boldsymbol{I})\boldsymbol{e}_{j}-2x_{j}\boldsymbol{y}^{T}\boldsymbol{A}\boldsymbol{e}_{j}+\frac{\rho_{j}}{\big|\mu_{j}^{(n_{0})}\big|}|x_{j}|$$

$$=\big|x_{j}\big|\Big(|x_{j}|(\boldsymbol{A}^{T}\boldsymbol{A}+\lambda\boldsymbol{I})_{jj}+\frac{\rho_{j}}{\big|\mu_{j}^{(n_{0})}\big|}\Big)$$

$$-2x_{j}\Big(\boldsymbol{y}^{T}\boldsymbol{A}\boldsymbol{e}_{j}-\boldsymbol{x}_{b}^{T}(\boldsymbol{A}^{T}\boldsymbol{A}+\lambda\boldsymbol{I})\boldsymbol{e}_{j}\Big)$$ (23)

We want to show that this difference is always positive except for $x_{j}=0$ which means $x_{j}$ must be zero in order for $f_{n_{0}+1}\big(\boldsymbol{x}^{(n_{0}+1)}\big)$ to be minimum. To do so, we show the following statements are true for nonzero $x_{j}$:

$$\Big|2x_{j}\Big(\boldsymbol{y}^{T}\boldsymbol{A}\boldsymbol{e}_{j}-\boldsymbol{x}_{b}^{T}(\boldsymbol{A}^{T}\boldsymbol{A}+\lambda\boldsymbol{I})\boldsymbol{e}_{j}\Big)\Big|<\big|x_{j}\big|\Big(|x_{j}|(\boldsymbol{A}^{T}\boldsymbol{A}+\lambda\boldsymbol{I})_{jj}+\frac{\rho_{j}}{\big|\mu_{j}^{(n_{0})}\big|}\Big)$$

$$\Leftrightarrow 2\Big|\boldsymbol{y}^{T}\boldsymbol{A}\boldsymbol{e}_{j}-\boldsymbol{x}_{b}^{T}\boldsymbol{A}^{T}\boldsymbol{A}\boldsymbol{e}_{j}+\lambda\boldsymbol{x}_{b}^{T}\boldsymbol{e}_{j}\Big| < |x_{j}|(\boldsymbol{A}^{T}\boldsymbol{A}+\lambda\boldsymbol{I})_{jj}+\frac{\rho_{j}}{\big|\mu_{j}^{(n_{0})}\big|}$$

$$\Leftrightarrow\penalty\ \penalty\ \penalty\ 2\Big|\boldsymbol{y}^{T}\boldsymbol{A}\boldsymbol{e}_{j}-\boldsymbol{x}_{b}^{T}\boldsymbol{A}^{T}\boldsymbol{A}\boldsymbol{e}_{j}\Big| < (1+\lambda)|x_{j}|+\frac{\rho_{j}}{\big|\mu_{j}^{(n_{0})}\big|}$$

$$\Leftrightarrow\penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ 2\Big|(\boldsymbol{y}-\boldsymbol{A}\boldsymbol{x}_{b})^{T}\boldsymbol{A}\boldsymbol{e}_{j}\Big| < (1+\lambda)|x_{j}|+\frac{\rho_{j}}{\big|\mu_{j}^{(n_{0})}\big|}$$ (24)

In the above derivations, we used the fact that $(\boldsymbol{A}^{T}\boldsymbol{A})_{jj}=1$ since columns of $\boldsymbol{A}$ have unity norm. On the other hand, Cauchy-Schwarz inequality implies that,

$$2\Big|(\boldsymbol{y}-\boldsymbol{A}\boldsymbol{x}_{b})^{T}\boldsymbol{A}\boldsymbol{e}_{j}\Big| \leq\penalty\ \penalty\ 2||\boldsymbol{y}-\boldsymbol{A}\boldsymbol{x}_{b}||.||\boldsymbol{A}\boldsymbol{e}_{j}|| =\penalty\ \penalty\ 2||\boldsymbol{y}-\boldsymbol{A}\boldsymbol{x}_{b}||$$

$$\leq\penalty\ \penalty\ \penalty\ 2(||\boldsymbol{y}||+||\boldsymbol{A}\boldsymbol{x}_{b}||) \leq\penalty\ \penalty\ 2(\sqrt{q}+p)$$

Last inequality holds because of the fact that we assumed that magnitude of $x_{i}$ and $y_{i}$ do not exceed one. Also since we assumed $\big|\mu_{j}^{(n_{0})}\big|<\alpha\rho_{j}$ and by definition of $\alpha$ we have:

$$(1+\lambda)|x_{j}|+\frac{\rho_{j}}{\big|\mu_{j}^{(n_{0})}\big|}\geq\frac{1}{\alpha}\geq 2(q+p)\geq 2(\sqrt{q}+p)$$

Therefore, (24) is always true, since the right hand side is always greater than the left hand side. This implies that (23) is positive for nonzero $x_{j}$ and, hence we must have $x_{j}$ = 0. Otherwise, it would contradict the fact that $f\big(\boldsymbol{x}^{(n_{0}+1)}\big)$ is the minimum value. Note that these are loose bounds and in practice they are easily satisfied. For example, $||\boldsymbol{y}-\boldsymbol{A}\boldsymbol{x}_{b}||$ is practically very small. ∎

Lemma 2.

If $\big|\mu_{j}^{(n)}\big|\geq\alpha\rho_{j}$ for all $n\geq n_{0}$, then there exists $N_{j}\geq n_{0}$ such that for all $n>N_{j}$ we have

$$\bigg|\frac{1}{\big|\mu_{j}^{(n+1)}\big|}-\frac{1}{\big|\mu_{j}^{(n)}\big|}\bigg|\leq\frac{c}{n+1}$$ (25)

where $c$ is some positive constant.

Proof.

Assume $\big|\mu_{j}^{(n)}\big|\geq\alpha\rho_{j}=\epsilon$. First, note that it is straightforward to see that the difference of consecutive average values has the following property: $-\frac{1}{n+1}\leq\big|\mu_{j}^{(n+1)}\big|-\big|\mu_{j}^{(n)}\big|\leq\frac{1}{n+1}$. Now, let $N_{j}=\frac{2}{\alpha\rho_{j}}$, then for all $n>N_{j}$ we have:

$$\big|\mu_{j}^{(n)}\big|-\frac{1}{n+1}\leq\big|\mu_{j}^{(n+1)}\big|\leq\big|\mu_{j}^{(n)}\big|+\frac{1}{n+1}$$ (26)

where the left hand side is positive, since

$$\big|\mu_{j}^{(n)}\big|-\frac{1}{n+1}\geq\alpha\rho_{j}-\frac{1}{N_{j}}=\frac{\alpha\rho_{j}}{2}=\delta>0$$

Using this fact and (26) we infer that:

$$\frac{1}{\big|\mu_{j}^{(n)}\big|+\frac{1}{n+1}}\leq\frac{1}{\big|\mu_{j}^{(n+1)}\big|}\leq\frac{1}{\big|\mu_{j}^{(n)}\big|-\frac{1}{n+1}}$$

$$\Rightarrow \frac{1}{\big|\mu_{j}^{(n)}\big|+\frac{1}{n+1}}-\frac{1}{\big|\mu_{j}^{(n)}\big|}\leq\frac{1}{\big|\mu_{j}^{(n+1)}\big|}-\frac{1}{\big|\mu_{j}^{(n)}\big|}\leq$$

$$\penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \frac{1}{\big|\mu_{j}^{(n)}\big|-\frac{1}{n+1}}-\frac{1}{\big|\mu_{j}^{(n)}\big|}$$

$$\Rightarrow \frac{-\frac{1}{n+1}}{\Big(\big|\mu_{j}^{(n)}\big|+\frac{1}{n+1}\Big)\big|\mu_{j}^{(n)}\big|}\leq\frac{1}{\big|\mu_{j}^{(n+1)}\big|}-\frac{1}{\big|\mu_{j}^{(n)}\big|}\leq$$

$$\penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \frac{\frac{1}{n+1}}{\Big(\big|\mu_{j}^{(n)}\big|-\frac{1}{n+1}\Big)\big|\mu_{j}^{(n)}\big|}$$

In the last expression, we have:

$$\textit{RHS}=\frac{1}{(n+1)\Big(\big|\mu_{j}^{(n)}\big|-\frac{1}{n+1}\Big)\big|\mu_{j}^{(n)}\big|}\leq\frac{1}{(n+1)\epsilon\delta}$$

$$\textit{LHS}=\frac{-1}{(n+1)\Big(\big|\mu_{j}^{(n)}\big|+\frac{1}{n+1}\Big)\big|\mu_{j}^{(n)}\big|}\geq\frac{-1}{(n+1)\epsilon\delta}$$

Therefore,

$$\bigg|\frac{1}{\big|\mu_{j}^{(n+1)}\big|}-\frac{1}{\big|\mu_{j}^{(n)}\big|}\bigg|\leq\frac{1}{(n+1)\epsilon\delta},\penalty\ \penalty\ \penalty\ \penalty\ n>N_{j}.$$

∎

Theorem 1.

After a sufficiently large $n$, the sequence of optimal cost function values obtained from ICR forms a Quasi-Cauchy sequence. i.e. $a_{n}=f_{n}(\boldsymbol{x}^{(n)})$ is a Quasi-Cauchy sequence of numbers.

$$\big|f_{n+1}(\boldsymbol{x}^{(n+1)})-f_{n}(\boldsymbol{x}^{(n)})\big|\leq\frac{c^{\prime}}{n}$$ (27)

Proof.

Before proving the theorem, note that we can assume for a sufficiently large $N_{0}$, if $n\geq N_{0}$, then $\big|\mu_{j}^{(n)}\big|$ is either always less than $\alpha\rho_{j}$ or always greater. Because according to Lemma 1, we know that if $\big|\mu_{j}^{(n)}\big|$ once becomes smaller than $\alpha\rho_{j}$ for some $n$, it will remain less than $\alpha\rho_{j}$ for all the following iterations. Therefore, let $n_{j},\penalty\ j=1...p$ be the iteration index that for all $n>n_{j}$, $\big|\mu_{j}^{(n)}\big|<\epsilon$. Note that some $n_{j}$’s may be equal to infinity which means they are never smaller than $\epsilon$. For those $j$ that $n_{j}=\infty$, let $N_{j}$ to be the same as $N_{j}$ defined in proof of Lemma 2. With these definitions, we now proceed to prove the Theorem. We first show that for $n>N_{0}=\max(\max_{j}{n_{j}},\max_{j}{N_{j}})$, the sequence of $f_{n}(\boldsymbol{x}^{(n)})$ has the following property:

$$\big|f_{n+1}(\boldsymbol{x}^{(n)})-f_{n}(\boldsymbol{x}^{(n)})\big|=\Big|\sum_{i=1}^{p}\rho_{i}\Big(\frac{1}{\big|\mu_{i}^{(n)}\big|}-\frac{1}{\big|\mu_{i}^{(n-1)}\big|}\Big)|x_{i}|\Big|$$ (28)

$$\leq \sum_{|\mu_{i}^{(n-1)}|<\epsilon}\rho_{i}\Bigg|\frac{1}{\big|\mu_{i}^{(n)}\big|}-\frac{1}{\big|\mu_{i}^{(n-1)}\big|}\Bigg||x_{i}|$$

$$\penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ \penalty\ +\sum_{|\mu_{i}^{(n-1)}|\geq\epsilon}\rho_{i}\Bigg|\frac{1}{\big|\mu_{i}^{(n)}\big|}-\frac{1}{\big|\mu_{i}^{(n-1)}\big|}\Bigg||x_{i}|$$

$$\leq \sum_{|\mu_{i}^{(n-1)}|\geq\epsilon}\rho_{i}\frac{c}{n}|x_{i}|\penalty\ \leq\penalty\ p\max{\{\rho_{i}\}}\frac{c}{n}\penalty\ \leq\penalty\ \frac{c^{\prime}}{n}.$$

This property also holds for $\boldsymbol{x}^{(n+1)}$. Finally, We show that for $n>N_{0}$, $a_{n}=f_{n}(\boldsymbol{x}^{(n)})$ is Quasi-Cauchy. Since the minimum value $f_{n+1}(\boldsymbol{x}^{(n+1)})$ is smaller than $f_{n+1}(\boldsymbol{x}^{(n)})$, we can write:

$$f_{n+1}(\boldsymbol{x}^{(n+1)})-f_{n}(\boldsymbol{x}^{(n)})\leq f_{n+1}(\boldsymbol{x}^{(n)})-f_{n}(\boldsymbol{x}^{(n)})\leq\frac{c^{\prime}}{n}$$

where we used (28) for $n>N_{0}$. With the same reasoning for $n>N_{0}$ we have:

$$f_{n+1}(\boldsymbol{x}^{(n+1)})-f_{n}(\boldsymbol{x}^{(n)})\geq f_{n+1}(\boldsymbol{x}^{(n+1)})-f_{n}(\boldsymbol{x}^{(n+1)})\geq-\frac{c^{\prime}}{n}$$

Therefore,

$$\big|f_{n+1}(\boldsymbol{x}^{(n+1)})-f_{n}(\boldsymbol{x}^{(n)})\big|\leq\frac{c^{\prime}}{n}$$

for $n>N_{0}$. ∎

Remark: Despite the fact that analytical results show a decay of order $\frac{1}{n}$ in difference between consecutive optimal cost function values, ICR shows much faster convergence in practice.