Natural Thresholding Algorithms for Signal Recovery with Sparsity

Given $(y,A)$, the gradient of the error metric $\psi(x)=(1/2)\|y-Ax\|_{2}^{2}$ is given by $\nabla\psi(x)=-A^{T}(y-Ax).$ At the iterate $x^{(p)},$ let $u^{(p)}$ be the vector generated by the classic gradient descent method, i.e.,

where $\lambda_{p}$ is a steplength. Using the operator ${\cal H}_{k}$ to generate an iterate satisfying the constraint of the problem (1) leads to the well-known IHT algorithm [15, 25, 26] and its enhanced version HTP [27]. A good feature of IHT and HTP is that they are easier to implement than other algorithms such as convex optimization [38, 39, 40, 41], orthogonal matching pursuit (OMP) [42, 43, 44], compressive sampling matching pursuit (CoSaMP)[45], and subspace pursuit (SP) [46].

Clearly, the $k$ terms of $u^{(p)}$ selected by ${\cal H}_{k}$ are not necessarily the best $k$ terms in the sense that they might not be the $k$ terms that best fit the linear measurements. As pointed out in [34], performing the hard thresholding of $u^{(p)}$ is independent of the error metric $\psi.$ Thus it may happen that the value of the error metric at the point $x^{+}={\cal H}_{k}(u^{(p)})$ is larger than that of the current iterate $x^{(p)},$ i.e., $\|y-Ax^{+}\|_{2}>\|y-Ax^{(p)}\|_{2}.$ This observation motivates the recent development of optimal $k$-thresholding algorithms in [34, 35].

