Sparse Signal Recovery via Generalized Entropy Functions Minimization

TL;DR

This paper introduces generalized entropy functions as sparsity-promoting regularizers for compressive sensing, leading to improved sparse signal recovery by adaptively balancing high and low-energy coefficients.

Contribution

It proposes novel nonconvex entropy-based regularizers and an efficient optimization framework for enhanced sparse signal recovery.

Findings

Outperforms existing methods in simulated data

Achieves state-of-the-art results on real data

Efficiently solves the nonconvex optimization problem

Abstract

Compressive sensing relies on the sparse prior imposed on the signal of interest to solve the ill-posed recovery problem in an under-determined linear system. The objective function used to enforce the sparse prior information should be both effective and easily optimizable. Motivated by the entropy concept from information theory, in this paper we propose the generalized Shannon entropy function and R\'{e}nyi entropy function of the signal as the sparsity promoting regularizers. Both entropy functions are nonconvex, non-separable. Their local minimums only occur on the boundaries of the orthants in the Euclidean space. Compared to other popular objective functions, minimizing the generalized entropy functions adaptively promotes multiple high-energy coefficients while suppressing the rest low-energy coefficients. The corresponding optimization problems can be recasted into a series of reweighted $l_1$-norm minimization problems and then solved efficiently by adapting the FISTA. Sparse signal recovery experiments on both the simulated and real data show the proposed entropy functions minimization approaches perform better than other popular approaches and achieve state-of-the-art performances.

Tables

Table 1. TABLE I: Prior entropy-based functionals.

$g_6\left(𝒙\right)$[22]	$g_7\left(𝒙\right)$[23, 24, 22, 25]	$g_8\left(𝒙\right)$[26]	$g_9\left(𝒙\right)$[22]
$-{\sum }_i\frac{\left|x_i\right|}{{\Vert 𝒙\Vert }_1}\log \frac{\left|x_i\right|}{{\Vert 𝒙\Vert }_1}$	$-{\sum }_i\frac{{\left|x_i\right|}^2}{{\Vert 𝒙\Vert }_2^2}\log \frac{{\left|x_i\right|}^2}{{\Vert 𝒙\Vert }_2^2}$	$-{\sum }_i\left|x_i\right|\log \left|x_i\right|$	$-{\sum }_i{\left|x_i\right|}^2\log {\left|x_i\right|}^2$
$𝒙\in {ℝ}^N$, $𝒙\ne 0$	$𝒙\in {ℝ}^N$, $𝒙\ne 0$	$𝒙\in {ℝ}^N$, ${\Vert 𝒙\Vert }_1$ is fixed.	$𝒙\in {ℝ}^N$, ${\Vert 𝒙\Vert }_2^2$ is fixed.

Equations

$$\mathbf{y}=\mathbf{A}\mathbf{x}+\mathbf{w},$$

$$P_0^{ϵ}(\mathbf{x}):\underset{\mathbf{x}}{\min }\Vert \mathbf{x}{\Vert }_0\text{subject to }\Vert \mathbf{y}-\mathbf{A}\mathbf{x}{\Vert }_2^2\le ϵ.$$

$$P_0^S(\mathbf{x}):\underset{\mathbf{x}}{\min }\Vert \mathbf{y}-\mathbf{A}\mathbf{x}{\Vert }_2^2\text{subject to }\Vert \mathbf{x}{\Vert }_0\le S.$$

$$P_0^{\lambda }(\mathbf{x}):\underset{\mathbf{x}}{\min }\Vert \mathbf{y}-\mathbf{A}\mathbf{x}{\Vert }_2^2+\lambda \Vert \mathbf{x}{\Vert }_0,$$

$${\mathbf{x}}^{(t+1)}=\text{H}\left({\mathbf{x}}^{(t)}+{\mathbf{A}}^{\text{T}}(\mathbf{y}-\mathbf{A}{\mathbf{x}}^{(t)})\right),$$

$$P_g(\mathbf{x}):\underset{\mathbf{x}}{\min }\Vert \mathbf{y}-\mathbf{A}\mathbf{x}{\Vert }_2^2+\lambda g(\mathbf{x}),$$

$$h_p(\mathbf{x})=-\sum _{i=1}^N\frac{|x_i{|}^p}{\Vert \mathbf{x}{\Vert }_p^p}\log \frac{|x_i{|}^p}{\Vert \mathbf{x}{\Vert }_p^p},$$

$$h_{p,\alpha }(\mathbf{x})=\frac{1}{1-\alpha }\log \left(\sum _{i=1}^N{\left(\frac{|x_i{|}^p}{\Vert \mathbf{x}{\Vert }_p^p}\right)}^{\alpha }\right),$$

$$(1-{\delta }_S)\Vert \mathbf{x}{\Vert }_2^2\le \Vert \mathbf{A}\mathbf{x}{\Vert }_2^2\le (1+{\delta }_S)\Vert \mathbf{x}{\Vert }_2^2.$$

$$\mathcal{H}(\mathcal{V})=-{\sum }_{i=1}^{|\mathcal{V}|}\mathcal{P}(v_i)\log \mathcal{P}(v_i).$$

$${\mathcal{H}}_{\alpha }(\mathcal{V})=\frac{1}{1-\alpha }\log \left({\sum }_{i=1}^{|\mathcal{V}|}\mathcal{P}(v_i{)}^{\alpha }\right),$$

$$\mathbf{x}\to \left[\frac{|x_1{|}^p}{\Vert \mathbf{x}{\Vert }_p^p},\frac{|x_2{|}^p}{\Vert \mathbf{x}{\Vert }_p^p},\cdots ,\frac{|x_N{|}^p}{\Vert \mathbf{x}{\Vert }_p^p}\right],$$

$$P_{h_p}(\mathbf{x}):$$

$$P_{h_{p,\alpha }}(\mathbf{x}):$$

$$\frac{{\mathbf{x}}_{(1)}}{\Vert {\mathbf{x}}_{(1)}{\Vert }_p}$$

$$\frac{\mathbf{y}}{\Vert {\mathbf{x}}_{(1)}{\Vert }_p}=\frac{\mathbf{A}{\mathbf{x}}_{(1)}}{\Vert {\mathbf{x}}_{(1)}{\Vert }_p}$$

$$\text{sign}({\tilde{\mathbf{x}}}_{(1)})\circ |{\tilde{\mathbf{x}}}_{(1)}{|}^p=\ddot{\mathbf{x}}=\text{sign}({\tilde{\mathbf{x}}}_{(2)})\circ |{\tilde{\mathbf{x}}}_{(2)}{|}^p,$$

$${\tilde{\mathbf{x}}}_{(1)}=\text{sign}({\tilde{\mathbf{x}}}_{(1)})\circ |{\tilde{\mathbf{x}}}_{(1)}|=\text{sign}({\tilde{\mathbf{x}}}_{(2)})\circ |{\tilde{\mathbf{x}}}_{(2)}|={\tilde{\mathbf{x}}}_{(2)}.$$

$$h_p(\mathbf{x})=g(\ddot{\mathbf{x}})=-{\sum }_{i=1}^N|{\ddot{x}}_i|\log |{\ddot{x}}_i|$$

$$h_{p,\alpha }(\mathbf{x})=g_{\alpha }(\ddot{\mathbf{x}})=\frac{1}{1-\alpha }\log \left({\sum }_{i=1}^N|{\ddot{x}}_i{|}^{\alpha }\right).$$

$$\frac{{\mathbf{x}}_1}{\Vert {\mathbf{x}}_1{\Vert }_p}={\tilde{\mathbf{x}}}_1={\tilde{\mathbf{x}}}_i={\tilde{\mathbf{x}}}_2=\frac{{\mathbf{x}}_2}{\Vert {\mathbf{x}}_2{\Vert }_p}.$$

$${\mathbf{x}}_o⟷{\tilde{\mathbf{x}}}_o\text{and}\rho ({\mathbf{x}}_o)⟷\rho ({\tilde{\mathbf{x}}}_o).$$

$${\tilde{\mathbf{x}}}_o⟷{\ddot{\mathbf{x}}}_o\text{and}\rho ({\tilde{\mathbf{x}}}_o)⟷\rho ({\ddot{\mathbf{x}}}_o).$$

$$\begin{array}{rl}\mathbf{x}& =\mathbf{x}-\eta \cdot \partial h_p(\mathbf{x})\\ & =\mathbf{x}-\eta \cdot \text{sign}(\mathbf{x})\circ \nabla h_p(|\mathbf{x}|),\end{array}$$

$$\begin{array}{rl}& \frac{\partial h_p(\mathbf{x})}{\partial |x_i|}=\frac{\partial }{\partial |x_i|}\left(-\frac{{\sum }_{i=1}^N|x_i{|}^p\log |x_i{|}^p}{\Vert \mathbf{x}{\Vert }_p^p}+\log \Vert \mathbf{x}{\Vert }_p^p\right)\\ & =-\frac{1}{\Vert \mathbf{x}{\Vert }_p^p}\left(p|x_i{|}^{p-1}\log |x_i{|}^p+p|x_i{|}^{p-1}\right)\\ & +\frac{p|x_i{|}^{(p-1)}{\sum }_{l=1}^N|x_l{|}^p\log |x_l{|}^p}{\Vert \mathbf{x}{\Vert }_p^{2p}}+\frac{p|x_i{|}^{p-1}}{\Vert \mathbf{x}{\Vert }_p^p}\\ & =-\frac{p|x_i{|}^{(p-1)}\log |x_i{|}^p}{\Vert \mathbf{x}{\Vert }_p^p}+\frac{p|x_i{|}^{(p-1)}{\sum }_{l=1}^N|x_l{|}^p\log |x_l{|}^p}{\Vert \mathbf{x}{\Vert }_p^{2p}}.\end{array}$$

