Iterative hard thresholding for compressed data separation

Abstract

We study the problem of reconstructing signals’ distinct subcomponents, which are approximately sparse in morphologically different dictionaries, from a small number of linear measurements. We propose an iterative hard thresholding algorithm adapted to dictionaries. We show that under the usual assumptions that the measurement system satisfies a restricted isometry property (adapted to a composed dictionary) condition and the dictionaries satisfy a mutual coherence condition, the algorithm can approximately reconstruct the distinct subcomponents after a fixed number of iterations.

Introduction

The basic perceptiveness of compressed sensing is that sparse signals can be reconstructed via efficient algorithms from a small number of linear measurements.

In standard compressed sensing, one observes $(\mathbf{A},\mathbf{y})$ with

$$\mathbf{y}=\mathbf{A}\mathbf{f}+\mathbf{e},$$

where $\mathbf{y}\in {\mathbb{R}}^m$, $\mathbf{A}\in {\mathbb{R}}^{m\times n}$ with $m\ll n$, $\mathbf{f}\in {\mathbb{R}}^n$ is a (approximately) $s$-sparse signal of interest, and $\mathbf{e}\in {\mathbb{R}}^m$ is a vector of measurement errors. The goal is to reconstruct the sparse signal $\mathbf{f}$ based on the measurement matrix $\mathbf{A}$ and the measurement vector $\mathbf{y}$ via an efficient algorithm. The constrained $l_1$-minimization has been shown to be very effective for these problems. See, e.g., Candès and Tao [9], [11], Donoho [15] and Donoho, Elad, and Temlyakov [16].

Though the constrained $l_1$-minimization can be solved efficiently using iterative algorithms from convex optimization, it may suffer substantial computation expense in large-scale applications [18], [33]. Thus it is necessary to use alternative iterative methods that are not based on optimization, such as orthogonal matching pursuit (OMP) [14], [30], stagewise OMP [19], regularized OMP [29], compressive sampling matching pursuit (CoSaMP) [28], iterative hard thresholding (IHT) [4], subspace pursuit [13] and many other variants [5], [21]. Recovery results based on a restricted isometry property (RIP) [9] condition ${\delta }_{cs}\le \theta$ have been well developed for these algorithms, see [4], [5], [7], [13], [21], [28], [34], [35]. The RIP condition ${\delta }_{cs}\le \theta$ is satisfied with high probability for several random matrix ensembles, such as the subgaussian ensembles and random partial Fourier transforms [3], [10], [31], provided one chooses $m\ge O(s{log}^{\beta }n∕{\theta }^2)$. Consequently, with high probability, the aforementioned algorithms can approximately recover every $s$-sparse vector with small or zero errors from $O\left(s{log}^{\beta }\left(n\right)\right)$ random measurements.

For signals which are sparse in some orthonormal basis the above techniques hold. However, in practical examples, there are numerous signals of interest which are not sparse in an orthonormal basis. Often, sparsity is expressed not in terms of an orthogonal basis but in terms of an overcomplete (and possibly coherent) dictionary, which means that our signal $\mathbf{f}$ is now expressed as $\mathbf{f}=\mathbf{D}\mathbf{x}$ where $\mathbf{D}\in {\mathbb{R}}^{n\times d}(d\ge n)$ is a redundant dictionary and $\mathbf{x}$ is (approximately) sparse, see e.g. [8], [12] and references therein.

The $l_1$-analysis approach is one of the effective approach for solving these problems. Analogous recovery results based on $\mathbf{D}$-RIP (RIP adapted to the dictionary $\mathbf{D}$) as those in the classical setting have been developed for these approaches [1], [8], [23], [27], [32] . Specially, when $\mathbf{D}$ is a tight frame, Candès et al. [8] showed that if the measurement matrix satisfies a $\mathbf{D}$-RIP condition ${\delta }_{2s}\le 0.08$, the solution $\hat{\mathbf{f}}$ of the analysis Basis Pursuit (ABP)

$$\underset{\tilde{\mathbf{f}}\in {\mathbb{R}}^d}{arg\; min}{\Vert {\mathbf{D}}^{\ast }\tilde{\mathbf{f}}\Vert }_1\text{s.t.}{\Vert \mathbf{A}\tilde{\mathbf{f}}-\mathbf{y}\Vert }_2\le ϵ$$

has an error bound

$${\Vert \hat{\mathbf{f}}-\mathbf{f}\Vert }_2\le c_0\frac{{\sigma }_s\left({\mathbf{D}}^{\ast }\mathbf{f}\right)}{\sqrt{s}}+c_1ϵ,$$