As defined in [34], the optimal $k$-thresholding of a vector $u$, denoted by $Z_{k}^{\#}(u),$ is given by $Z_{k}^{\#}(u)=u\otimes w^{*}$ which is the Hadamard product of $u$ and $w^{*},$ where $w^{*}$ is the solution to the problem

where e is the vector of ones, and the vector $w$ is $k$-sparse and binary with exactly $k$ entries equal to 1. Given such a binary vector $w,$ the Hadamard product of $u$ and $w$ retains $k$ entries of $u$ only and sets other entries to zeros. Based on the operator $Z_{k}^{\#}(\cdot),$ the conceptual optimal $k$-thresholding (OT) method and the optimal $k$-thresholding pursuit (OTP) were proposed in [34]. The two methods share the same first step except for the second step at which OTP performs an orthogonal projection.

We would rather treat OT and OTP as conceptual algorithms due to the fact that directly solving the binary problem (3) might be expensive when the dimension of the problem is high. A more practical way is to solve a convex relaxation of (3), which can be solved by existing interior-point methods. This consideration leads to the relaxed optimal $k$-thresholding (ROT) and relaxed optimal $k$-thresholding pursuit (ROTP) introduced in [34, 35].

While the optimization problem (4) can be solved efficiently by existing convex optimization solvers, the computational cost remains high for large-scaled signal recovery problems, compared with the algorithms which involve only orthogonal projection, hard thresholding or linear programming solving. This stimulates the natural thresholding (NT) method described in the next subsection.

The model (4) is a convex relaxation of (3). However, such a relaxation might be too loose, and thus its solution is usually not binary. To overcome such a drawback of (4), we may introduce a certain regularization term into the objective of (4) so that the regularized model can approximate the model (3) appropriately and the chance for the solution of regularized model being binary is enhanced. We now describe such a regularized model which is essential for the development of new algorithms in this paper. In fact, given vector $u^{(p)},$ consider the model of the form:

where $\alpha>0$ is a parameter and $\phi(w)$ is referred to as a binary regularization function which is formally defined as follows.

Given a vector $u^{(p)},$ we use $f$ and $g_{\alpha}$ to denote the functions

The function $g_{\alpha}(w)$ is the objective of the problem (II-C). The parameter $\alpha$ is referred to as the regularization parameter. The solution of (II-C) depends on the parameter $\alpha.$ Thus for any given $\alpha>0,$ we use $w(\alpha)$ to denote the solution of (II-C). $\phi(w)$ can be seen as the penalty for the variable $w$ violating the 0-1 restriction. From the classic optimization theory for penalty functions (see, e.g., [47]), it is not difficult to see that as $\alpha\to\infty,$ the value of $\phi(w(\alpha))$ will decrease to its minimum in the feasible set of (II-C). Such minimum attains only at 0-1 vectors according to the definition of binary regularization. This implies that $w(\alpha)$ tends to a 0-1 point in $D$ as $\alpha\to\infty,$ and any accumulation point of $w(\alpha)$ must be a solution to the problem (3). We summarize this fact in the following lemma whose proof is given in Appendix A for completeness.

The above lemma shows that the model (II-C) can approximate the model (3) to any expected level of accuracy provided that the parameter $\alpha$ is taken large enough. From this perspective, the model (II-C) as an approximation of (3) is better than the convex relaxation (4). By the definition, a binary regularization is nonconvex. As shown by (6), it is very easy to construct a binary regularization which is concave. As a result, the objective function $g_{\alpha}$ of (II-C) can be naturally made to be concave by properly choosing the parameter $\alpha.$ Thus for the convenience of discussion, we assume that $g_{\alpha}$ is concave throughout the remainder of this section and Section III. The concavity of the function $g_{\alpha}(w)$ can be guaranteed if $\phi$ and $\alpha$ are chosen properly. Note that if $\phi$ is twice continuously differentiable, then the Hessian of $g_{\alpha}$ is given as

from which it can be seen that when $\nabla^{2}\phi(w)$ is negative definite, then $\nabla^{2}g_{\alpha}(w)$ will be negative semi-definite (and thus $g_{\alpha}$ is concave) provided that $\alpha$ is large enough.

Take (6) as an example. It is easy to verify that the Hessian of (6) is $\nabla^{2}\phi(w)=-2I,$ where $I$ is the identity matrix. Thus when (6) is used, the Hessian of $g_{\alpha}(w)$ is given as

where $U^{(p)}=\textrm{diag}(u^{(p)}).$ Note that $\nabla^{2}g_{\alpha}(w)$ is negative semi-definite over the whole space provided that $\alpha I\succeq U^{(p)}A^{T}AU^{(p)}$ which can be guaranteed when $\alpha$ is not smaller than the largest eigenvalue of the matrix $U^{(p)}A^{T}AU^{(p)},$ i.e., $\alpha\geq\alpha^{*}:=\lambda_{\max}(U^{(p)}A^{T}AU^{(p)}).$ Therefore, the function $g_{\alpha}(w)$ is always concave over $\mathbb{R}^{n}$ provided that $\alpha\geq\alpha^{*}.$

Note that (3) is to find the optimal $k$-thresholding of $u^{(p)}.$ The above discussion indicates that performing the optimal $k$-thresholding of a vector can be approximately achieved by solving the problem (II-C) which from the above discussion can be made to be a concave minimization problem which remains NP-hard in general [49, 50]. The fundamental idea of the traditional MM method for solving a nonconvex optimization problem is to replace the ‘difficult’ objective function of the problem with a certain relatively easy surrogate convex function that dominates the original objective function (see, e.g., [51, 52]). For a concave objective, the simplest surrogate function is its first-order approximation at any given point. Thus the first-order approximation of the concave function can be used to develop an algorithm with significantly low computational cost.

Specifically, suppose that $w^{-}$ is the current iterate. At $w^{-},$ by the concavity of $g_{\alpha}(w)$ over an open neighborhood of $D,$ denoted by $O,$ one has

for any $w\in O$, where $\nabla g_{\alpha}(w^{-})$ denotes the gradient of $g_{\alpha}$ at $w^{-}.$ Specifically, if $\phi$ is given by (6), the gradient of $g_{\alpha}$ at $w^{-}$ is given as

The relation (8) shows that the first-order approximation of $g_{\alpha}(w)$ dominates $g_{\alpha}(w)$ on the whole open set $O.$ By the idea of MM methods, the problem (II-C), i.e., $\min\{g_{\alpha}(w):w\in P\},$ where

can be approximated by its surrogate counterpart which, in this case, is the linear programming (LP) problem

Removing constant terms in the objective above, the above LP problem is reduced to

Let $w^{+}$ be the solution to this problem. We may replace $w^{-}$ with $w^{+}$ and repeat the above process until a certain stopping criterion is met. This is the basic idea for the development of a new algorithm.

We now point out that there is no need to solve the LP problem (9) by an algorithm since an optimal solution of (9) can be obtained explicitly due to the nature of the polytope $P.$ In fact, according to the classic LP theory, the problem (9) must have an optimal solution, denoted by $w^{+},$ which is an extreme point of the polytope $P.$ As shown in [34], any extreme point of $P$ is $k$-sparse and binary. Hence the optimal value $\left[\nabla g_{\alpha}(w^{-})\right]^{T}w^{+}$ is actually the summation of $k$ entries of the vector $\nabla g_{\alpha}(w^{-}).$ Clearly, the smallest value of the function $\left[\nabla g_{\alpha}(w^{-})\right]^{T}w$ over the polytope $P$ is the sum of the smallest $k$ entries of $\nabla g_{\alpha}(w^{-}).$ Thus the optimal solution $w^{+}$ of the LP problem (9) can be explicitly given as follows. It is determined by the support of the $k$ smallest entries of $\nabla g_{\alpha}(w^{-})$, i.e., $w^{+}_{i}=1$ if $i\in{\cal S}_{k}(\nabla g_{\alpha}(w^{-}))$, otherwise, $w^{+}_{i}=0$ for $i=1,\dots,n,$ where ${\cal S}_{k}(\cdot)$ is the index set of the $k$ smallest entries of a vector. We now describe the algorithms which are referred to as the natural thresholding (NT) and the natural thresholding pursuit (NTP) in this paper.

The only difference between the two algorithms is that NT does not perform an orthogonal projection, but NTP does. At the beginning of every iteration of NT and NTP, the initial point $w^{-}$ can be any given $k$-sparse and binary vector. In particular, we may use the binary vector (10) as an initial point determined by $L_{k}(u^{(p)})$.

which together with the concavity of $g_{\alpha}$ implies that $f(w^{+})\leq f(w^{-}).$ (This fact will be used in the proof of Lemma III.2 later.) In fact, the concavity of $g_{\alpha}(w)$ and (12) implies that

i.e.,

Since $w^{+}$ and $w^{-}$ are 0-1 vectors at which $\phi$ achieves the minimum value $\phi(w^{+})=\phi(w^{-})$ over $D.$ Thus the above inequality reduces to $f(w^{+})\leq f(w^{-}).$ This means the iterate $w^{+}$ generated by NT or NTP is never worse than the one generated by IHT or HTP from the perspective of the error metric $f(w).$ If $w^{-}$ is not optimal to (9), then $\left[\nabla g_{\alpha}(w^{-})\right]^{T}w^{+}<\left[\nabla g_{\alpha}(w^{-})\right]^{T}w^{-},$ and hence $f(w^{+})<f(w^{-}).$ Thus in this case the algorithm finds the $k$ terms of $u^{(p)}$ which is better than the $k$ terms selected by the hard thresholding.

The NT and NTP algorithms perform linearization of $g_{\alpha}$ only once in every iteration. We may perform linearization more than once within their step S1, i.e., replacing the role of $w^{-}$ with $w^{+}$ and repeating the step S1 of NT more than once so that a point better than $w^{+}$ might be found. Such iterations within S1 are referred to as inner iterations. Repeating the inner iterations up to $q\geq 1$ times yields the algorithms NT_q and NTP_q.

Before we analyze the guaranteed performance of NT-type algorithms, it is necessary to construct more examples of concave binary regularization functions. We still construct such functions based on the kernel univariate one $\tau(t)=(t+1/2)(3/2-t)$ which is concave and positive over an open neighborhood of the interval $[0,1].$ Assume $g$ is another concave and strictly increasing univariate function over the interval $(0,\infty).$ Then the composite function $h(t)=g(\tau(t))$ must be a concave function over an open neighborhood of the interval $[0,1]$ due to the fact

where the last inequality follows from the fact $g^{\prime\prime}(\tau(t))\leq 0,\tau^{\prime\prime}(t)\leq 0$ and $g^{\prime}(\tau(t))\geq 0.$ We may use such a function to construct a binary regularization as follows:

The function constructed like this can be used in NT-type algorithms. We now give two such specific functions.

• •

  Let $g(t):=\log(1+t)$ which is concave and strictly increasing over the interval $(0,\infty).$ Then

  $$\phi(w)=\sum_{i=1}^{n}\log\left(1+(w_{i}+1/2)(3/2-w_{i})\right)$$

  (14)

  is a concave binary regularization over an open neighborhood of the region $D.$

• •

  Let $g(t):=\frac{t}{1+t}$ which is concave and strictly increasing over the interval $(0,\infty).$ Then

  $$\phi(w)=\sum_{i=1}^{n}\frac{(w_{i}+1/2)(3/2-w_{i})}{1+(w_{i}+1/2)(3/2-w_{i})}$$

  (15)

  is a concave binary regularization over an open neighborhood of the region $D.$

Similar to $f(w)$ which is related to $u^{(p)}$, the regularization $\phi$ can be also made related to $u^{(p)}.$ Such a regularization will be discussed in more detail in Section IV.