$$\displaystyle\frac{\partial h_{p}({\bm{x}})}{\partial|x_{i}|}\left\{\begin{array}[]{l}<0\\ =0\\ >0\end{array}\quad\begin{array}[]{l}\textnormal{if $|x_{i}|>\nu$}\\ \textnormal{if $|x_{i}|=\nu$}\\ \textnormal{if $|x_{i}|<\nu$}\,.\end{array}\right.$$

$$\begin{array}{rl}\frac{\partial h_{p,\alpha }(\mathbf{x})}{\partial |x_i|}& =\frac{1}{1-\alpha }\times \frac{1}{{\sum }_{l=1}^N{\left(\frac{|x_l|}{\Vert \mathbf{x}{\Vert }_p}\right)}^{p\alpha }}\times \frac{p\alpha }{\Vert \mathbf{x}{\Vert }_p^{p(\alpha +1)}}\\ & \times \left[|x_i{|}^{p\alpha -1}\Vert \mathbf{x}{\Vert }_p^p-|x_i{|}^{p-1}\Vert \mathbf{x}{\Vert }_{p\alpha }^{p\alpha }\right].\end{array}$$

$$\displaystyle\frac{\partial h_{p,\alpha}({\bm{x}})}{\partial|x_{i}|}\left\{\begin{array}[]{l}<0\\ =0\\ >0\end{array}\quad\begin{array}[]{l}\textnormal{if $|x_{i}|>\nu$}\\ \textnormal{if $|x_{i}|=\nu$}\\ \textnormal{if $|x_{i}|<\nu$}\,.\end{array}\right.$$

$$\sqrt{1-\delta (i)}\Vert {\mathbf{x}}_i{\Vert }_2\le \Vert {\mathbf{A}}_i{\mathbf{x}}_i{\Vert }_2\le \sqrt{1+\delta (i)}\Vert {\mathbf{x}}_i{\Vert }_2.$$

$${\sigma }_i=\Vert \mathbf{y}-{\mathbf{x}}_{\mathbf{I}}-{\mathbf{A}}_i{\mathbf{x}}_i{\Vert }_2,\forall i\in [1,\cdots ,C].$$

$$c=\arg {\min }_i{\sigma }_i.$$

$$\Vert \mathbf{y}-{\mathbf{x}}_{\mathbf{I}}-{\mathbf{A}}_1{\mathbf{x}}_1-{\mathbf{A}}_2{\mathbf{x}}_2-\overline{{\mathbf{A}}_{1,2}}\overline{{\mathbf{x}}_{1,2}}{\Vert }_2\le \sqrt{ϵ}.$$

Full text

I Introduction

Nowadays, there is an increasing amount of digital information generated constantly from every aspect of our life: the data that we work with grow in both size and variety. Fortunately, most of the data have inherent sparse structures. Compressive sensing [1, 2, 3] offers us an efficient framework to not only collect those data, but also process and analyze them in a timely fashion. Various compressive sensing tasks eventually boil down to the sparse signal recovery problem in an under-determined linear system:

[TABLE]

where ${\bm{y}}\in\mathbb{R}^{M}$ is the linear measurement, ${\bm{A}}\in\mathbb{R}^{M\times N}$ is the sensing matrix, ${\bm{x}}\in\mathbb{R}^{N}$ is the signal of interest with mostly zero entries, and ${\bm{w}}\in\mathbb{R}^{M}$ is the measurement noise.

In compressive sensing, we aim to recover the sparse signal ${\bm{x}}$ given $\{{\bm{y}},{\bm{A}}\}$, with $M\ll N$. In this case, the sensing matrix ${\bm{A}}$ contains more columns than rows, and there are more than one solutions satisfying the constraint $\|{\bm{y}}-{\bm{A}}{\bm{x}}\|_{2}^{2}\leq\epsilon$, where $\epsilon\geq 0$ is an upper bound on the measurement noise. This makes the recovery of ${\bm{x}}$ an ill-posed problem. On the other hand, since the signal of interest itself is sparse, the most straight-forward way to obtain a solution is to pick the one that also shares this “sparse” property. Following the well known Occam’s razor, we can simply use the $l_{0}$ norm as the criterion and choose the sparsest (simplest) one:

[TABLE]

This rather naïve attempt is actually backed up by the theoretical work of Candès et al. in [4, 5, 6, 7]. Under noiseless conditions, it can be shown that the sparsest solution is indeed the true signal when ${\bm{x}}$ is sufficiently sparse and ${\bm{A}}$ satisfies the corresponding restricted isometry property [4, 7].

Alternatively, we could impose an upper bound $S\in\mathbb{N}^{+}$ on the sparsity level of ${\bm{x}}$, and minimize the data fidelity term:

[TABLE]

Both $P_{0}^{\epsilon}({\bm{x}})$ and $P_{0}^{S}({\bm{x}})$ are nonconvex NP-hard problems in general [8, 9]. The $l_{0}$ norm can also be used as a regularizer [10], which gives rise to the following unconstrained problem:

[TABLE]

where $\lambda>0$ is the regularization parameter.

In practice, the following two approaches are usually employed to solve the sparse signal recovery problem:

Minimization of the $l_{0}$ norm $\|{\bm{x}}\|_{0}$.

• •

Various pursuit methods directly search for the nonzero entries of ${\bm{x}}$ incrementally or collectively.

Let $\widetilde{{\bm{x}}}$ denote the nonzero entries and $\widetilde{{\bm{A}}}$ denote the corresponding support set, matching pursuit (MP) [11] greedily searches for the sparsest representation of ${\bm{y}}$ by adding one column of ${\bm{A}}$ to $\widetilde{{\bm{A}}}$ at a time until the residue is no larger than $\epsilon$ or $\|{\bm{x}}\|_{0}$ reaches $S$. The corresponding coefficient is also added to $\widetilde{{\bm{x}}}$ one by one during the process. Orthogonal matching pursuit [12, 13] is an extension of MP, it updates all the coefficients in $\widetilde{{\bm{x}}}$ with the orthogonal projection of ${\bm{y}}$ onto $\widetilde{{\bm{A}}}$ in every iteration.

If the sparsity level $S$ is known, the support set $\widetilde{{\bm{A}}}$ can be further updated collectively during every iteration using CoSaMP [14] or subspace pursuit [15].

• •

In an effort to solve (3) and (4), Blumensath and Davis [10, 16, 17] proposed the iterative hard thresholding (IHT) to refine the solution ${\bm{x}}^{(t+1)}$ iteratively:

[TABLE]

where $\textnormal{H}({\bm{x}})$ is a nonlinear operator defined as follows:

1. a)

$\textnormal{H}_{S}({\bm{x}})$ preserves the largest $S$ entries of ${\bm{x}}$ in magnitude and sets the rest to [math]. 2. b)

$\textnormal{H}_{\lambda}({\bm{x}})$ preserves the entries whose magnitude is no less than $\lambda$ and sets the rest to [math].

The convergence of IHT is guaranteed provided that the operator norm $\|{\bm{A}}\|_{2}$ is less than $1$ [16]. IHT can be generalized to recover a row-sparse matrix ${\bm{X}}$ [18].

The various pursuit methods and IHT can be further modified under the oblique pursuits framework [19] to work with a wider class of measurement matrices. 2. 2.

Relaxation of the $l_{0}$ norm $\|{\bm{x}}\|_{0}$.

The relaxed alternative objective functions should be able to quantify the sparsity level of ${\bm{x}}$ and easy to optimize [20]. A detailed introduction is deferred to section I-A1.

Here we focus on studying the “relaxation” approach that solves the following unconstrained recovery problem:

[TABLE]

where $\lambda>0$ is the parameter that balances the trade-off between the data fidelity term $\|{\bm{y}}-{\bm{A}}{\bm{x}}\|_{2}^{2}$ and the sparsity regularizer $g({\bm{x}})$. The sparse prior is enforced via the regularizer $g({\bm{x}})$. A proper $g({\bm{x}})$ is crucial to the success of the sparse signal recovery task: it should favor sparse solutions, and the corresponding problem $P_{g}({\bm{x}})$ can be solved efficiently.

Built upon our earlier work in [21], in this paper we propose the following generalized Shannon entropy function and Rényi entropy function as new sparsity-promoting regularizers:

Shannon entropy function:

[TABLE]

where $p>0$. 2. 2.

Rényi entropy function:

[TABLE]

where $p>0$, $\alpha>0$ and $\alpha\neq 1$.

Both $h_{p}({\bm{x}})$ and $h_{p,\alpha}({\bm{x}})$ are nonconvex functions, scale-independent, and their local minimums occur at the boundary of each orthant in the Euclidean space $\mathbb{R}^{N}$, i.e. the axes. They share similar structures in $\mathbb{R}^{N}$. Take the Shannon entropy function for example, the 2D level plots of $h_{p}({\bm{x}})$ with $p=\{1,\,2\}$ are shown in Fig. 1. In this paper we show their advantages over other popular regularizers in promoting sparse solutions with both theoretical analyses and experimental evaluations.

I-A Prior Work

I-A1 Sparsity-promoting regularizers

Various regularizers have been proposed as the relaxations of the $l_{0}$ norm.

1. a)

Most popular among them are the $l_{p}$-norm based functions, i.e. the convex $l_{1}$ norm $\|{\bm{x}}\|_{1}$ [24, 26, 27] and the nonconvex $p$-th power of the $l_{p}$ norm $\|{\bm{x}}\|_{p}^{p}$, $0<p<1$ [25, 28, 29, 30]:

• •

$g_{1}({\bm{x}})=\|{\bm{x}}\|_{1}=\sum_{i}|x_{i}|$, where ${\bm{x}}\in\mathbb{R}^{N}$.

• •

$g_{2}({\bm{x}})=\|{\bm{x}}\|_{p}^{p}=\sum_{i}|x_{i}|^{p}$, where ${\bm{x}}\in\mathbb{R}^{N}$, $p\in(0,1)$.

Both $\|{\bm{x}}\|_{1}$ and $\|{\bm{x}}\|_{p}^{p}$ are additive separable functions, and strict error bounds can be established for the recovered solutions from $l_{1}$-minimization [4, 5, 6, 7] and $l_{p}^{p}$-minimization problems [28, 31, 29]. Compared to $l_{1}$-minimization, $l_{p}^{p}$-minimization has a tighter error bound and better sparse recovery performances. 2. b)

When ${\bm{x}}$ lies in the probability simplex, i.e. $\|{\bm{x}}\|_{1}=1$ and $x_{i}\geq 0,\forall i$, minimizing $\|{\bm{x}}\|_{1}$ would not work in this case. [32] proposed a lower bound $g_{3}({\bm{x}})$ on $\|{\bm{x}}\|_{0}$, and showed that the optimization problem can be efficiently solved using convex programming:

• •

$g_{3}({\bm{x}})=\frac{1}{\|{\bm{x}}\|_{\infty}}=\frac{1}{\max_{j}x_{j}}\leq\|{\bm{x}}\|_{0}$,

where ${\bm{x}}>\bm{0}$, $\|{\bm{x}}\|_{1}=1$.

[33] later proposed the following sparsity regularizer111The sparsity measure proposed in [33] is actually $-g_{4}({\bm{x}})$. The sign is flipped because we are minimizing the regularizer here.:

• •

$g_{4}({\bm{x}})=\frac{\|{\bm{x}}\|_{1}}{N\cdot\|{\bm{x}}\|_{\infty}}-1=\frac{\|{\bm{x}}\|_{1}}{N\cdot\max_{j}|x_{j}|}-1$,

where ${\bm{x}}\in\mathbb{R}^{N}$, ${\bm{x}}\neq\bm{0}$.

When $\|{\bm{x}}\|_{1}=1$ and ${\bm{x}}>0$, minimizing $g_{4}({\bm{x}})$ is equivalent to minimizing $g_{3}({\bm{x}})$. The function $g_{3}({\bm{x}})$ could be considered to be a modified special case of $g_{4}({\bm{x}})$. For readability purposes, we shall refer to $g_{4}({\bm{x}})$ as the “$\textnormal{L}_{1}/\textnormal{L}_{\infty}$” function in this paper.

Although the $\textnormal{L}_{1}/\textnormal{L}_{\infty}$ function and the proposed $h_{p}({\bm{x}})$, $h_{p,\alpha}({\bm{x}})$ are all scale-independent (see Fig. 1 for a 2D example), they differ in the following aspects:

• –

The proposed entropy functions are continuously differentiable w.r.t. the modulus $|{\bm{x}}|$, while the $\textnormal{L}_{1}/\textnormal{L}_{\infty}$ function is not.

• –

The proposed entropy functions have adaptive energy-promoting behaviors, while the $\textnormal{L}_{1}/\textnormal{L}_{\infty}$ function has a rigid energy-promoting behavior. The detailed analysis is presented in section II-D and Appendix B. 3. c)

