Sparse Recovery with Oversampling Ratio Greater than One Half

TL;DR

This paper demonstrates the theoretical feasibility of sparse recovery with an oversampling ratio greater than one half and introduces two modified algorithms to achieve such performance.

Contribution

It presents the first theoretical proof and practical algorithms for sparse recovery when the oversampling ratio exceeds one half.

Findings

Theoretical proof of recovery feasibility beyond half oversampling ratio.

Two modified HTP algorithms capable of high oversampling recovery.

Empirical validation of the proposed algorithms' effectiveness.

Abstract

Solving sparse recovery problem with high oversampling ratio is hard. We show that it is theoretically possible and we propose two modified HTP algorithms with such performances.

Equations

$$(P_x)\{\begin{array}{ll}{\min }_{z\in {\mathbb{R}}^N}& f(z)=\frac{1}{2}\Vert Az-b{\Vert }_2^2\\ \text{mbox}s.t.& \Vert z{\Vert }_0\le s,\end{array}$$

$$\Phi=\left\{x\in\mathbb{R}^{N}\big{|}\ \|x\|_{0}\leq s,\mbox{($P_{x}$) has multiple solutions}\right\}.$$

$$\Phi =\bigcup _{\begin{array}{l}{S}^{\prime },S^{\prime \prime }\subseteq \{1,\dots ,N\}\\ |S^{\prime }|,|S^{\prime \prime }|\le s\end{array}}\Phi (S^{\prime },S^{\prime \prime }),$$

$$\Phi (S^{\prime },S^{\prime \prime })=$$

$$\displaystyle\qquad\big{\{}x^{\prime}\in\mathbb{R}^{N}\big{|}\mathrm{supp}(x)=S^{\prime},\exists x^{\prime\prime}\neq x^{\prime},\mathrm{supp}(x^{\prime\prime})=S^{\prime\prime},A(x^{\prime}-x^{\prime\prime})=0\big{\}},$$

$$\dim\big{(}\mathrm{span}(\Phi(S^{\prime},S^{\prime\prime}))\big{)}\leq 2s-m<s.$$

$$\Phi =$$

$$=$$

$$A(x^{\prime }-x^{\prime \prime })\equiv Av=0.$$

$$\dim (\ker (A_T))=s^{\prime }+s^{\prime \prime }-l-m>0.$$

$$(v_T{)}_j=\{\begin{array}{ll}{x}_j^{\prime }& \text{mbox}ifj\in S^{\prime }\backslash S^{\prime \prime }\\ x_j^{\prime }-x_j^{\prime \prime }& \text{mbox}ifj\in S^{\prime }\cup S^{\prime \prime }\\ -x_j^{\prime \prime }& \text{mbox}ifj\in S^{\prime \prime }\backslash S^{\prime },\end{array}$$

$$\Psi (S^{\prime },S^{\prime \prime })=$$

$$\displaystyle\quad\quad\big{\{}z_{S^{\prime}\backslash S^{\prime\prime}}\in\mathbb{R}^{N}\big{|}\ z\in\mathbb{R}^{N},\mathrm{supp}(z)=T,z_{T}\in\ker(A_{T}),T=S^{\prime}\cup S^{\prime\prime}\big{\}}.$$

$$\dim (\Psi (S^{\prime },S^{\prime \prime }))\le \dim (\ker (A_T)).$$

$$x^{\prime}_{S^{\prime}\cap S^{\prime\prime}}\in\Omega(S^{\prime},S^{\prime\prime})=\big{\{}z_{S^{\prime}\cap S^{\prime\prime}}\big{|}\ z\in\mathbb{R}^{N}\big{\}},$$

$$x^{\prime }\in \Psi (S^{\prime },S^{\prime \prime })\oplus \Omega (S^{\prime },S^{\prime \prime }),$$

$$\dim (\Phi (S^{\prime },S^{\prime \prime }))$$

$$\le s^{\prime }+s^{\prime \prime }-l-m+l=s^{\prime }+s^{\prime \prime }-m\le 2s-m.$$

$$\mathrm{supp}(x^{\prime })=S^{\prime }=(S\backslash \{j\})\cup \{j^{\prime }\},$$

$$v\ne 0,A_{T}v=0.$$

$$x^{\prime }=x-v,$$

$$(P_{(A_2,x_2)}^{\prime })\{\begin{array}{ll}{\min }_{z\in {\mathbb{R}}^{N-s_1}}& \frac{1}{2}\Vert b-A_{2}z{\Vert }_2^2\\ \text{mbox}s.t.& \Vert z{\Vert }_0\le s-s_1,\end{array}$$

$$(P_{(A_2,x_2)}^{\prime \prime })\{\begin{array}{ll}{\min }_{z\in {\mathbb{R}}^{N-s_1}}& \frac{1}{2}\Vert b-A_1{\hat{x}}_1-A_{2}z{\Vert }_2^2\\ \text{mbox}s.t.& \Vert z{\Vert }_0\le s-s_1,\end{array}$$

$$P_1=I-A_1(A_1^T{A}_1{)}^{-1}{A}_1^T,$$

$$P_1{A}_1{x}_1+P_1{A}_2{x}_2=P_{1}b.$$

$$(P_{(A_2,x_2)})\{\begin{array}{ll}{\min }_{z\in {\mathbb{R}}^{N-s_1}}& \frac{1}{2}\Vert P_{1}b-P_1{A}_{2}z{\Vert }_2^2=\frac{1}{2}\Vert b-A_1{\hat{x}}_1-P_1{A}_{2}z{\Vert }_2^2\\ \text{mbox}s.t.& \Vert z{\Vert }_0\le s-s_1.\end{array}$$

$$(1-{\delta }_k^{(n)})x^{T}x\le x^T{A}^{(n)T}{A}^{(n)}x\le (1+{\delta }_k^{(n)})x^{T}x,\forall \Vert x{\Vert }_0\le k,$$

$$(1-{\delta }_{k-1}^{(n+1)}){\tilde{x}}^T\tilde{x}\le {\tilde{x}}^T{A}^{(n+1)T}{A}^{(n+1)}\tilde{x}\le (1+{\delta }_{k-1}^{(n+1)}){\tilde{x}}^T\tilde{x},\forall \Vert \tilde{x}{\Vert }_0\le k-1,$$

$${\delta }_{k-1}^{(n+1)}\le {\delta }_k^{(n)}.$$

$${\tilde{B}}_2=P^{(n)}{B}_2=(I-U_1{U}_1^T)B_2=U_2{U}_2^T{B}_2.$$

$$B=U{U}^{T}B=\left[\begin{array}{cc}{U}_1& U_2\end{array}\right]\left[\begin{array}{cc}{U}_1^T{B}_1& U_1^T{B}_2\\ 0& U_2^T{B}_2\end{array}\right],$$

$$\left[\begin{array}{cc}{U}_1^T{B}_1& U_1^T{B}_2\\ 0& U_2^T{B}_2\end{array}\right],$$

Full text

Sparse Recovery with Oversampling Ratio Greater than One Half

Wenbin Zhang

1 Introduction

Sparse recovery technique (or compressive sensing  [2, 6] technique) offers to recover signals under sparsity prior with fewer measurements than Nyquist compression would require. Mathematically, the standard sparse recovery problem can be formulated as the following:

[TABLE]

where sensing matrix $A\in\mathbb{R}^{m\times N}$, sparse signal $x\in\mathbb{R}^{N}(\|x\|_{0}\leq s)$, dummy variable $z\in\mathbb{R}^{N}$, and measurement vector $b=Ax\in\mathbb{R}^{m}$. Various algorithms have been developed to tackle $(P_{x})$ or varieties of $(P_{x})$ such as orthogonal matching pursuit [12], iterative hard thresholding [1], CoSaMP [11], hard thresholding pursuit [8], subspace pursuit [4], etc. Successful sparse recovery algorithms possess a list of properties [11] including minimal number of measurements required, adaptability to different sampling schemes, noise robusty, optimal error guarantees, and efficient computational resource usage.

In this paper, however, we are interested in designing algorithms that can successfully solve $(P_{x})$ with number $m$ of measurements as small as possible. This leads to the solvability of $\ell_{0}$ sparse recovery problem $(P_{x})$ with respect to oversampling ratio $\rho=s/m$ [5, 7]. Unfortunately, if we require $(P_{x})$ to be solvable for ALL $\|x\|_{0}\leq s$, then there is a theoretic upper-bound for oversampling ratio $\rho$. which is stated in the following theorem:

Theorem 1.1**.**

Given a sensing matrix $A\in\mathbb{R}^{m\times N}$ ($m\leq N$). The sparse recovery problem ($P_{x}$) has unique solution for all $\|x\|_{0}\leq s$ if and only if all the $m$-by-$2s$ sub-matrices of $A$ have full column ranks, implying that the oversampling ratio $s/m\leq 1/2$.

• Proof.

See [9], Theorem 2.13.∎

However, we observe that the phase transition behaviour of a certain modified Hard Thresholding Pursuit algorithm (named DP-HTP) seems to contradict theorem 1.1, see Fig. 1. Note that the large area in the phase transition map of DP-HTP with oversampling ratio $s/m\geq 1/2$ means DP-HTP algorithm successfully recovers the true sparse signals.

Here, we attempt to answer the following two questions:

1. 1

Is it possible for $\ell_{0}$ sparse recovery problem $(P_{x})$ to have unique solution, in some sense, with oversampling ratio $s/m$ larger than one half? 2. 2

How to explain the phase transition behaviour of DP-HTP algorithm?

The remainder of the paper is organized as follows. Section 2 gives a somehow positive answer to question 1 as well as a full theoretical result on uniqueness of solution to $(P_{x})$ with respect to oversampling ratio $\rho=s/m$. Section 3 presents the algorithm DP-HTP and justifies its phase-transition behaviour. Section 4 presents extended numerical experiments of HTP and DP-HTP with respect to different sparse signal $x$’s as well as different sensing matrix $A$s.

2 Uniqueness of solution when oversampling ratio greater than one half

The existence of solution to $(P_{x})$ is guaranteed since we generate this problem using a sparse vector $x$ and thus $x$ itself must be a solution. Our concern is the uniqueness.

It is trivial to observe that when oversampling ratio $\rho=s/m>1$, the solution to $(P_{x})$ is never unique since sub-matrix $A_{S}\in\mathbb{R}^{m\times s}$ has more columns than rows, where index set $S$ denotes the support of the true sparse signal $x$, and even if we know the true support, there are still infinitely many solutions to $(P_{x})$.

To prevent the above situation from happening, we introduce a basic assumption on the sensing matrix $A\in\mathbb{R}^{m\times N}$:

Assumption 2.1**.**

$A_{T}$* is non-singular for any index set $T\subset\{1,\ldots,N\}$ and $|T|\leq m$.*

Conditions on sub-matrix of the sensing matrix is quite common in compressive sensing community, such as the restricted isometry property(RIP) [3].Assumption 2.1 is just a weaker version of RIP. Also note that sensing matrix $A$ has overwhelming large probability to satisfy the basic Assumption 2.1 if $A$ is generated using usual random generating techniques such as random gaussian or random DCT  [9].

Since we already have negative answer for the uniqueness when $\rho=s/m>1$, the rest situations to be discussed are situation $1/2<\rho<1$ and $\rho=1$.

2.1 Situation $1/2<\rho<1$

We characterize the following set, which is the set of $s$-sparse $x$ such that ($P_{x}$) has multiple solutions:

[TABLE]

Theorem 2.2**.**

Given a fixed sensing matrix $A\in\mathbb{R}^{m\times N}$ and we suppose the basic Assumption 2.1 holds on $A$. If $1/2<\rho<1$, or $s<m<2s$, then the non-uniqueness set $\Phi$ has the following decomposition:

[TABLE]

where

[TABLE]

and

[TABLE]

• Proof.

We decompose $\Phi$ first to the support of $x^{\prime}$ and then further to the support of another solution:

[TABLE]

For any $x^{\prime}\in\Phi(S^{\prime},S^{\prime\prime})$, there exists an $x^{\prime\prime}\neq x^{\prime}$ with $\mathrm{supp}(x^{\prime\prime})=S^{\prime\prime}$ and

[TABLE]

We first let $s^{\prime}=|S^{\prime}|$, $s^{\prime\prime}=|S^{\prime\prime}|$, $T=S^{\prime}\cup S^{\prime\prime}$, $l=|S^{\prime}\cap S^{\prime\prime}|$ and then analyze the constraints imposed on $x^{\prime}$. We will only discuss indexes in support $T$ and this may help to eliminate the ambiguity of some symbols.

Of course, $|T|=s^{\prime}+s^{\prime\prime}-l$ and the existence of $x^{\prime\prime}$ implies

[TABLE]

The vector $v_{T}$ has a special structure in regard to $x^{\prime}$ and $x^{\prime\prime}$:

[TABLE]

and accordingly, we separate $x^{\prime}$ into two parts: $x^{\prime}_{S^{\prime}\backslash S^{\prime\prime}}$ and $x^{\prime}_{S^{\prime}\cup S^{\prime\prime}}$.

For index set $S^{\prime}\backslash S^{\prime\prime}$, we let

[TABLE]

Since $\Psi(S^{\prime},S^{\prime\prime})$ is some projection image of $\ker(A_{T})$, its dimension is no larger than that of $\ker(A_{T})$:

[TABLE]

Of course, $x^{\prime}_{S^{\prime}\backslash S^{\prime\prime}}\in\Psi(S^{\prime},S^{\prime\prime})$.

Now we consider index set $S^{\prime}\cap S^{\prime\prime}$. Although we have constraints on $x^{\prime}-x^{\prime\prime}$, $x^{\prime}$ itself on the joint support $S^{\prime}\cap S^{\prime\prime}$ is actually arbitrary. Thus we have

[TABLE]

of which the dimension is just $l$: $\dim(\Omega(S^{\prime},S^{\prime\prime}))=l$.

Since the full $x^{\prime}=x^{\prime}_{S^{\prime}\backslash S^{\prime\prime}}+x^{\prime}_{S^{\prime}\cup S^{\prime\prime}}$, we have

[TABLE]

and

[TABLE]

As long as $s<m$, the dimension of $\Phi(S^{\prime},S^{\prime\prime})$ is less than $s$, as is stated in the theorem. ∎

• Remark

Those “bad” $x$’s are contained in a finite union of subspaces of dimension smaller than $s$. If we use random generators such gaussian or uniform, for true sparse signal $x$, the probability of encountering a “bad” $x$ is actually zero.

2.2 Situation $\rho=1$

This situation is a little bit tricky: we can not simply assert that every $(P_{x})$ has multiple solutions since $A_{T}$ is square and nonsingular according to the basic Assumption 2.1; the technique in situation $1/2<\rho<1$ does not apply here neither since it is only an inequality. But we can still draw a definitive conclusion with a more careful investigation, shown in below.

Theorem 2.3**.**

Given a fixed sensing matrix $A\in\mathbb{R}^{m\times N}$ and we suppose the basic Assumption 2.1 holds on $A$. If $\rho=1$, or $s=m$, then every sparse recovery problem $(P_{x})$ with $\|x\|_{0}\leq s$ has at lease two different solutions.

• Proof.

We first consider $\|x\|_{0}=s$. For $x$ with support $S$, we construct another solution $x^{\prime}$ to $(P_{x})$ of the following special form:

[TABLE]

with some $j\in S$ and $j^{\prime}\in\{1,\ldots,N\}\backslash S$.

As usual, we let $T=S\cup S^{\prime}=S\sqcup\{j^{\prime}\}$ and $|T|=s+1$. We fix $j^{\prime}\in\{1,\ldots,N\}\backslash S$ and wish to find a $j\in S$ that satisfies our demand. Note that $T$ will not change any more once $j^{\prime}$ is fixed.

$A_{T}\in\mathbb{R}^{m\times(m+1)}$ has a one-dimensional null space and we denote $v$ as a representative:

[TABLE]

We now face two situations: there is at least one nonzero entry in $v$ with index lying in $S$ or the only nonzero entry is $v_{j^{\prime}}$. We claim that the latter situation never happens. If that happens, we have $v_{i}=0$ for all $i\in S$ and $v_{j}\neq 0$. This can only lead to $j^{\prime}$-th column of $A$ equalling to zero, which contradicts with the basic assumption on $A$.

Now we have $j\in S$ and $v_{j}\neq 0$. Without loss of generality, we assume $v_{j}=x_{j}\neq 0$, otherwise we can just scale $v$ to be so. Let

[TABLE]

then $x^{\prime}_{j}=0$, which means $\|x^{\prime}\|_{0}\leq s$ and $x^{\prime}\neq x$.

If $\|x\|_{0}<s$, we can artificially complement $\mathrm{supp}(x)$ to have length $s$, and the rest of the proof is exactly the same as above. ∎

2.3 Summary on uniqueness with respect to oversampling ratio

If we fix number $s$ of non-zero entries in true sparse signal $x$, the situation gets worse as the number $m$ of measurements decreases. First, the strongest condition $m\geq 2s$ yields the “strongest” uniqueness: the solution to $(P_{x})$ is unique for ALL $\|x\|_{0}\leq s$. Next, as the condition gets weaker ($1/2<s/m<1$), the set of “bad” $x$’s becomes a union of some low-dimensional subspaces and as number of measurements $m$ decreases, their dimensions get larger (Theorem 2.2). Then $m$ decreases to be the same as $s$ ($s/m=1$) and every $(P_{x})$ has multiple solutions (Theorem 2.3), though we still need to construct a different solution using a different support. Finally, when $s/m>1$, even if we know the support $S$ of true sparse signal $x$, there are still infinitely many solutions to ($P_{x}$) since sub-matrix $A_{S}$ itself is a “fat” matrix and thus column rank deficient.

3 Deflation-and-projection-HTP (DP-HTP) algorithm

In Figure 1, we observe that the success region of DP-HTP algorithm is significantly larger than that of pure HTP algorithm. In the last section, theoretical analysis shows that it is possible for an algorithm to success even when $1/2<\rho<1$ (the sparse recovery problem ($P_{x}$) has unique solution for almost all $x$).

In this section, we state deflation-and-projection-HTP algorithm and the rationale behind. Then we make attempts to justify the transition map behaviour of DP-HTP algorithm via “fast descendent” condition on true sparse signal.

The HTP algorithm [8] serves as a critical ingredient in modified algorithm DP-HTP and thus we state it here in Algorithm 1 for later use.

3.1 Statement of DP-HTP algorithm

If we know part of the solution, we may be able to modify the original problem into an easier one. Deflation is a term for this philosophy. In order to apply this idea, we need to investigate the sparse recovery problem from a different angle: instead of solving $x$ in problem ($P_{x}$), we only need to find its support $S=\mathrm{supp}(x)$. Now if we know in advance some indices in support $S$, how do we find the rest ones?

Formally, let $A=[A_{1},A_{2}]$, $x=[x_{1}^{T},x_{2}^{T}]^{T}$, $Ax=A_{1}x_{1}+A_{2}x_{2}=b$, $x_{1}\in\mathbb{R}^{s_{1}},s_{1}<s$, and suppose $A_{1}$ is known. (Note that we do not suppose $x_{1}$ is known.) Also we assume that all the elements in $x_{1}$ are non-zero (thus only a few entries in $x_{2}$ are non-zero since $x$ itself is sparse).

We can not just solve

[TABLE]

since $y\neq A_{2}x_{2}$. Nor can we solve

[TABLE]

where $\hat{x}_{1}=A_{1}^{\dagger}b$, since $y\neq A_{1}\hat{x}_{1}+A_{2}x_{2}$.

But we can do the following: let $P_{1}\in\mathbb{R}^{m\times m}$ be the orthogonal projection onto the orthogonal complement of $\mbox{span}(A_{1})$, i.e.,

[TABLE]

and then multiply $P_{1}$ on both sides of $A_{1}x_{1}+A_{2}x_{2}=b$ to get

[TABLE]

Since $P_{1}A_{1}\equiv 0$, we obtain the following size-reduced sparse recovery problem:

[TABLE]

Note that $P_{1}b$ is just the residual after least squares solving on $A_{1}$.

The difference between ($P^{\prime\prime}_{(A_{2},x_{2})}$) and ($P_{(A_{2},x_{2})}$) is that the sensing matrix in the latter problem is $P_{1}A_{2}$, a projected one, rather than $A_{2}$.

Motivated by the above analysis, we design a three-step iteration below:

• 1

Try HTP on the sparse recovery problem;

• 2

Select an index based on the result given by HTP;

• 3

Use projection to shrink the problem into a smaller one.

Formally, we write the above iteration steps as the DP-HTP algorithm:

3.2 Justification of DP-HTP algorithm

For terminology simplification, we denote from now on true sparse signal by $x^{\ast},\|x^{\ast}\|\leq s$, its support $S^{\ast}=\mathrm{supp}(x^{\ast})$, and measurement vector $b=Ax^{\ast}$. The DP-HTP algorithm calls HTP algorithm in each iteration, picks up the largest entry in the result and then projects the problem into an “one-order-smaller”(using smaller is better?) one.

There are two conditions for DP-HTP algorithm to work: that the HTP step always returns a result (though it may not be the correct solution) and that the index-selection step always picks up a correct index. Although our numerical experiments show that the iterate in HTP algorithm never become periodic in practice, there seem to be no theoretical characterisation for the situation here. Thus we leave the former condition as an assumption. For the latter one, out attempt is to impose assumptions on both the sensing matrix and the distribution of the non-zero entries in the true sparse signal.

For the sensing matrix part, we have the following theoretical result:

Theorem 3.1**.**

The $(k-1)$-th order restricted isometry constant of $A^{(n+1)}$ can be controlled by $k$-th order RIC of $A^{(n)}$, i.e., if there exists $\delta_{k}^{(n)}<1$ s.t.,

[TABLE]

then there exists $\delta_{k-1}^{(n+1)}<1$ such that

[TABLE]

and

[TABLE]

• Proof.

Verifying restricted isometry property of $A^{(n+1)}$ involves estimating the singular values of sub-matrices with $k-1$ columns.

Choose any $k-1$ columns of $A^{(n+1)}$, denoted by $\tilde{B}_{2}\in\mathbb{R}^{m\times(k-1)}$. Find the corresponding columns in $A^{(n)}$, and name them as $B_{2}\in\mathbb{R}^{m\times(k-1)}$. Through the definition of $A^{(n+1)}$, we have $\tilde{B}_{2}=P^{(n)}B_{2}$. Let $B_{1}=A_{i^{(n)}}^{(n)}$ and $B=[B_{1},B_{2}]$. We are to bound the singular values of $\tilde{B_{2}}$ by those of $B$.

Let $[U_{1},U_{2}]$ be a set of orthogonal basis of $\mbox{span}(B)$ with $U_{1}=B_{1}/\|B_{1}\|_{2}$. Since $P^{(n)}$ is the projection onto the orthogonal complement of $\mbox{span}(B_{1})$,

[TABLE]

Write

[TABLE]

and we have

• –

The $k$ singular values of $B$ are exactly the same as those of

[TABLE]

• –

The $(k-1)$ singular values of $\tilde{B}_{2}=U_{2}U_{2}^{T}B_{2}$ are exactly the same as those of $U_{2}^{T}B_{2}$.

By Cauchy interlacing theorem [10], the ordered singular values of $U_{2}^{T}B_{2}$ can be inserted into the order singular values of $B$, which implies

[TABLE]

and thus implies $\delta_{k-1}^{(n+1)}\leq\delta_{k-1}^{(n+1)}$. ∎

Theorem 3.2**.**

The $(k-n)$-th order RIC $\delta_{k-n}^{(n)}$ of $A^{(n)}$ can be controlled by the $k$-order RIC $\delta_{k}$ of $A$, i.e., $\delta_{k-n}^{(n)}\leq\delta_{k}$.

Theorem 3.2 asserts that the smaller sensing matrix is never worse than the previous one in the sense of restricted isometry property. Our next mission is to derive conditions that ensure the index of the largest entry in the result given by HTP algorithm corresponds to some index lying in the support of true solution $\mathrm{supp}(x^{\ast})$. Note that it is impossible to develop any algorithm that can work for ALL $x^{\ast},\|x^{\ast}\|\leq s$ when $m<2s$, since the sparse recovery problem itself may have multiple solutions then (see Theorem 2.2). Thus it is natural to assume constraints on the true sparse signal $x^{\ast}$. Our approach is to analysis one(the first) iteration in DP-HTP algorithm and see whether it is compatible with mathematical induction.

Let $x=\mbox{HTP}(A,b,s)$ be the result given by HTP algorithm (the first iteration in DP-HTP algorithm), and $S=\mathrm{supp}(x)$ be its support. Furthermore, we give three conditions below:

Assumption 3.3**.**

$A$* has $(s+2)$-th order restricted isometry constant $\delta=\delta_{s+2}<1$. Its detailed bound will be derived later.*

Assumption 3.4**.**

HTP algorithm on $(A,b,s)$ converges (may not converge to true sparse signal $x^{\ast}$).

Assumption 3.5**.**

The largest entry in $x^{\ast}$ is dominant by the following definition:

[TABLE]

Our aim is to derive some condition involving restricted isometry constant $\delta$ and dominance factor $\gamma$ that can ensure $l\in S^{\ast}$.

Let $u=x+A^{T}(b-Ax)$. According to Assumption 3.4, we only need to ensure

[TABLE]

By the definition of HTP algorithm, the nonzero entries $x_{S}$ are solved via least squares:

[TABLE]

Also we denote the orthogonal project onto the orthogonal complement of $\mathrm{span}(A_{S})$ by

[TABLE]

for latter use.

Proposition 3.6**.**

Let $B=[B_{1},B_{2}]$, $B_{1}\in\mathbb{R}^{m\times s}$, $B_{2}\in\mathbb{R}^{m\times 1}$ ($m\geq s+1$), and orthogonal projection $P=I-B_{1}(B_{1}^{T}B_{1})^{-1}B_{1}^{T}$. Suppose $\|B^{T}B-I\|_{2}\leq\delta$, or equivalently, singular values of $B$ are bounded by $[\sqrt{1-\delta},\sqrt{1+\delta}]$. Then $1-\delta\leq\langle B_{2},PB_{2}\rangle\leq 1+\delta$.

• Proof.

Using Gram-Schmidt orthogonalization, we can assume that the columns of $U\in\mathbb{R}^{m\times(s+1)}$ are the orthogonal basis of $\mathrm{span}(B)$, $U=[U_{1},U_{2}]$, and columns of $U_{1}\in\mathbb{R}^{m\times s}$ are the orthogonal basis of $B_{1}$, which yields $U_{2}^{T}B_{1}=0$. Write

[TABLE]

$U^{T}_{2}B_{2}$ (a scalar) is the $2$-by-$2$ block, and by the Cauchy interlacing theorem, it can be bounded by the singular values of $B$: $\sqrt{1-\delta}\leq U^{T}_{2}B_{2}\leq\sqrt{1+\delta}$. Thus $\langle B_{2},PB_{2}\rangle=B_{2}^{T}U_{2}U_{2}^{T}B_{2}\in[1-\delta,1+\delta]$. ∎

Proposition 3.7**.**

Let $B=[B_{1},B_{2},B_{3}]$, $B_{1}\in\mathbb{R}^{m\times s}$, $B_{2},B_{3}\in\mathbb{R}^{m\times 1}$ ($m\geq s+2$), and orthogonal projection $P=I-B_{1}(B_{1}^{T}B_{1})^{-1}B_{1}^{T}$. Suppose $\|B^{T}B-I\|_{2}\leq\delta$, or equivalently, singular values of $B$ are bounded by $[\sqrt{1-\delta},\sqrt{1+\delta}]$. Then $|\langle B_{3},PB_{2}\rangle|\leq\delta$.

• Proof.

Similar to the proof in proposition 3.6, we assume there exists orthonormal $U=[U_{1},U_{23}]\in\mathbb{R}^{m\times(s+1)}$, $U_{1}\in\mathbb{R}^{m\times s}$, $\mathrm{span}(U)=\mathrm{span}(B)$, and $\mathrm{span}(U_{1})=\mathrm{span}(B_{1})$. Write identity

[TABLE]

$C=U_{23}^{T}B_{23}\in\mathbb{R}^{2\times 2}$ is the 2-by-2 block and thus has bound $\|C^{T}C-I\|_{2}\leq\delta$. Again we separate blocks of $C^{T}C-I$ by

[TABLE]

then

[TABLE]

since the 2-norm of a sub-matrix is bounded by that of the full matrix.∎

Suppose $l\in S^{\ast}\backslash S$. We are to derive a constraint involving $\delta$ and $\gamma$ that will yield $|u|_{l}>|u|_{j},\ \forall j\in\overline{S\cup S^{\ast}}$, which will in turn lead to contradiction. To separate the term containing $x^{\ast}_{l}$, we write

[TABLE]

By Proposition 3.6, term $\langle A_{l},P_{S}A_{l}\rangle$ in (36) can bounded from below by

[TABLE]

then $|u_{l}|$ has a lower bound

[TABLE]

The last inequality is derived via Assumption 3.5.

Using similar techniques and proposition 3.7, we derive the upper bound of $|u_{j}|$, $j\in\overline{S\cup S^{\ast}}$:

[TABLE]

If the dominance factor $\gamma$ satisfies

[TABLE]

then $|u_{l}|>|u_{j}|$, $\forall j\in\overline{S\cup S^{\ast}}$, implying contradiction on the supposition $l\in S^{\ast}\backslash S$.

Now we turn to the case when $l\in S\cap S^{\ast}$, and we are to find conditions to assure $|u_{l}|>|u_{i}|$, $\forall i\in S\backslash S^{\ast}$. Similarly, we derive lower bound for $u_{l}$ and upper bound for $u_{i}$. On one hand,

[TABLE]

Considering the fact that the singular values of $(A_{S}^{T}A_{S})^{-1}A_{S}^{T}$ are the reciprocal of those of $A_{S}$, we have

[TABLE]

On the other hand, $\forall i\in S\backslash S^{\ast}$, we have similarly,

[TABLE]

Again, as long as the dominance factor $\gamma$ satisfies

[TABLE]

then $|u_{l}|>|u_{i}|$, $\forall i\in S\backslash S^{\ast}$. Combining inequalities (43) and (50) together, we conclude that the index of the largest entry in $x$, though may not be $l$, is assured to lie in the support $S^{\ast}$ of the true sparse signal, if the following inequality holds:

[TABLE]

• Remark

If we simply relax the factor $\|x^{\ast}_{S^{\ast}\backslash(S\cup\{l\})}\|_{2}/\|x^{\ast}_{S^{\ast}\backslash\{l\}}\|_{2}$ to be one, the dominance factor $\gamma$ will need to be larger than $2(1+\delta)/(1-2\delta)\geq 2$ to guarantee validity, according to (51). This makes the assumption 3.5 too strong to satisfy. However, if the support $S$ given by HTP algorithm already contains a large part of the true support $S^{\ast}$, the factor $\|x^{\ast}_{S^{\ast}\backslash(S\cup\{l\})}\|_{2}/\|x^{\ast}_{S^{\ast}\backslash\{l\}}\|_{2}$ can be expected to be far smaller than one, making it possible for a random sparse signal with Gaussian distribution to satisfy Assumption 3.5. The heuristic observations on factor $\|x^{\ast}_{S^{\ast}\backslash(S\cup\{l\})}\|_{2}/\|x^{\ast}_{S^{\ast}\backslash\{l\}}\|_{2}$ are discussed in numerical examples in the next section. (need completion)

• Remark

Comparison with OMP-type strategy. The DP-HTP algorithm utilizes the result given by HTP algorithm. One may wonder whether HTP is necessary: after all, each HTP requires 10-20 iterations and each iteration involves a least squires problem solving. Here, we discuss a simplified version: using an Orthogonal Matching Pursuit (OMP)-type iteration instead of HTP iteration in each DP-HTP step (see step 2 in algorithm 2). Since OMP uses $A^{T}Ab$ as the support indicator, we estimate

[TABLE]

and $\forall j\neq l$,

[TABLE]

Thus the condition for ensuring $|u_{l}|>|u_{j}|$ should be

[TABLE]

The fundamental difference between (51) and (58) is the factor

[TABLE]

which is essential for Assumption 3.5 to be satisfied on randomly generated true signal $x^{\ast}$, with, for example, Gaussian distribution.

• Remark

The effectiveness of the DP-HTP algorithm heavily relies on the result given by the HTP algorithm in every step, as can be seen from the definition of DP-HTP and the analysis above.

4 Numerical Experiments

4.1 Gaussian Setting

We set the experiment as follows:

• •

Signal length $N=800$. Number of measurements $m$ and number of nonzeros $s$ are set such that both oversampling ratio $\rho=s/m$ and undersampling ratio $\delta=m/N$ vary from $0.2$ to $0.8$.

• •

Both the entries of sensing matrix $A$ and the nonzero entries of true solution $x$ are $i.i.d$ Gaussians with unit standard variation.

• •

The experient is repeated 1000 times independently.

• •

The algorithms in comparison are OMP, Subspace Pursuit and HTP, as well as modified HTP algorithms proposed in this paper.

The overall phase transition map for the five algorithms in comparison can be found in Fig.2 and Fig.3 shows a detailed comparison for 95% success rate performance. The two modified HTP algorithms TR-HTP and DP-HTP perform far better on phase transition behavior than the rest three greedy algorithms do. Also when undersampling ratio $\delta$ is large (close to $0.8$), TR-HTP and DP-HTP manage to solve $P_{x}$ with oversampling ratio $\rho$ much larger than $0.5$.