provided that ${\Vert \mathbf{e}\Vert }_2\le ϵ$. Here, ${\sigma }_s\left(\mathbf{x}\right)$ denotes the best $s$-term approximation error in $l_1$-norm:

$${\sigma }_s\left(\mathbf{x}\right)=\underset{{\Vert \mathbf{u}\Vert }_0\le s}{min}{\Vert \mathbf{u}-\mathbf{x}\Vert }_1.$$

Under the assumption that the measurement matrix $\mathbf{A}$ satisfies the $\mathbf{D}$-RIP condition ${\delta }_{3s}\le 1∕2$, Lin and Li [23] showed that the solution $\hat{\mathbf{f}}$ of the analysis Dantzig selector (ADS)

$$\underset{\tilde{\mathbf{f}}\in {\mathbb{R}}^d}{arg\; min}{\Vert {\mathbf{D}}^{\ast }\tilde{\mathbf{f}}\Vert }_1\text{s.t.}{\Vert {\mathbf{D}}^{\ast }{\mathbf{A}}^{\ast }(\mathbf{A}\tilde{\mathbf{f}}-\mathbf{y})\Vert }_{\infty }\le \lambda$$

can recover the signal with an error bound

$${\Vert \hat{\mathbf{f}}-\mathbf{f}\Vert }_2\le c_0\frac{{\sigma }_s\left({\mathbf{D}}^{\ast }\mathbf{f}\right)}{\sqrt{s}}+c_1\sqrt{s}\lambda ,$$

provided that ${\Vert {\mathbf{D}}^{\ast }{\mathbf{A}}^{\ast }\mathbf{e}\Vert }_{\infty }\le \lambda$. Recall that $\mathbf{D}$-RIP [8] is defined as follows.

Definition 1

The measurement matrix $\mathbf{A}$ satisfies restricted isometry property adapted to $\mathbf{D}\in {\mathbb{R}}^{n\times d}$ (abbreviated as $\mathbf{D}$-RIP) of order $s$ with constant $\delta \in (0,1)$, if for all $s$-sparse vectors $\mathbf{z}\in {\mathbb{R}}^d$,

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

The restricted isometry constant adapted to $\mathbf{D}$ (abbreviated as $\mathbf{D}$-RIC) of order $s$ is the smallest number such that (1.1) holds for all $s$-sparse vectors $\mathbf{z}\in {\mathbb{R}}^d$, and it is denoted as ${\delta }_s$.

As noted in [8], the $\mathbf{D}$-RIP condition ${\delta }_{cs}\le \theta$ is satisfied with high probability by $m\times n$ matrices populated by i.i.d subgaussian entries with variance $m^{-1}$ provided $m\ge c\left(\delta \right)sln(ed∕s)$. It is also fulfilled with high probability for random partial Fourier matrices after sign randomization of their columns [22].

In signal processing and statistics, it is common to assume that the noise vector is a Gaussian noise, $\mathbf{e}\sim \mathcal{N}(0,{\sigma }^2{\mathbf{I}}_m)$. The Gaussian noise is essentially bounded (e.g. [11], [23]). Thus, all derived results mentioned above for the ABP and the ADS can be applied directly to the Gaussian noise. In this case, the ABP and the ADS provide very similar guarantees, but there are certain circumstances that the ADS is preferable since the ADS yields a bound that is adaptive to the unknown level of sparsity of the object signal and thus providing a stronger guarantee when $s$ is small [23] (this is also noted by Candès and Tao [11] for the classical compressed sensing).

Most recently, Foucart [20] studied IHT adapted to a dictionary and showed that under a $\mathbf{D}$-RIP condition, IHT provides the same theoretical guarantees as ABP considering an $l_2$-bounded noise. In this paper, as a byproduct of our analysis, we will show that IHT has the same theoretical guarantees as ADS, yielding preferable error bounds (that are adaptive to the unknown level of sparsity of the object signal) than those from [20] for the Gaussian noise.

The main goal of this paper is to study the problem of compressed data separation, i.e., reconstruction of signal’s different sparse components from compressed measurements. In the latter case, the signal is composed of two different components, i.e. $\mathbf{f}={\mathbf{f}}_1+{\mathbf{f}}_2$. More specifically, we observe the data $\mathbf{y}\in {\mathbb{R}}^m$ following the linear measurement model

$$\mathbf{y}=\mathbf{A}({\mathbf{f}}_1+{\mathbf{f}}_2)+\mathbf{e}.$$