Another popular class of regularizers is constructed from the logarithm function $\log(\cdot)$ [24, 34, 35, 25]. The following “logarithm of energy” function was proposed in [24] as a measure of sparsity:

• •

$g_{5}({\bm{x}})=\sum_{i:|x_{i}|\neq 0}\log|x_{i}|^{2}=2\sum_{i:|x_{i}|\neq 0}\log|x_{i}|$, where ${\bm{x}}\in\mathbb{R}^{N}$, ${\bm{x}}\neq\bm{0}$.

Candès et al. [36] later proposed the reweighted $l_{1}$-minimization algorithm as a way to enhance the sparsity, which is essentially the iterative minimization of the first-order approximation of $g_{5}({\bm{x}})$: $\sum_{i}\log(|x_{i}|+\epsilon)$, where $\epsilon>0$ is some small positive constant. 4. d)

As listed in Table I, entropy-based functionals have also been widely used to promote sparsity in ${\bm{x}}$. Specifically, $g_{6}({\bm{x}})$ and $g_{7}({\bm{x}})$ are special cases of the proposed generalized Shannon entropy function $h_{p}({\bm{x}})$ in (7) when $p=1$ and $2$. There is some imprecision in [22]’s analysis on the local minimums of $g_{7}({\bm{x}})$, which were believed to “occur just shy of the boundaries defined by the coordinate axes”. In section II-C, we can prove that the local minimums actually occur exactly on the coordinate axes. $g_{8}({\bm{x}})$ only promotes sparsity if the signal ${\bm{x}}$ has fixed $l_{1}$ norm, while $g_{9}({\bm{x}})$ works for signals with fixed $l_{2}$ norm. 5. e)

The functions $g_{6}({\bm{x}})$, $\cdots$, $g_{9}({\bm{x}})$ are all nonconvex. [37] later proposed the following convex entropy function as an approximation to the $l_{1}$-norm $\|{\bm{x}}\|_{1}$:

• •

$g_{10}({\bm{x}})=\sum_{i}\left[(|x_{i}|+e^{-1})\log(|x_{i}|+e^{-1})\,+\,e^{-1}\right]$,

where ${\bm{x}}\in\mathbb{R}^{N}$.

This way $g_{10}({\bm{x}})$ maintains the strictly convex property of $\|{\bm{x}}\|_{1}$ and is continuously differentiable in $\mathbb{R}^{N}$. However, we should note that $g_{10}({\bm{x}})$ only produces concentrated but not truly sparse solutions.

I-A2 $l_{1}$-norm minimization algorithms

When $g({\bm{x}})=\|{\bm{x}}\|_{1}$, the convex $l_{1}$-norm minimization problem can be efficiently solved using various readily available algorithms such as GPSR [38], ISTA [39, 40, 41], TwIST[42], FISTA [43, 44], NESTA [45, 46], ADM[47], etc. Detailed reviews of popular $l_{1}$-norm minimization algorithms can be found in [48, 49].

In this paper the optimization problem (6) is recasted into a series of reweighted $l_{1}$ problems [36, 30]. We choose to adapt the FISTA to solve (6) for its easy implementation and fast convergence.

I-A3 Uniqueness of the solution

The restricted isometry property (RIP) of the $M\times N$ sensing matrix ${\bm{A}}$ was introduced in [4] to provide guarantees for exact recovery. Let $\delta_{S}$ be the smallest number such that:

[TABLE]

holds for all vectors ${\bm{x}}$ with $\|{\bm{x}}\|_{0}\leq S$, then the matrix ${\bm{A}}$ is said to satisfy RIP with $S$-restricted constant $\delta_{S}$. In the noiseless case [7], when $\delta_{2S}<\sqrt{2}-1$, minimizing the convex relaxation $\|{\bm{x}}\|_{1}$ is equivalent to minimizing the nonconvex $\|{\bm{x}}\|_{0}$, and has a unique $S$-sparse solution. In the noisy case [7], a strict upper bound on the $l_{2}$ error $\|\hat{{\bm{x}}}-{\bm{x}}\|$ of the recovered $\|{\bm{x}}\|_{1}$-minimization solution $\hat{{\bm{x}}}$ can be established.

It was shown in [7] that the RIP can be met with high probability for Gaussian random matrix with i.i.d entries and $M>CS\log(N/S)$ for some constant $C$. Chartrand et al. [28, 31] later gave the uniqueness conditions for $\|{\bm{x}}\|_{p}^{p}$-minimization in terms of the RIP of ${\bm{A}}$. Both $\|{\bm{x}}\|_{1}$ and $\|{\bm{x}}\|_{p}^{p}$ are additive separable functions in that they can be written as a sum of identical functions involving individual entries of ${\bm{x}}$. However, the proposed entropy functions are non-separable, and different approaches are needed to find the properties of ${\bm{A}}$ that could guarantee the uniqueness of the solution.

I-B Main Contributions

Compared to previously adapted entropy functions [23, 24, 22, 25] where only $l_{1}$ and $l_{2}^{2}$ normalized ${\bm{x}}$ were studied, the proposed Shannon entropy function $h_{p}({\bm{x}})$ and Rényi entropy function $h_{p,\alpha}({\bm{x}})$ in (7, 8) are more general. They are defined with respect to the probability distribution in (12) where $p>0$ can be any positive number. This gives us more freedom in constructing the proper entropy functions for the sparse signal recovery task. As evident from the experiments on both simulated and real data, the ability to choose suitable $p$ values enables us to fully exploit the sparsity-promoting power of the generalized entropy functions and to achieve better performances over other state-of-the-art regularized approaches.

Previous works [22, 25] focused on the study of the Schur concavity of the generalized entropy functions with respect to ${\bm{x}}$ where $p=1,\,2$. They believed that the Shannon entropy function produced truly sparse solutions only when $p=1$. Here we show that it’s actually the (Schur) concavity with respect to the distribution $\mathscr{P}(\cdot)$ that really matters. The Shannon entropy function with $p=1$ is not the only case where truly sparse solutions can be obtained. In fact, $\forall\,p>0$ and $0<\alpha\neq 1$, we can prove in section II-C that the local minimums of the generalized entropy functions only occur at the boundaries of the orthants in $\mathbb{R}^{N}$. Hence minimizing $h_{p}({\bm{x}})$ or $h_{p,\alpha}({\bm{x}})$ in said orthant $\mathbb{O}$ will lead us to the solutions on its boundary, i.e. sparser solutions.

Additionally, minimizing the generalized entropy functions adaptively promotes multiple high-energy coefficients while suppressing the rest low-energy coefficients. Combining the proximal approximation [50, 51] of the data fidelity term and the first order approximations of the generalized entropy functions, we can recast the nonconvex entropy functions minimization problems into a series of simple reweighted $l_{1}$ problems, and solve them efficiently via the accelerated inexact proximal gradient method [43, 44, 52].

II Sparsity-Regularization Entropy Functions

II-A Entropy in Information Theory

Here we first introduce the entropy concepts in information theory [53, 54]. Both the Shannon entropy and Rényi entropy are defined with respect to the probability distribution $\bm{\mathscr{P}}(\mathcal{V})$ of some random variable $\mathcal{V}$. We give the definitions in terms of discrete probability distribution222For continuous distributions, the sum $\sum$ in (10,11) should be replaced with integration $\int$. as follows:

• •

Shannon entropy333The “$\log$” in this paper is by default natural logarithm, i.e. base $e$:

[TABLE]

$\mathcal{H}(\mathcal{V})$ is strictly concave with respect to the probability distribution $\bm{\mathscr{P}}(\mathcal{V})=\{\mathscr{P}(v_{1}),\cdots,\mathscr{P}(v_{|\mathcal{V}|})\}$.

• •

Rényi entropy:

[TABLE]

where $\alpha\geq 0$ and $\alpha\neq 1$. When $\alpha\in(0,1)$, $\mathcal{H}_{\alpha}(\mathcal{V})$ is strictly concave with respect to $\bm{\mathscr{P}}(\mathcal{V})$ [55]; when $\alpha\in(1,\infty)$, $\mathcal{H}_{\alpha}(\mathcal{V})$ is strictly Schur concave with respect to $\bm{\mathscr{P}}(\mathcal{V})$ [56].

We should make it clear that Shannon entropy $\mathcal{H}(\mathcal{V})$ is not a special case of the Rényi entropy, but the limiting value of the Rényi entropy $\mathcal{H}_{\alpha}(\mathcal{V})$ as $\alpha\rightarrow 1$ [57]. Hence we will discuss them separately in this paper.

II-B Generalized Entropy Functions of the Sparse Signal

Entropy measures the uncertainty about the random variable $\mathcal{V}$ with cardinality $|\mathcal{V}|=N$. The lower the entropy is, the more predictable the variable $\mathcal{V}$ will be. This corresponds to a skewed distribution $\bm{\mathscr{P}}(\mathcal{V})$. The idea of a skewed distribution could translate naturally to a sparse probability vector $\bm{\mathscr{P}}_{\mathcal{V}}=\left[\mathscr{P}(v_{1}),\cdots,\mathscr{P}(v_{N})\right]^{\textnormal{T}}$ in the sense that only a few probability entries of $\bm{\mathscr{P}}_{\mathcal{V}}\in\mathbb{R}^{N}$ are significant. In other words, the entropy can be used as a measure of how sparse the probability vector $\bm{\mathscr{P}}_{\mathcal{V}}$ is. This observation motivates us to adapt the concept of entropy as a sparsity-measure for the general signal ${\bm{x}}$, and use it as a regularizer in the sparse signal recovery task.

As mentioned before, the entropy is defined with respect to some probability distribution $\bm{\mathscr{P}}(\cdot)$. Here we can construct the following discrete probability distribution out of the signal of interest ${\bm{x}}\in\mathbb{R}^{N}$:

[TABLE]

where $p>0$. The adaptation from the classical entropy to the generalized entropy function is then pretty straightforward: we simply plug (12) into (10) and (11) to obtain the generalized Shannon entropy function (7) and the generalized Rényi entropy function (8).

Here we should not confuse the “entropy function” of ${\bm{x}}$ with the “entropy” of ${\bm{x}}$. Take the Shannon entropy of ${\bm{x}}$ for example: $\mathcal{H}({\bm{x}})=-\int_{{\bm{x}}}\mathscr{P}({\bm{x}})\log\mathscr{P}({\bm{x}})\,d{\bm{x}}$. $\mathcal{H}({\bm{x}})$ measures the uncertainty of the random variable ${\bm{x}}$, it is different from $h_{p}({\bm{x}})$ in (7) that measures the sparsity of ${\bm{x}}$.

In order to promote sparsity in the recovered solutions, we are going to minimize the generalized entropy functions. The sparse signal recovery problems in the form of (6) using the Shannon entropy function (SEF) $h_{p}({\bm{x}})$ and the Rényi entropy function (REF) $h_{p,\alpha}({\bm{x}})$ as regularizers then become:

[TABLE]

II-C Sparsity-promoting Analysis

We next show that $h_{p}({\bm{x}})$ and $h_{p,\alpha}({\bm{x}})$ can be used as sparsity regularizers in the following sense: minimizing them in an orthant $\mathbb{O}$ of the Euclidean space $\mathbb{R}^{N}$ leads to solutions on the boundary of said orthant, i.e. sparser solutions.

• ❑ Noiseless recovery:

In this case we are minimizing $h_{p}({\bm{x}})$ or $h_{p,\alpha}({\bm{x}})$ subject to the constraint ${\bm{y}}={\bm{A}}{\bm{x}}$, ${\bm{y}}\neq 0$. We first show that there is a one to one mapping in each orthant between ${\bm{x}}=\left[x_{1},\cdots,x_{N}\right]^{\textrm{T}}$ and $\ddot{{\bm{x}}}=\left[\ddot{x}_{1},\cdots,\ddot{x}_{N}\right]^{\textrm{T}}$, where $\ddot{x}_{i}=\textrm{sign}(x_{i})\cdot\frac{|x_{i}|^{p}}{\|{\bm{x}}\|^{p}_{p}}$. This will be done in two steps: Proposition 1 and Proposition 2.

Proposition 1**.**

If ${\bm{x}}$ is the solution to ${\bm{y}}={\bm{A}}{\bm{x}}$, ${\bm{y}}\neq 0$, then there is a one to one mapping in each orthant between ${\bm{x}}$ and $\tilde{{\bm{x}}}=\frac{{\bm{x}}}{\|{\bm{x}}\|_{p}}$.

Proof.

We need to prove ${\bm{x}}\longleftrightarrow\tilde{{\bm{x}}}$:

• –

It is easy to verify that ${\bm{x}}\rightarrow\tilde{{\bm{x}}}$.

• –

Suppose that $\tilde{{\bm{x}}}$ can be mapped to two solutions of ${\bm{y}}={\bm{A}}{\bm{x}}$: ${\bm{x}}_{(1)},{\bm{x}}_{(2)}$ in the same orthant. We then have:

[TABLE]

which tells us $\frac{{\bm{y}}}{\|{\bm{x}}_{(1)}\|_{p}}=\frac{{\bm{y}}}{\|{\bm{x}}_{(2)}\|_{p}}$. Since ${\bm{y}}\neq 0$, we have $\|{\bm{x}}_{(1)}\|_{p}=\|{\bm{x}}_{(2)}\|_{p}$. Using (15), we get ${\bm{x}}_{(1)}={\bm{x}}_{(2)}$. Hence ${\bm{x}}\leftarrow\tilde{{\bm{x}}}$.

∎

Proposition 2**.**

There is a one to one mapping in each orthant between $\tilde{{\bm{x}}}$ and $\ddot{{\bm{x}}}$.

Proof.

We need to prove $\tilde{{\bm{x}}}\longleftrightarrow\ddot{{\bm{x}}}$:

• –

We can rewrite $\ddot{{\bm{x}}}$ in terms of $\tilde{{\bm{x}}}$: $\ddot{{\bm{x}}}=\textrm{sign}(\tilde{{\bm{x}}})\circ|\tilde{{\bm{x}}}|^{p}$, where $\textrm{sign}(\tilde{{\bm{x}}})$ is a vector of the signs of every entry in $\tilde{{\bm{x}}}$. Hence $\ddot{{\bm{x}}}\leftarrow\tilde{{\bm{x}}}$.

• –

Suppose that $\ddot{{\bm{x}}}$ can be mapped to two points $\tilde{{\bm{x}}}_{(1)},\tilde{{\bm{x}}}_{(2)}$ in the same orthant. We then have:

[TABLE]

where “$\circ$” is the element-wise Hadamard product. Since $\textrm{sign}(\tilde{{\bm{x}}}_{(1)})=\textrm{sign}(\tilde{{\bm{x}}}_{(2)})$, we can get $|\tilde{{\bm{x}}}_{(1)}|=|\tilde{{\bm{x}}}_{(2)}|$.

[TABLE]

Hence $\tilde{{\bm{x}}}_{(1)}=\tilde{{\bm{x}}}_{(2)}$, and $\ddot{{\bm{x}}}\rightarrow\tilde{{\bm{x}}}$.

∎

Combining Proposition 1 and Proposition 2, we have ${\bm{x}}\longleftrightarrow\ddot{{\bm{x}}}$, as shown in Fig. 2.

Let $\mathcal{X}=\left\{{\bm{x}}_{1},{\bm{x}}_{2},\cdots\right\}$ be the solutions of ${\bm{y}}={\bm{A}}{\bm{x}}$, ${\bm{y}}\neq 0$ in one of the orthants $\mathbb{O}$. Specifically, $\mathcal{X}=\mathcal{X}_{1}\cup\mathcal{X}_{2}$ and $\mathcal{X}_{1}\cap\mathcal{X}_{2}=\emptyset$, where $\mathcal{X}_{1}$ contains solutions on the boundary of $\mathbb{O}$ and $\mathcal{X}_{2}$ contains the rest solutions that are not on the boundary. The solution ${\bm{x}}$ is then mapped to $\ddot{{\bm{x}}}$ one by one, producing the corresponding sets $\ddot{\mathcal{X}}_{1},\ddot{\mathcal{X}}_{2}$ on the simplex $\|\ddot{x}\|_{1}=1$. We can verify that the solutions in $\mathcal{X}_{1}$ are always sparser than those in $\mathcal{X}_{2}$. Additionally, we have the following Proposition 3:

Proposition 3**.**

In the noiseless case where ${\bm{y}}={\bm{A}}{\bm{x}}$, ${\bm{y}}\neq 0$, for every solution ${\bm{x}}\in\mathcal{X}_{2}$, there is a solution ${\bm{x}}^{*}\in\mathcal{X}_{1}$ on the boundary of the orthant $\mathbb{O}$ such that $h_{p}({\bm{x}}^{*})<h_{p}({\bm{x}})$ and $h_{p,\alpha}({\bm{x}}^{*})<h_{p,\alpha}({\bm{x}})$.

Proof.

By definition we have:

[TABLE]

For the SEF $h_{p}({\bm{x}})$, we first study the local minimums on the simplex $\|\ddot{{\bm{x}}}\|_{1}=1$ in the orthant $\mathbb{O}$. $g(\ddot{{\bm{x}}})$ is strictly concave with respect to $\ddot{{\bm{x}}}$, and the local minimums of $g(\ddot{{\bm{x}}})$ are on the boundary of $\mathbb{O}$. Hence for every $\ddot{{\bm{x}}}\in\ddot{\mathcal{X}}_{2}$, there is a $\ddot{{\bm{x}}}^{*}\in\ddot{\mathcal{X}}_{1}$ such that $g(\ddot{{\bm{x}}}^{*})<g(\ddot{{\bm{x}}})$.

For the REF $h_{p,\alpha}({\bm{x}})$, when $\alpha\in(0,1)$, $g_{\alpha}(\ddot{{\bm{x}}})$ is strictly concave with respect to $\ddot{{\bm{x}}}$, the local minimums of $g_{\alpha}(\ddot{{\bm{x}}})$ are on the boundary of the orthant $\mathbb{O}$. When $\alpha\in(1,\infty)$, $g_{\alpha}(\ddot{{\bm{x}}})$ is strictly Schur concave [56], since the boundary of the orthant $\mathbb{O}$ majorizes the $\ddot{{\bm{x}}}$ inside $\mathbb{O}$. The local minimums of $g_{\alpha}(\ddot{{\bm{x}}})$ are also on the boundary of $\mathbb{O}$. Hence for every $\ddot{{\bm{x}}}\in\ddot{\mathcal{X}}_{2}$, there also exists a $\ddot{{\bm{x}}}^{*}\in\ddot{\mathcal{X}}_{1}$ such that $g_{\alpha}(\ddot{{\bm{x}}}^{*})<g_{\alpha}(\ddot{{\bm{x}}})$ for $\alpha\in(0,1)\cup(1,\infty)$.

Proposition 1 and 2 guarantee a one to one mapping in $\mathbb{O}$ between ${\bm{x}}$ and $\ddot{{\bm{x}}}$: ${\bm{x}}\longleftrightarrow\ddot{{\bm{x}}}$. Since $h_{p}({\bm{x}})=g(\ddot{{\bm{x}}})$ and $h_{p,\alpha}({\bm{x}})=g_{\alpha}(\ddot{{\bm{x}}})$, for every ${\bm{x}}\in\mathcal{X}_{2}$, there is an ${\bm{x}}^{*}\in\mathcal{X}_{1}$ such that $h_{p}({\bm{x}}^{*})<h_{p}({\bm{x}})$ and $h_{p,\alpha}({\bm{x}}^{*})<h_{p,\alpha}({\bm{x}})$. ∎

From Proposition 3 we can see that minimizing $h_{p}({\bm{x}})$ or $h_{p,\alpha}({\bm{x}})$ in the orthant $\mathbb{O}$ leads us to sparser solutions in $\mathcal{X}_{1}$.

• ❑ Noisy recovery:

We can show similarly that minimizing $h_{p}({\bm{x}})$ or $h_{p,\alpha}({\bm{x}})$ subject to the constraint $\|{\bm{y}}-{\bm{A}}{\bm{x}}\|_{2}^{2}\leq\epsilon$ in an orthant $\mathbb{O}$ of the Euclidean space $\in\mathbb{R}^{N}$ also produces sparse solutions. First, we have the following Proposition 4:

Proposition 4**.**

Let $\mathcal{X}^{\epsilon}=\left\{{\bm{x}}_{1},{\bm{x}}_{2},\cdots\right\}$ are the nonzero solutions satisfying the constraint $\|{\bm{y}}-{\bm{A}}{\bm{x}}\|_{2}^{2}\leq\epsilon,{\bm{y}}\neq 0$ such that: $\forall{\bm{x}}_{i}\neq{\bm{x}}_{j}$, ${\bm{x}}_{i}=\tau{\bm{x}}_{j}$ for some $\tau>0$. Pick any ${\bm{x}}_{i}\in\mathcal{X}^{\epsilon}$, there is a one to one mapping in each orthant between the set $\mathcal{X}^{\epsilon}$ and $\tilde{{\bm{x}}}_{i}=\frac{{\bm{x}}_{i}}{\|{\bm{x}}_{i}\|_{p}}$.

Proof.

We need to prove $\mathcal{X}^{\epsilon}\longleftrightarrow\tilde{{\bm{x}}}_{i}$:

• –

$\forall{\bm{x}}_{j}\in\mathcal{X}^{\epsilon}\backslash{\bm{x}}_{i}$, we have $\tilde{{\bm{x}}}_{i}=\frac{{\bm{x}}_{i}}{\|{\bm{x}}_{i}\|_{p}}=\frac{\tau{\bm{x}}_{j}}{\|\tau{\bm{x}}_{j}\|_{p}}=\frac{{\bm{x}}_{j}}{\|{\bm{x}}_{j}\|_{p}}=\tilde{{\bm{x}}}_{j}$. We can verify that ${\bm{x}}_{j}\rightarrow\tilde{{\bm{x}}}_{j}=\tilde{{\bm{x}}}_{i}$. Hence $\mathcal{X}^{\epsilon}\rightarrow\tilde{{\bm{x}}}_{i}$.

• –

Suppose that there are two sets $\mathcal{X}^{\epsilon}_{(1)},\mathcal{X}^{\epsilon}_{(2)}$ in the same orthant being mapped to the same $\tilde{{\bm{x}}}_{i}$. Let ${\bm{x}}_{1}\in\mathcal{X}^{\epsilon}_{(1)}$ and ${\bm{x}}_{2}\in\mathcal{X}^{\epsilon}_{(2)}$, we have:

[TABLE]