The goal is to reconstruct the unknown constituents ${\mathbf{f}}_1$ and ${\mathbf{f}}_2$ based on the measurement vector $\mathbf{y}$ and the measurement matrix $\mathbf{A}$. Refer to [25] and references therein for further details on compressed data separation. As in [17], [25], we assume that the different components ${\mathbf{f}}_1$ and ${\mathbf{f}}_2$ are sparse or approximately sparse in terms of two different tight frames ${\mathbf{D}}_1\in {\mathbb{R}}^{n\times d_1}$ and ${\mathbf{D}}_2\in {\mathbb{R}}^{n\times d_2}$, respectively, i.e., ${\mathbf{f}}_1={\mathbf{D}}_1{\mathbf{x}}_1$ and ${\mathbf{f}}_2={\mathbf{D}}_2{\mathbf{x}}_2$ for some (approximately) sparse vectors ${\mathbf{x}}_1\in {\mathbb{R}}^{d_1}$ and ${\mathbf{x}}_2\in {\mathbb{R}}^{d_2}$.¹ The morphological difference between components ${\mathbf{f}}_1$ and ${\mathbf{f}}_2$ is measured in terms of the mutual coherence, defined as follows.

Definition 2

The mutual coherence between two dictionaries ${\mathbf{D}}_1$ and ${\mathbf{D}}_2$ is defined as

$$\mu =\underset{i,j}{max}\left|〈{\left[{\mathbf{D}}_1\right]}_i,{\left[{\mathbf{D}}_2\right]}_j〉\right|,$$

where ${\left[{\mathbf{D}}_1\right]}_i$ and ${\left[{\mathbf{D}}_2\right]}_j$ are the $i$th column of ${\mathbf{D}}_1$ and $j$th column of ${\mathbf{D}}_2$.

The $l_1$-split analysis approach was proposed in [2], [8], which finds ${\mathbf{f}}_1$ and ${\mathbf{f}}_2$ via the $l_1$-minimization

$$\underset{{\tilde{\mathbf{f}}}_1,{\tilde{\mathbf{f}}}_2}{arg\; min}{\Vert {\mathbf{D}}_1^{\ast }{\tilde{\mathbf{f}}}_1\Vert }_1+{\Vert {\mathbf{D}}_2^{\ast }{\tilde{\mathbf{f}}}_2\Vert }_1\text{s.t.}{\Vert \mathbf{A}({\tilde{\mathbf{f}}}_1+{\tilde{\mathbf{f}}}_2)-\mathbf{y}\Vert }_2\le ϵ.$$

The $l_1$-split analysis [24], [25] can approximately reconstruct the different two components ${\mathbf{f}}_1$ and ${\mathbf{f}}_2$ with an error bound

$${\Vert {\hat{\mathbf{f}}}_1-{\mathbf{f}}_1\Vert }_2+{\Vert {\hat{\mathbf{f}}}_2-{\mathbf{f}}_2\Vert }_2\le C_0\left(\frac{{\sigma }_{s_1}\left({\mathbf{D}}_1^{\ast }{\mathbf{f}}_1\right)+{\sigma }_{s_2}\left({\mathbf{D}}_2^{\ast }{\mathbf{f}}_2\right)}{\sqrt{s_1+s_2}}\right)+C_1ϵ,$$

provided that the measurement matrix satisfies a $\mathit{\Phi }$-RIP condition ${\delta }_{cs}\le C$ where $s=s_1+s_2$ and $\mathit{\Phi }=\left[{\mathbf{D}}_1\right|{\mathbf{D}}_2]$, and the dictionaries ${\mathbf{D}}_1$ and ${\mathbf{D}}_2$ satisfy a mutual coherence condition $\mu s\le c$.

In this paper, we propose an IHT type algorithm for compressed data separation. We show that under essentially the same conditions as $l_1$-split analysis, the algorithm after a finite number of iterations has the following error bound:

$${\Vert {\mathbf{f}}_1^t-{\mathbf{f}}_1\Vert }_2+{\Vert {\mathbf{f}}_2^t-{\mathbf{f}}_2\Vert }_2\le C_0\left(\frac{{\sigma }_{s_1}\left({\mathbf{D}}_1^{\ast }{\mathbf{f}}_1\right)+{\sigma }_{s_2}\left({\mathbf{D}}_2^{\ast }{\mathbf{f}}_2\right)}{\sqrt{s_1+s_2}}\right)+C_1\Delta ,$$

where