We then have ${\bm{x}}_{1}=\frac{\|{\bm{x}}_{1}\|_{p}}{\|{\bm{x}}_{2}\|_{p}}{\bm{x}}_{2}$, which means that ${\bm{x}}_{1},{\bm{x}}_{2}$ belongs to the same set, i.e. $\mathcal{X}^{\epsilon}_{(1)}=\mathcal{X}^{\epsilon}_{(2)}$. Hence $\mathcal{X}\leftarrow\tilde{{\bm{x}}}_{i}$.

∎

We can see that the solutions in the same $\mathcal{X}^{\epsilon}$ all have the same entropy function value. Combining Proposition 4 and Proposition 2, we have $\mathcal{X}^{\epsilon}\longleftrightarrow\ddot{{\bm{x}}}_{i}$, as illustrated in Fig. 3.

Let $\mathcal{X}$ be the set defined in Proposition 4. Specifically, $\mathcal{X}^{\epsilon}=\mathcal{X}^{\epsilon}_{1}\cup\mathcal{X}^{\epsilon}_{2}$ and $\mathcal{X}^{\epsilon}_{1}\cap\mathcal{X}^{\epsilon}_{2}=\emptyset$, where $\mathcal{X}^{\epsilon}_{1}$ contains solutions on the boundary of the orthant $\mathbb{O}$, and $\mathcal{X}^{\epsilon}_{2}$ contains the rest solutions that are not on the boundary. The solution ${\bm{x}}$ is then mapped to $\ddot{{\bm{x}}}$ one by one, producing the corresponding sets $\ddot{\mathcal{X}}_{1},\ddot{\mathcal{X}}_{2}$ on the simplex $\|\ddot{x}\|_{1}=1$. Similarly we have:

Proposition 5**.**

In the noisy case where $\|{\bm{y}}-{\bm{A}}{\bm{x}}\|_{2}^{2}\leq\epsilon$, for every solution ${\bm{x}}\in\mathcal{X}^{\epsilon}_{2}$, there is a solution ${\bm{x}}^{*}\in\mathcal{X}^{\epsilon}_{1}$ on the boundary of the orthant $\mathbb{O}$ such that $h_{p}({\bm{x}}^{*})<h_{p}({\bm{x}})$ and $h_{p,\alpha}({\bm{x}}^{*})<h_{p,\alpha}({\bm{x}})$.

Proof.

From proposition 3, we have that for every $\ddot{{\bm{x}}}\in\ddot{\mathcal{X}}_{2}$, there is a $\ddot{{\bm{x}}}^{*}\in\ddot{\mathcal{X}}_{1}$ such that $g(\ddot{{\bm{x}}}^{*})<g(\ddot{{\bm{x}}})$.

Proposition 4 and 2 guarantee a one to one mapping in $\mathbb{O}$ between $\mathcal{X}^{\epsilon}$ and $\ddot{{\bm{x}}}$: $\mathcal{X}^{\epsilon}\longleftrightarrow\ddot{{\bm{x}}}$. Since $h_{p}({\bm{x}})=g(\ddot{{\bm{x}}})$ and $h_{p,\alpha}({\bm{x}})=g_{\alpha}(\ddot{{\bm{x}}})$, for every ${\bm{x}}\in\mathcal{X}^{\epsilon}_{2}$, there is an ${\bm{x}}^{*}\in\mathcal{X}^{\epsilon}_{1}$ such that $h_{p}({\bm{x}}^{*})<h_{p}({\bm{x}})$ and $h_{p,\alpha}({\bm{x}}^{*})<h_{p,\alpha}({\bm{x}})$. ∎

We can see that minimizing $h_{p}({\bm{x}})$ or $h_{p,\alpha}({\bm{x}})$ in the orthant $\mathbb{O}$ also leads to sparser solutions in $\mathcal{X}^{\epsilon}_{1}$ in the noisy case.

• ❑ Local minimums:

Here we formally introduce the following Proposition 6 on the local minimums of the generalized entropy functions:

Proposition 6**.**

The local minimums of the generalized entropy functions $h_{p}({\bm{x}})$ and $h_{p,\alpha}({\bm{x}})$ only occur at the boundaries of each orthant in $\mathbb{R}^{N}$.

Proof.

Without loss of generality, suppose ${\bm{x}}_{o}\in\mathbb{R}^{N}$ is a local minimum of the generalized entropy function inside some orthant $\mathbb{O}$, as illustrated in Fig. 4. There exists a small neighborhood $\rho({\bm{x}}_{o})$ surrounding ${\bm{x}}_{o}$ such that $\forall{\bm{x}}\in\rho({\bm{x}}_{o})$, $h_{p}({\bm{x}})\geq h_{p}({\bm{x}}_{o})$ or $h_{p,\alpha}({\bm{x}})\geq h_{p,\alpha}({\bm{x}}_{o})$.

As done in Proposition 1 and 4, we project ${\bm{x}}_{o}$ along with its neighborhood $\rho({\bm{x}}_{o})$ respectively to $\tilde{{\bm{x}}}_{o}$ and $\rho(\tilde{{\bm{x}}}_{o})$ on the sphere $\|\tilde{{\bm{x}}}\|_{p}^{p}=1$ inside the same orthant $\mathbb{O}$:

[TABLE]

$\tilde{{\bm{x}}}_{o}$ and $\rho(\tilde{{\bm{x}}}_{o})$ are further projected to $\ddot{{\bm{x}}}_{o}$ and $\rho(\ddot{{\bm{x}}}_{o})$ on the simplex $\|\ddot{{\bm{x}}}\|_{1}=1$ inside the same orthant $\mathbb{O}$.

[TABLE]

Consequently, $\forall\ddot{{\bm{x}}}\in\rho(\ddot{{\bm{x}}}_{o})$, we have $h_{p}(\ddot{{\bm{x}}})\geq h_{p}(\ddot{{\bm{x}}}_{o})$ or $h_{p,\alpha}(\ddot{{\bm{x}}})\geq h_{p,\alpha}(\ddot{{\bm{x}}}_{o})$. $\ddot{{\bm{x}}}_{o}$ is also a corresponding local minimum $\ddot{{\bm{x}}}_{o}$ on the simplex $\|\ddot{{\bm{x}}}\|_{1}=1$ inside the orthant $\mathbb{O}$.

However, we can show that such a local minimum $\ddot{{\bm{x}}}_{o}$ does not exist on the simplex $\|\ddot{{\bm{x}}}\|_{1}=1$ inside the orthant $\mathbb{O}$:

SEF $h_{p}(\ddot{{\bm{x}}})$ with $p>0$: It is strictly concave with respect to $\ddot{{\bm{x}}}$, there are no local minimums on the simplex $\|\ddot{{\bm{x}}}\|_{1}=1$ inside the orthant $\mathbb{O}$, i.e. $\ddot{{\bm{x}}}_{o}$ does not exist. 2. 2.

REF $h_{p,\alpha}(\ddot{{\bm{x}}})$ with $p>0$, $\alpha>0$ and $\alpha\neq 1$: When $\alpha\in(0,1)$, $h_{p,\alpha}(\ddot{{\bm{x}}})$ is strictly concave with respect to $\ddot{{\bm{x}}}$. When $\alpha\in(1,\infty)$, $h_{p,\alpha}(\ddot{{\bm{x}}})$ is strictly Schur concave with respect to $\ddot{{\bm{x}}}$ [56]. Hence $\ddot{{\bm{x}}}_{o}$ does not exist either.

We can see that there are no local minimums inside each orthant $\mathbb{O}$. Furthermore, the local minimums of $h_{p}(\ddot{{\bm{x}}})$ and $h_{p,\alpha}(\ddot{{\bm{x}}})$ occur at the boundaries of the simplex $\|\ddot{{\bm{x}}}\|_{1}=1$. Hence the local minimums of $h_{p}({\bm{x}})$ and $h_{p,\alpha}({\bm{x}})$ also only occur at the boundary of each orthant $\mathbb{O}$. ∎

II-D Energy-promoting Analysis

The sparsity-promoting property of the various regularizers already encourages the energy of the signal $\|{\bm{x}}\|_{2}^{2}$ to be concentrated in a few significant coefficients. Here we can show that minimizing the generalized entropy functions $h_{p}({\bm{x}})$ and $h_{p,\alpha}({\bm{x}})$ further promotes multiple high-energy coefficients while suppressing the rest low-energy coefficients within the recovered solution.

• ❑ SEF:

The energy-promoting property can be best illustrated from an optimization perspective. Take the SEF $h_{p}({\bm{x}})$ for example, should we use the subderivative descent to minimize it, the solution ${\bm{x}}$ is updated as follows:

[TABLE]

where $\eta>0$ is some suitable step size, and $\nabla h_{p}(|{\bm{x}}|)$ is the gradient with respect to the magnitude $|{\bm{x}}|$:

[TABLE]

It’s easy to verify the following remark:

Remark**.**

Let $\nu=\exp\left(\frac{\sum_{l}|x_{l}|^{p}\log|x_{l}|^{p}}{p\|{\bm{x}}\|_{p}^{p}}\right),\nu>0$, we have:

[TABLE]

As illustrated in Fig. 5, for the entries with relatively large magnitudes $|x_{i}|>\nu$, we can see that the update in (24) makes their magnitudes even larger. Conversely, if the entries have relatively small magnitudes $|x_{i}|<\nu$, the update makes their magnitudes even smaller. In this way minimizing $h_{p}({\bm{x}})$ promotes large-magnitude entries while suppressing small-magnitude entries in the recovered signal.

• ❑ REF:

For the REF $h_{p,\alpha}({\bm{x}})$, the derivative with respect to the magnitude $|{\bm{x}}|$ is:

[TABLE]

Similarly we have the following remark:

Remark**.**

Let $\nu=\exp\left(\frac{1}{p\alpha-p}\log\frac{\|{\bm{x}}\|_{p\alpha}^{p\alpha}}{\|{\bm{x}}\|_{p}^{p}}\right),\nu>0$, we have:

[TABLE]

We can see that minimizing $h_{p,\alpha}({\bm{x}})$ also promotes large-magnitude entries while suppressing small-magnitude entries in the recovered signal.

As discussed in Appendix B, we can see that minimizing the other regularizers introduced in section I-A1 do not have such energy-promoting behaviors.

II-D1 Advantage in SRC

This energy-promoting behavior becomes handy in applications such as sparse representation classification (SRC) [58, 59]. Specifically, suppose the subspace of interest is spanned by the columns of ${\bm{A}}=[{\bm{A}}_{1}\,{\bm{A}}_{2}\cdots{\bm{A}}_{C}\,{\bm{I}}]$, where ${\bm{A}}_{i}$ is a block matrix whose columns are the collected samples belonging to the $i$-th class, and ${\bm{I}}$ is the identity matrix used as the dictionary for the noise. The ${\bm{x}}$ can be written accordingly as ${\bm{x}}^{\textnormal{T}}=[{\bm{x}}_{1}^{\textnormal{T}}\,\,{\bm{x}}_{2}^{\textnormal{T}}\,\,\cdots{\bm{x}}_{C}^{\textnormal{T}}\,\,{\bm{x}}_{{\bm{I}}}^{\textnormal{T}}]$, where ${\bm{x}}_{i}$ contains the coefficients corresponding to the $i$-th class, and ${\bm{x}}_{{\bm{I}}}$ is the noise. Additionally, we assume the block matrix ${\bm{A}}_{i}$ satisfies the following restricted isometry property for some $\delta(i)\in(0,1)$ [7]:

[TABLE]

After the solution ${\bm{x}}$ that satisfies the constraint $\|{\bm{y}}-{\bm{A}}{\bm{x}}\|_{2}\leq\sqrt{\epsilon}$ is computed, the following residues are computed for every class:

[TABLE]

The class label of ${\bm{y}}$ is then:

[TABLE]

Without loss of generality, suppose ${\bm{y}}$ belongs to the $1$st class, and we are comparing $\sigma_{1}$ and $\sigma_{2}$. Let $\overline{{\bm{A}}_{1,2}}$ denote the matrix obtained after removing ${\bm{A}}_{1}$ and ${\bm{A}}_{2}$ from ${\bm{A}}$, while $\overline{{\bm{x}}_{1,2}}$ denotes the vector obtained after removing ${\bm{x}}_{1}$ and ${\bm{x}}_{2}$ from ${\bm{x}}$. We then have:

[TABLE]

Using the reverse triangle inequality, we have:

[TABLE]

Using the reverse triangle inequality, we have:

[TABLE]

We can see that (46) gives an upper bound on $\sigma_{1}$, and we would like to keep it small. Hence $\|{\bm{x}}_{i}\|_{2}$, $\forall i\neq 1$ should be small. While (48) gives a lower bound on $\sigma_{2}$, and we would like to keep it large. Hence $\|{\bm{x}}_{1}\|_{2}$ should be large, while $\|{\bm{x}}_{i}\|$, $\forall i\neq\{1,2\}$ should be small. In other words, we would like to promote the energy of ${\bm{x}}_{1}$ and suppress the energy of ${\bm{x}}_{i}$, $\forall i\neq 1$ at the same time. The energy-promoting property of the generalized entropy functions is well suited for this purpose.

III Entropy Function Minimization

Here we present the algorithm to perform the sparse signal recovery tasks formalized in (13),(14). The data fidelity term $f({\bm{x}})$ is convex and smooth, the sparsity regularizers $h_{p}({\bm{x}}),h_{p,\alpha}({\bm{x}})$ are nonconvex and nonsmooth. FISTA [43] is adapted as the accelerated inexact proximal gradient method [52] to minimize the nonconvex objective functions.

III-A Inexact Proximal Gradient Method

Using the generalized Shannon entropy function $h_{p}({\bm{x}})$ as an example, we first introduce the inexact proximal gradient method to solve the nonconvex problem (13):

• Step 1:

Proximal regularization using the data fidelity term $f({\bm{x}})=\|{\bm{y}}-{\bm{A}}{\bm{x}}\|_{2}^{2}$.

In the $(t+1)$-th iteration, we use the proximal regularization [50] of $f({\bm{x}})$ at ${\bm{x}}^{(t)}$, as done in [43]:

[TABLE]

where $o({\bm{x}}^{(t)})$ is a relative constant depending on the previous $t$-th iteration’s solution ${\bm{x}}^{(t)}$ and can thus be ignored; $\nabla f({\bm{x}}^{(t)})=2({\bm{A}}^{\textrm{T}}{\bm{A}}{\bm{x}}^{(t)}-{\bm{A}}^{\textrm{T}}{\bm{y}})$; the $\kappa$ is the Lipschitz constant of the gradient $\nabla f({\bm{x}}^{(t)})$ [60]. The smallest value $\kappa$ can take is twice the largest eigenvalue of ${\bm{A}}^{\textrm{T}}{\bm{A}}$, this is to ensure that $f({\bm{x}})$ in (49) is upper-bounded by the proximal regularization at the right-hand side. When $\kappa$ is unknown or difficult to compute, we can use the backtracking strategy describe in [43] to find it.

Plug (49) into (13), we need to solve the following proximal operator of the nonconvex function $h_{p}({\bm{x}})$:

[TABLE]

• Step 2:

The nonconvex proximal operator (50) is inexactly solved.

Since (50) is nonconvex, it is still an open-problem to obtain closed-form or exact solution [52, 61, 62]. Our goal here is to approximately solve it so that the solution ${\bm{r}}$:

[TABLE]

As shown in Appendix A, this would ensure the solution converges to a critical point.

The nonconvex function $h_{p}({\bm{x}})$ is replaced with its first-order approximation:

[TABLE]

where $\nabla h_{p}(|{\bm{x}}|)$ is the gradient with respect to the modulus $|x_{i}|$ given in (25). When computing $\nabla h_{p}(|x_{i}^{(t)}|)$, we add a small positive value $\epsilon=10^{-12}$ to $|x_{i}^{(t)}|$ to avoid $\log 0$. Ignoring the relative constant terms in (52) that depend only on ${\bm{x}}^{(t)}$, the proximal operator (50) can be approximated:

[TABLE]

where ${\bm{s}}^{(t)}={\bm{x}}^{(t)}-\frac{1}{\kappa}\nabla f({\bm{x}}^{(t)})$. This is a simple reweighted $l_{1}$-norm minimization problem that consists of a series of independent one-dimensional problems. The solution $r_{i}$ to (53) can be obtained using the soft shrinkage operator:

[TABLE]

where $\Gamma_{\tau}(\cdot)$ is defined as follows:

[TABLE]

Conventionally, the soft shrinkage operator solves a convex problem and requires the threshold $\tau$ to be positive. However, the gradient $\nabla h_{p}(|{\bm{x}}^{(t)}|)$ could be negative, which makes (53) nonconvex. In Appendix C we show that the global optimal solution can still be obtained using the soft thresholding operator in (59), yet with a different derivation process.

We further update ${\bm{r}}$ as follows to ensure (51) holds:

[TABLE]

III-B Accelerated Inexact Proximal Gradient Method

The proximal gradient method [50, 51] was proposed to solve convex optimization problems, and the error from the optimum objective function is guaranteed to be within $O(\frac{1}{t})$ after $t$ iterations. Beck and Teboulle [43, 44] later proposed the accelerated proximal gradient method, a.k.a. FISTA, that has a convergence rate of $O(\frac{1}{t^{2}})$. Li and Lin [52] later extended the accelerated proximal gradient (APG) to handle nonconvex and nonsmooth problems. In this case, the APG is guaranteed to converge to a critical point of the objective function [52]. Empirically, the APG could speed up the convergence by a great deal. However, it is still an open problem and a nontrival task to perform convergence rate analysis for general nonconvex objective functions [52, 61, 62].

Here we adapt the nonconvex acceleration approach by [52] to solve (13). Let $F_{p}({\bm{x}})=f({\bm{x}})+\lambda\,\,h_{p}({\bm{x}})$, we have:

[TABLE]

The iterations are stopped either when the maximum number of iterations is reached or when ${\bm{x}}^{(t+1)}$ reaches convergence, i.e. $\frac{\|{\bm{x}}^{(t+1)}-{\bm{x}}^{(t)}\|_{2}}{\|{\bm{x}}^{(t+1)}\|_{2}}\leq\tau_{c}$. The convergence threshold $\tau_{c}$ is usually chosen to be $0<\tau_{c}\leq 10^{-3}$.

Finally, the proposed entropy function minimization approach is summarized in Algorithm 1:

III-C Regularization Parameter

Naturally, choosing a proper regularization parameter $\lambda$ is the key to the success of sparse signal recovery.

III-C1 Noiseless recovery

The solution can be obtained by solving a series of minimization problems characterized by a decreasing $\lambda$ sequence: $\lambda\rightarrow 0$ [49]. Take the Shannon entropy function $h_{p}({\bm{x}})$ for example, we have the following:

1:Start with a relatively large $\lambda_{0}$, and set ${\bm{x}}^{\lambda_{0}}$ to the $l_{1}$-minimization solution.

2:for $l=\{1,\cdots,T\}$ do

3:   Use ${\bm{x}}^{\lambda_{l-1}}$ as the initializer, solve the minimization problem characterized by $\lambda_{l}$:

[TABLE]

4:   Update $\lambda_{l+1}=\rho\cdot\lambda_{l}$.

5:   if $\frac{\left\|{\bm{x}}^{\lambda_{l+1}}-{\bm{x}}^{\lambda_{l}}\right\|_{2}}{\|{\bm{x}}^{\lambda_{l+1}}\|_{2}}\leq\tau_{c}$ then break.

6:end for

To ensure the best performance, $\rho$ is chosen in $[0.9,1)$.

III-C2 Noisy recovery

The optimal $\lambda$ depends on the noise level and is usually unknown. An optimal fixed $\lambda$ can be pre-tuned on some development set and later used in practice.

III-D Constructing Suitable Entropy Functions

The choice of different $p$ and $\alpha$ values also greatly affects the recovery performance.

• •

For the SEF $h_{p}({\bm{x}})$, where $p>0$. The optimal $p$ value is usually around $1$. In practice, we can focus on tuning $p$ in the range $[0.8,\,1.3]$.

• •

For the REF $h_{p,\alpha}({\bm{x}})$, where $p>0$ and $0<\alpha\neq 1$. Both optimal $p$ and $\alpha$ values are usually around $1$. Empirically, we can focus on tuning $p$ in the range $[0.8,\,1.3]$, and $\alpha$ in the range $[0.9,\,1)\cup(1,\,1.1]$.

As an example, here we compare the noiseless recovery performances of different entropy functions. For the experiments, $100$ random sparse signals ${\bm{x}}$ of length $N=1000$ and sparsity level $S=200$ are generated: the nonzero entries are independently and identically drawn from $\mathcal{N}(0,1)$. The sensing matrices ${\bm{A}}$s are random Gaussian matrices with i.i.d entries following $\mathcal{N}(0,1)$, with the columns normalized and centralized. The number of measurements $M$ is varied in $[350,\,550]$, and Algorithm 1 is used to obtain the recovered signal $\hat{{\bm{x}}}$. The recovery is considered to be a success if $\frac{\|\hat{{\bm{x}}}-{\bm{x}}\|_{2}}{\|{\bm{x}}\|_{2}}<10^{-3}$. The success rates out of the $100$ recoveries for different choices of $M,p,\alpha$ are shown in Fig. 6.

III-E Initialization Strategy

A proper initialization is needed to ensure good performances for the nonconvex problems. Following the same experimental setup as section III-D, we choose $p=1.1$ for the SEF $h_{p}({\bm{x}})$, and $(p=1.1,\,\alpha=1.1)$ for the REF $h_{p,\alpha}({\bm{x}})$. We compare the following initialization strategies:

The $l_{1}$-minimization solution. 2. 2.

The $l_{2}$-minimization solution, i.e. the Moore-Penrose pseudo-inverse of ${\bm{y}}$: ${\bm{A}}^{\textnormal{T}}({\bm{A}}{\bm{A}}^{\textnormal{T}})^{-1}{\bm{y}}$. 3. 3.

Back projection: ${\bm{A}}^{\textnormal{T}}{\bm{y}}$. 4. 4.

All zeros: $[0,\cdots,0]^{\textnormal{T}}$.

The success rates using the above initializers are shown in Fig. 6. We can see that the solution from $l_{1}$-norm minimization usually leads to the best performances. It is thus used as the initializer in the proposed Algorithm 1.

IV Experimental Results