$$\Delta =\left\{\begin{array}{cc}2\sqrt{s_1+s_2}\lambda ,\hfill & \text{if }{\Vert {\mathit{\Phi }}^{\ast }{\mathbf{A}}^{\ast }\mathbf{e}\Vert }_{\infty }\le \lambda ,\hfill \\ \sqrt{2(1+{\delta }_{4s})}ϵ,\hfill & \text{if }{\Vert \mathbf{e}\Vert }_2\le ϵ.\hfill \end{array}\right.$$

As a corollary, our results show that if $\mathbf{e}$ is a Gaussian noise, $\mathbf{e}\sim \mathcal{N}(0,{\sigma }^2{\mathbf{I}}_m)$, then with high probability

$${\Vert {\mathbf{f}}_1^t-{\mathbf{f}}_1\Vert }_2+{\Vert {\mathbf{f}}_2^t-{\mathbf{f}}_2\Vert }_2\le C_0\left(\frac{{\sigma }_{s_1}\left({\mathbf{D}}_1^{\ast }{\mathbf{f}}_1\right)+{\sigma }_{s_2}\left({\mathbf{D}}_2^{\ast }{\mathbf{f}}_2\right)}{\sqrt{s_1+s_2}}\right)+C_1\sigma \sqrt{(s_1+s_2)log(d_1+d_2)}.$$

The derived error bound is adaptive to the unknown level of sparsity of the object components and thus IHT provides a stronger guarantee (when $s_1+s_2$ is small) than those from [25] for the $l_1$-split analysis.

Rewriting the measurement model (1.2) as

$$\mathbf{y}=\mathbf{A}({\mathbf{D}}_1{\mathbf{x}}_1+{\mathbf{D}}_2{\mathbf{x}}_2)+\mathbf{e},$$

one can naturally use standard compressed sensing techniques to first obtain an estimate $[{\hat{\mathbf{x}}}_1;{\hat{\mathbf{x}}}_2]$ of the sparse coefficient vectors via an $l_1$-minimization or standard IHT and then reconstruct the signal’s components by a synthesis operation ${\hat{\mathbf{f}}}_j={\mathbf{D}}_j{\hat{\mathbf{x}}}_j$. Recovery guarantees for such approaches may be achieved under a condition on the coherence of $\left[\mathbf{A}{\mathbf{D}}_1\right|\mathbf{A}{\mathbf{D}}_2]$, or under a condition on the coherence of $\mathbf{A}{\mathbf{D}}_j$ and the mutual coherence between ${\mathbf{D}}_1$ and ${\mathbf{D}}_2$. However, as noted in [8] it is very hard for $\mathbf{A}{\mathbf{D}}_j$ to satisfy the coherence condition when ${\mathbf{D}}_j$ is highly correlated ($j=1,2$). These are different from our setting here, where we impose no incoherence property on the dictionaries ${\mathbf{D}}_j$ themselves. We shall restrict this work to the setting of real valued signals. Notice that as in [20], our results can be extended to the complex valued signals case.

For a vector $\mathbf{v}\in {\mathbb{R}}^d$, ${\Vert \mathbf{v}\Vert }_0$ is the number of nonzero entries of $\mathbf{v}$. For any $q\in (0,\infty )$, denote ${\Vert \mathbf{v}\Vert }_q={\left({\sum }_{j=1}^d{\left|{\mathbf{v}}_j\right|}^q\right)}^{1∕q}$ and ${\Vert \mathbf{v}\Vert }_{\infty }={max}_j\left|{\mathbf{v}}_j\right|$. For $d\in \mathbb{N}$, we write $\left[d\right]$ to mean $\{1,2,\dots ,d\}$. Given an index set $T\subset \left[d\right]$ and a matrix $\mathbf{D}\in {\mathbb{R}}^{n\times d}$, $T^c$ is the complement of $T$ in $\left[d\right]$, ${\mathbf{D}}_T$ (or ${\left[\mathbf{D}\right]}_T$) is the submatrix of $\mathbf{D}$ formed from the columns of $\mathbf{D}$ indexed by $T$.² Write ${\mathbf{D}}^{\ast }$ to mean the conjugate transpose of a matrix $\mathbf{D}$, ${\mathbf{D}}_T^{\ast }$ to mean ${\left({\mathbf{D}}_T\right)}^{\ast }$, and $\Vert \mathbf{D}\Vert$ to mean the spectral norm of $\mathbf{D}$. For a vector $\mathbf{x}\in {\mathbb{R}}^d$, ${\mathbf{x}}_{\left[s\right]}$ denotes the vector consisting of the $s$ largest entries of $\mathbf{x}$ in magnitude. $C>0$ (or $c$, $c_1$) denotes a universal constant that might be different in each occurrence.