We compare the proposed Shannon entropy function (SEF) and Rényi entropy function (REF) regularizers with the relaxed sparsity-promoting regularizers introduced in section I-A1, OMP [12], CoSaMP [14] and the iterative hard thresholding (IHT) approach [10, 16, 17] that directly uses the $\|{\bm{x}}\|_{0}$ as a regularizer on both simulated and real datasets.

The accelerated inexact proximal gradient method in Section III-B can be adapted to solve (6) when $g({\bm{x}})$ is some other lower semi-continuous nonconvex/nonsmooth regularizers. The first-order approximation (52) needs to be updated accordingly. For best performances and a fair comparison, the $l_{1}$-minimization solution obtained via FISTA is used as the initializer for all the other approaches.

Take the regularizers introduced in section I-A1 for example, the function $g_{3}({\bm{x}})$ is a special case of $g_{4}({\bm{x}})$ for signals with fixed $l_{1}$ norm; the functions $g_{6}({\bm{x}})$, $\cdots$, $g_{9}({\bm{x}})$ are special cases of the proposed entropy functions; the function $g_{10}({\bm{x}})$ is not a sparsity-promoting regularizer. For the experiments, we thus compare the proposed entropy functions with the rest regularizers.

Specifically, The $l_{1}$ norm $g_{1}({\bm{x}})=\|{\bm{x}}\|_{1}$ is minimized using FISTA. The $l_{p}$ norm $g_{2}({\bm{x}})$ and the “$\textnormal{L}_{1}/\textnormal{L}_{\infty}$” function $g_{4}({\bm{x}})$ are minimized using Algorithm 1, i.e. adapted FISTA. We show that $g_{2}({\bm{x}})$, $g_{4}({\bm{x}})$ are both lower semi-continuous and give their first-order approximations in Appendix D. The “logarithm of energy” $g_{5}({\bm{x}})$ is not lower semi-continuous, its corresponding optimization problem is solved using the regularized FOCUSS algorithm [63, 64]: a reweighted $l_{2}$ norm minimization algorithm.

The detailed experimental settings and code files can be found at https://github.com/shuai-huang/GEFM.

IV-A Simulated Sparse Signal Recovery

IV-A1 Noiseless recovery

For the noiseless sparse signal recovery experiments, we fix $N=1000$ and vary the sampling ratio $\sigma=\frac{M}{N}\in[0.05,0.1,\cdots,0.95]$ and the sparsity ratio $\rho=\frac{S}{M}\in[0.05,0.1,\cdots,0.95]$, where $S$ is the sparsity level of the signal, i.e. the number of nonzero coefficients. As mentioned in section I-A3, generally the sensing matrix ${\bm{A}}$ needs to satisfy the RIP to guarantee the uniqueness of the recovered solution $\hat{{\bm{x}}}$. For each combination of $\sigma$ and $\rho$, we randomly generate $100$ pairs of $\{{\bm{x}},{\bm{A}}\}$: ${\bm{A}}$ is a $M\times N$ random Gaussian matrix with normalized and centralized columns. The nonzero entries of the sparse signal ${\bm{x}}\in\mathbb{R}^{N}$ are i.i.d. generated according to the Gaussian distribution $\mathcal{N}(0,1)$.

Given the measurement vector ${\bm{y}}={\bm{A}}{\bm{x}}$ and the sensing matrix ${\bm{A}}$, we try to recover the signal ${\bm{x}}$. If $\epsilon=\|\hat{{\bm{x}}}-{\bm{x}}\|_{2}/\|{\bm{x}}\|_{2}<10^{-3}$, the recovered $\hat{{\bm{x}}}$ is considered to be a success. A decreasing $\lambda$ sequence was used to approach the optimal $\lambda=0$ as explained in section III-C. The parameters were tuned to obtain best performance for each method: for the SEF minimization approach, $p=1.1$; for the REF minimization approach, $p=1.1,\,\alpha=1.1$; for the Lp minimization approach, $p=0.5$. The $\textnormal{L}_{1}/\textnormal{L}_{\infty}$ function has a rigid energy-promoting behavior. When it is minimized, only the highest-energy coefficient(s) is/are promoted while the rest coefficients are suppressed. In our experiments, we found that the recovered $\hat{{\bm{x}}}$ generally converged to its initializer, i.e. the $\textnormal{L}_{1}$-min solution. For the $\log$-energy minimization approach, as the lambda sequence went to [math], the regularized FOCUSS algorithm was equivalent to the original FOCUSS algorithm [64]. For the IHT approach, we scaled the sensing matrix ${\bm{A}}$ by its operator norm to make $\|{\bm{A}}\|_{2}<1$ as suggested by [10] to make the algorithm converge. For the OMP and CoSaMP, the sparsity level $S$ is tuned to obtain best performances. The measurement ${\bm{y}}$ was also scaled accordingly. Based on the $100$ trials, we computed the success recovery rate for each combination of $\{\sigma,\,\rho\}$, and plotted the phase transition curves (PTC) in Fig. 7.

As shown in Fig. 7, the PTC is the contour that corresponds to the 0.5 success rate in the domain $\{\sigma,\rho\}\in\{0,1\}^{2}$, it divides the domain into a “success” phase (lower right) and a “failure” phase (upper left). We can see that the proposed SEF minimization and REF minimization approaches generally perform equally well, and they both perform better than the other approaches.

IV-A2 Noisy recovery

We next try to recover the sparse signal ${\bm{x}}$ from a noisy measurement vector ${\bm{y}}$. Specifically, we fix $S=100,N=1000$ and increase $M$ gradually. The measurement ${\bm{y}}\in\mathbb{R}^{M}$ is generated as follows:

[TABLE]

where $\nu>0$ controls the amount of noise added to ${\bm{y}}$, the entries of ${\bm{w}}$ are i.i.d Gaussian $\mathcal{N}(0,1)$. We choose $\nu=0.05$, this creates a measurement ${\bm{y}}$ with signal to noise ratio (SNR) around $25$ dB. We randomly generate $100$ triples of $\{{\bm{x}},{\bm{A}},{\bm{w}}\}$. The average SNRs of the recovered signals $\hat{{\bm{x}}}$ are shown in Fig. 7. We can see that when $\sigma<0.4$, the SEF and REF minimization approaches outperform the other approaches.

IV-B Real Image Recovery

Real images are considered to be approximately sparse under some proper basis, such as the DCT basis, wavelet basis, etc. Here we compare the recovery performances of the aforementioned sparsity regularization approaches based on varying noiseless and noisy measurements of the $4$ real images shown in Fig. 9: Barbara, Boat, Lena, Peppers. Specifically, in order to reveal the sparse coefficients ${\bm{x}}$ of the real images ${\bm{s}}$, we use the sparsity averaging method by [65] to construct an over-complete wavelet basis by concatenating Db1-Db4 [66]:

[TABLE]

where $\Psi_{\textnormal{Db}i}$ is the corresponding operator of the Db$i$ wavelet. There is little performance improvement when adding more than $4$ wavelet operators. The $\frac{1}{2}$ is introduced in (72) to ensure ${\bm{s}}=\Psi({\bm{x}})$, and ${\bm{x}}=\Psi^{\dagger}({\bm{s}})$.

The sampling matrix ${\bm{U}}$ is constructed using the structurally random matrix approach by [67]: ${\bm{U}}={\bm{D}}{\bm{F}}{\bm{R}}$, where ${\bm{R}}$ is a uniform random permutation matrix that scrambles the signal’s sample locations globally while a diagonal matrix of Bernoulli random variables flips the signal’s sample signs locally, ${\bm{F}}$ is an orthonormal DCT matrix that computes fast transforms, ${\bm{D}}$ is a sub-sampling matrix that randomly selects a subset of the rows of the matrix ${\bm{F}}{\bm{R}}$.

The noiseless measurements ${\bm{y}}$ of the image ${\bm{s}}$ are obtained:

[TABLE]

The noisy measurements ${\bm{y}}$ are obtained:

[TABLE]

The entries of the noise ${\bm{w}}$ are generated using i.i.d. Gaussian distribution $\mathcal{N}(0,1)$, $\nu$ is chosen to be $0.02$ so that the SNR of the measurement vector ${\bm{y}}$ is around $30$ dB444When $\nu$ is set to other values, the relative performances of the four methods are similar..

Take the SEF minimization for example, the corresponding recovery problem is a cosparse problem [68] with respect to the image ${\bm{s}}$:

[TABLE]

The inexact proximal operator in the $(t+1)$-th iteration is:

[TABLE]

where $\kappa=2$ for the chosen sampling matrix ${\bm{U}}$. The inexact proximal operator (76) can be efficiently solved using the alternating split bregman shrinkage algorithm by [69].

Since the real images are only approximately sparse, both the noiseless and noisy recovery experiments are performed using a fixed $\lambda$. The parameters were tuned to obtain best performance for each approach555Since no direct matrix-vector product form is available for 72, OMP and CoSaMP are not compared in the image recovery experiments.. For the SEF minimization approach, $p=1$, $\lambda=5000$; for the REF minimization approach, $p=0.9$, $\alpha=1.1$, $\lambda=10^{4}$; for the $\textnormal{L}_{1}$ minimization approach, $\lambda=0.1$; for the Lp minimization approach, $p=0.8$, $\lambda=10^{-2}$; for the $\textnormal{L}_{1}/\textnormal{L}_{\infty}$ minimization approach, $\lambda=\frac{1}{10}N$; for the $\log$-energy minimization approach, $\lambda=1$; for the IHT approach, $\lambda=5$. The peak signal to noise ratios (PSNR) of the noiseless and noisy recovery experiments are shown in Fig. 8. We can see that the proposed SEF and REF minimization approaches perform equally well, and they generally have the best performances in terms of PSNR (dB). In the noisy recovery of “Lena” and “Peppers”, IHT performs slightly better than the entropy functions minimizations when $\sigma\geq 0.4$.

IV-C Face Recognition via Sparse Representation Classification

As introduced earlier in Section II-D1, SRC was proposed in [58, 59] as an effective approach to perform the classification task based on the sparse representation of a test sample ${\bm{y}}$ in the linear subspace spanned by the labeled training samples ${\bm{A}}$. The sparse representation ${\bm{x}}$ is computed by solving the following problem [58]:

[TABLE]

where $g({\bm{x}})$ is some sparsity-promoting regularizer. [58] used the $l_{1}$-norm $\|{\bm{x}}\|_{1}$ as the regularizer and tested the SRC model on the face recognition task. Significant and robust improvements over other popular approaches were observed across multiple benchmark datasets .

Following the same experimental setup as [58], we choose the challenging task of face recognition under random pixel corruption and compare the performances of the aforementioned sparsity regularization approaches. Specifically, each image in the Extended Yale B Face Database [70] is resized to $96\times 84$ pixels and then converted to a vector via column-concatenation. Each test image is randomly corrupted by replacing the pixel value with a number uniformly chosen from [0, 255], as show in Fig. 10.

In this case, the dictionary ${\bm{A}}$ is composed of two sub-dictionaries: ${\bm{A}}=[{\bm{A}}_{tr}\quad{\bm{I}}]$. ${\bm{A}}_{tr}$ is the collection of noiseless training samples that contributes to the “uncorrupted” test image, while ${\bm{I}}$ is the identity matrix that contributes to the random noise. ${\bm{x}}$ can then be separated into two parts accordingly: ${\bm{x}}^{T}=[{\bm{x}}_{tr}^{T}\quad{\bm{x}}_{{\bm{I}}}^{T}]$. Based on the observation that the respective distributions of the entries in ${\bm{x}}_{tr}$ and ${\bm{x}}_{{\bm{I}}}^{T}$ are quite different, we propose to regularize them using different regularization parameters as follows:

[TABLE]

As done before, the parameters were tuned for each approach separately to obtain best performances. For the SEF minimization approach, $p=1.2$, $\lambda=10^{-3}$, $\mu=75$; for the REF minimization approach, $p=1.15$, $\alpha=1.1$, $\lambda=10^{-3}$, $\mu=75$; for the $\textnormal{L}_{1}$ minimization approach, $\lambda=10^{-3}$, $\mu=2.5$; for the Lp minimization approach, $p=0.5$, $\lambda=10^{-4}$, $\mu=5$; for the $\textnormal{L}_{1}/\textnormal{L}_{\infty}$ minimization approach, $\lambda=10^{-6}$, $\lambda=2.5$; for the $\log$-energy minimization approach, $\lambda=10^{-6}$, $\mu=5$; for the IHT approach, $\lambda=10^{-3}$, $\mu=6$. For the OMP and CoSaMP, the sparsity constraint on the training samples is set to $S_{tr}\leq 20$, corresponding to the average number of samples per class in the experiments; the sparsity constraint on the noise is set to $S_{{\bm{I}}}\leq 7000$. The recognition rates of various approaches under increasing percentages of random pixel corruptions are shown in Fig. 11. We can see that the sparse representations obtained by minimizing the proposed SEF and REF lead to the highest recognition rates. The performance only starts to drop significantly when more than $90\%$ of the pixels are corrupted.

V Conclusion and Future Work

In this paper we proposed the generalized Shannon entropy function $h_{p}({\bm{x}})$ and Rényi entropy function $h_{p,\alpha}({\bm{x}})$ as the sparsity regularizers for the sparse signal recovery task. Regardless of the values of $p\in(0,\infty]$ and $\alpha\in(0,1)\cup(1,\infty]$, the local minimums of the generalized entropy functions occur on the boundaries of the orthants in $\mathbb{R}^{N}$, and minimizing them produces sparse solutions. Both $h_{p}({\bm{x}})$ and $h_{p,\alpha}({\bm{x}})$ are nonconvex functions, the corresponding minimization problems (13),(14) can be solved via the adapted FISTA until convergence. Compared to other sparsity regularizers, minimizing the generalized entropy functions adaptively promotes multiple high-energy coefficients while suppressing the rest low-energy coefficients of the recovered solution.

Sparse signal recovery experiments on both the simulated and real data show the proposed entropy function regularizations perform better than OMP, CoSaMP and other popular sparsity regularizers, and achieve state-of-the-art performances. This motivates us to explore sufficient conditions for the uniqueness of the solution, seek matrices that satisfy such uniqueness conditions, and establish error bounds on the recovered signal $\hat{{\bm{x}}}$ in our future work. Additionally, we would like to apply the generalized entropy functions minimization approaches to other Compressive Sensing applications such as RPCA [71, 72, 73, 74], dictionary learning[75, 76, 77], etc.

Appendix A Monotonic Objective Function Minimization

Take the Shannon entropy function minimization (SEF) problem $P_{h_{p}}({\bm{x}})$ in (13) for example, we have the following:

[TABLE]

The first inequality “(a)” is obtained using (49). The second inequality “(b)” holds because we enforce (51) when the nonconvex proximal operator (50) is inexactly solved. Combined with (69), the proposed Algorithm 1 is monotonic.

Additionally, we have:

• •

$h_{p}({\bm{x}})$ is continuous in $\mathbb{R}^{N}$, hence it is also lower semi-continuous.

• •

The data fidelity term $f({\bm{x}})=\|{\bm{y}}-{\bm{A}}{\bm{x}}\|_{2}^{2}\geq 0$ is bounded from below. $f({\bm{x}})\rightarrow\infty$ when $\|{\bm{x}}\|_{2}\rightarrow\infty$, hence $f({\bm{x}})$ is coercive. $h_{p}({\bm{x}})\geq 0$ is also bounded from below. Hence $F({\bm{x}})=f({\bm{x}})+\lambda h_{p}({\bm{x}})$ is coercive.

• •

$f({\bm{x}})$ has Lipschitz-continuous gradients [43].

Using Theorem $1$ of [52], we can get that the proposed Algorithm 1 converges to a critical point.

Appendix B Energy-promoting Analysis of other regularization functions

Here we discuss the energy-promoting behaviors of the other sparsity regularizers666$g_{10}({\bm{x}})$ is not a sparsity regularizer. introduced in section I-A1 from an optimization perspective. Specifically, the regularizers are minimized via subderivative descent:

[TABLE]

• •

The subderivatives of $g_{1}({\bm{x}})$, $g_{2}({\bm{x}})$, $g_{5}({\bm{x}})$ are as follows:

[TABLE]

The update (24) would decrease the magnitude of every coefficient.

• •

The subderivatives of $g_{3}({\bm{x}})$, $g_{4}({\bm{x}})$ are:

[TABLE]

where $\delta_{\infty}(i)$ is defined as follows:

[TABLE]

The update (24) would only increase the magnitude(s) of the highest-energy coefficient(s) for $g_{3}({\bm{x}})$ and $g_{4}({\bm{x}})$, while decreasing the magnitudes of the rest coefficients for $g_{4}({\bm{x}})$.

We should note that the energy-promoting behaviors of $g_{3}({\bm{x}})$ and $g_{4}({\bm{x}})$ are different from those of the proposed entropy functions. For the entropy functions, the thresholds $\nu$ in (32), (40) are adaptive in that they are updated every iteration along with the solution, as opposed to the rigid criterion induced by $g_{3}({\bm{x}})$ and $g_{4}({\bm{x}})$.

Appendix C Nonconvex Soft Shrinkage Thresholding

Take the optimization problem in (53) for example, let $\tau_{i}^{(t)}=\frac{\lambda}{\kappa}\nabla h_{p}(|x_{i}^{(t)}|)$, $s_{i}^{(t)}=x_{i}^{(t)}-\frac{1}{\kappa}\nabla f(x_{i}^{(t)})$, we have the following problem for the $i$-th entry $x_{i}$:

[TABLE]

• ❑ Positive threshold:

When $\tau_{i}^{(t)}\geq 0$, (89) is a convex problem. Its solution is given by applying the shrinkage operator given in (59) on $s_{i}^{(t)}$ with the threshold $\tau_{i}^{(t)}$.

• ❑ Negative threshold:

When $\tau_{i}^{(t)}<0$, (89) is a not necessarily a convex problem. Luckily this is a simple one dimensional problem, its global optimal solution can be still found as follows:

When $s_{i}^{(t)}<\tau_{i}^{(t)}$, as shown in Fig. A12:

• –

For $x_{i}\geq 0$, we have:

[TABLE]

Since $\tau_{i}^{(t)}-s_{i}^{(t)}>0$, the $x_{i}$ that minimizes (90) is [math].

• –

For $x_{i}<0$, we have:

[TABLE]

Since $-\tau_{i}^{(t)}-s_{i}^{(t)}>0$, the $x_{i}$ that minimizes (91) is $\tau_{i}^{(t)}+s_{i}^{(t)}$.

(89) is continuous at the point $x_{i}=0$. Hence the global minimum of (89) is obtained by $x_{i}=\tau_{i}^{(t)}+s_{i}^{(t)}$. 2. 2.

When $\tau_{i}^{(t)}\leq s_{i}^{(t)}<0$, as shown in Fig. A12:

• –

For $x_{i}\geq 0$, we have (90). Since $\tau_{i}^{(t)}-s_{i}^{(t)}<0$, the $x_{i}$ that minimizes (90) is $-\tau_{i}^{(t)}+s_{i}^{(t)}$.

• –

For $x_{i}<0$, we have (91). Since $-\tau_{i}^{(t)}-s_{i}^{(t)}>0$, the $x_{i}$ that minimizes (91) is $\tau_{i}^{(t)}+s_{i}^{(t)}$.

It’s easy to verify that the minimum of (91) is smaller than the minimum of (90). Hence the global minimum of (89) is obtained by $\tau_{i}^{(t)}+s_{i}^{(t)}$. 3. 3.

When $0\leq s_{i}^{(t)}<-\tau_{i}^{(t)}$, as shown in Fig. A12:

• –

For $x_{i}\geq 0$, we have (90). Since $\tau_{i}^{(t)}-s_{i}^{(t)}<0$, the $x_{i}$ that minimizes (90) is $-\tau_{i}^{(t)}+s_{i}^{(t)}$.

• –

For $x_{i}<0$, we have (91). Since $-\tau_{i}^{(t)}-s_{i}^{(t)}>0$, the $x_{i}$ that minimizes (91) is $\tau_{i}^{(t)}+s_{i}^{(t)}$.

It’s easy to verify that the minimum of (91) is larger than the minimum of (90). Hence the global minimum of (89) is obtained by $-\tau_{i}^{(t)}+s_{i}^{(t)}$. 4. 4.

When $s_{i}^{(t)}\geq-\tau_{i}^{(t)}$, as shown in Fig. A12:

• –

For $x_{i}\geq 0$, we have (90). Since $\tau_{i}^{(t)}-s_{i}^{(t)}<0$, the $x_{i}$ that minimizes (90) is $-\tau_{i}^{(t)}+s_{i}^{(t)}$.

• –

For $x_{i}<0$, we have (91). Since $-\tau_{i}^{(t)}-s_{i}^{(t)}<0$, the $x_{i}$ that minimizes (91) is [math].

(89) is continuous at the point $x_{i}=0$. Hence the global minimum of (89) is obtained by $x_{i}=-\tau_{i}^{(t)}+s_{i}^{(t)}$.

Combining the above $4$ different scenarios, we have the following results:

When $s_{i}^{(t)}\geq 0$, the solution to (89) is $s_{i}^{(t)}-\tau_{i}^{(t)}$. 2. 2.

When $s_{i}^{(t)}<0$, the solution to (89) is $s_{i}^{(t)}+\tau_{i}^{(t)}$.

This is exactly the soft shrinkage operator given in (59) on $s_{i}^{(t)}$ with the threshold $\tau_{i}^{(t)}$.

Appendix D First-Order Approximations of Other Regularization Functions

• •

$g_{2}({\bm{x}})=\sum_{i=1}^{N}|x_{i}|^{p}$, where ${\bm{x}}\in\mathbb{R}^{N}$. The function $g_{2}({\bm{x}})$ is continuous in its domain, hence it is also lower semi-continuous. The first-order approximation in (52) is updated accordingly:

[TABLE]

where $0<p<1$. The gradient with respect to $|x_{i}|$ is:

[TABLE]

• •

$g_{4}({\bm{x}})=\frac{\|{\bm{x}}\|_{1}}{N\cdot\max_{j}|x_{j}|}-1$, where ${\bm{x}}\in\mathbb{R}^{N}$ and ${\bm{x}}\neq\bm{0}$. The function $g_{4}({\bm{x}})$ is continuous in its domain. Hence it is also lower semi-continuous.

Since $g_{4}({\bm{x}})$ is not differentiable with respect to the modulus $|{\bm{x}}|$, we use its subderivative to get the first order approximation:

[TABLE]

The subderivative with respect to $|x_{i}|$ is:

[TABLE]

where $\delta_{\infty}(i)$ is given in (88).