An efficient algorithm for compression-based compressed sensing

TL;DR

This paper introduces a new low-complexity iterative algorithm called C-GD that leverages generic compression codes for efficient compressed sensing recovery, bridging the gap between simple structures and advanced compression-based signals.

Contribution

The paper proposes the C-GD algorithm, enabling the use of complex compression codes in compressed sensing recovery with theoretical guarantees and robustness to noise.

Findings

C-GD achieves state-of-the-art imaging performance with JPEG2000.

Theoretical analysis shows convergence and measurement bounds based on rate-distortion.

C-GD is robust to additive white Gaussian noise.

Abstract

Modern image and video compression codes employ elaborate structures existing in such signals to encode them into few number of bits. Compressed sensing recovery algorithms on the other hand use such signals' structures to recover them from few linear observations. Despite the steady progress in the field of compressed sensing, structures that are often used for signal recovery are still much simpler than those employed by state-of-the-art compression codes. The main goal of this paper is to bridge this gap through answering the following question: Can one employ a given compression code to build an efficient (polynomial time) compressed sensing recovery algorithm? In response to this question, the compression-based gradient descent (C-GD) algorithm is proposed. C-GD, which is a low-complexity iterative algorithm, is able to employ a generic compression code for compressed sensing and therefore elevates the scope of structures used in compressed sensing to those used by compression codes. The convergence performance of C-GD and its required number of measurements in terms of the rate-distortion performance of the compression code are theoretically analyzed. It is also shown that C-GD is robust to additive white Gaussian noise. Finally, the presented simulation results show that combining C-GD with commercial image compression codes such as JPEG2000 yields state-of-the-art performance in imaging applications.

Tables

Table 1. TABLE I: PSNR of 128 × 128 128 128 128\times 128 reconstructions with no measurement noise - Sampled by a random Gaussian measurement matrix.

Method	$m/n$	Boat	House	Barbara	Dog	Panda	Snake
BM3D-AMP	30%	29.66	39.71	31.3	21.30	23.90	20.87
	50%	34.19	43.70	33.70	24.35	26.76	22.76
JPEG-GD	30%	23.77	30.61	24.34	18.01	19.35	18.23
	50%	26.46	33.25	27.01	22.78	23.11	19.98
JP2K-GD	30%	30.68	35.22	29.96	22.45	24.00	20.44
	50%	35.28	40.18	34.67	26.35	27.13	23.03

Table 2. TABLE II: PSNR of reconstruction of 128 × 128 128 128 128\times 128 test images with Gaussian measurement noise with various SNR values - sampled by a random Gaussian measurement matrix

		Barbara		Boat		Panda
Method	$m/n$	SNR=10	SNR=30	SNR=10	SNR=30	SNR=10	SNR=30
BM3D-AMP	30%	19.15	28.20	20.12	28.50	13.67	18.82
	50%	21.38	30.16	22.72	33.65	18.44	21.17
JPEG-CG	30%	14.87	22. 71	13.50	22.01	10.50	18.79
	50%	18.44	24.60	19.21	25.36	15.83	22.01
JP2K-CG	30%	16.82	26.23	21.79	28.43	15.63	22.93
	50%	20.78	30.93	24.82	34.13	20.40	25.85

Table 3. TABLE III: PSNR of 512 × 512 512 512 512\times 512 reconstructions with no noise - sampled by a random partial-DCT measurement matrix.

Method	$m/n$	Boat	House	Barbara	Dog	Panda	Snake
NLR-CS	10%	23.06	27.26	20.34	19.53	21.61	18.20
	30%	26.38	30.74	23.67	23.04	25.60	21.80
JPEG-CG	10%	18.38	24.11	16.36	16.30	17.00	15.10
	30%	24.70	30.51	20.37	21.10	22.01	21.63
JP2K-CG	10%	20.75	26.30	18.64	19.74	18.24	18.36
	30%	27.73	38.07	24.89	24.82	25.70	24.37

Table 4. TABLE IV: PSNR of 512 × 512 512 512 512\times 512 reconstructions with Gaussian measurement noise with various SNR value - sampled by a random partial-DCT measurement matrix.

		MRI		Barbara		Snake
Method	$m/n$	SNR=10	SNR=30	SNR=10	SNR=30	SNR=10	SNR=30
NLR-CS	10%	11.66	24.14	12.10	19.83	10.50	18.75
	30%	12.60	26.84	13.32	24.05	11.98	24.82
JPEG-GD	10%	14.34	20.50	15.60	18.60	12.33	15.67
	30%	19.20	24.70	18.17	22.89	14.40	22.37
JP2K-GD	10%	17.33	25.40	16.53	21.65	18.00	23.12
	30%	21.56	35.38	21.82	28.19	21.06	29.30

Equations

$$\mathbf{y}=\mathrm{A}\mathbf{x}+\mathbf{z},$$

$${\Vert X\Vert }_{{\psi }_2}\triangleq \inf \left\{L>0:\mathrm{E}\left[\exp \left(\frac{|X{|}^2}{L^2}\right)\right]\le 2\right\}<\infty .$$

$${\Vert X\Vert }_{{\psi }_1}\triangleq \inf \left\{L>0:\mathrm{E}\left[{\mathrm{e}}^{\frac{|X|}{L}}\right]\le 2\right\}<\infty .$$

$$f:\mathcal{Q}\to \left\{1,2,\dots ,2^r\right\},$$

$$g:\left\{1,2,...,2^r\right\}\to \mathcal{Q}.$$

$$\mathcal{C}$$

$$\delta =\underset{\mathbf{x}\in \mathcal{Q}}{\sup }{\Vert \mathbf{x}-g(f(\mathbf{x}))\Vert }_2.$$

$$\delta (r)=\underset{\mathbf{x}\in \mathcal{Q}}{\sup }{\Vert \mathbf{x}-g_r(f_r(\mathbf{x}))\Vert }_2.$$

$$r(\delta )=\inf \{r:\delta (r)\le \delta \}.$$

$$\hat{\mathbf{x}}=\underset{\mathbf{u}\in {\mathcal{C}}_r}{\mathrm{argmin}}{\Vert \mathbf{y}-\mathrm{A}\mathbf{u}\Vert }_2^2.$$

$$m=\frac{\zeta r}{\log \frac{1}{\mathrm{e}\delta }}.$$

$$P({\Vert \mathbf{x}-\hat{\mathbf{x}}\Vert }_2\ge \theta {\delta }^{1-\frac{1+\nu }{\zeta }})\le {\mathrm{e}}^{-0.8m}+{\mathrm{e}}^{-0.3\nu r},$$

$$\alpha =\lim \underset{\delta \to 0}{\sup }\frac{r(\delta )}{\log (1/\delta )}.$$

$${\mathbf{x}}^{k+1}\leftarrow$$

$${\mathcal{P}}_{{\mathcal{C}}_r}\left(\mathbf{z}\right)=\underset{\mathbf{u}\in {\mathcal{C}}_r}{\mathrm{argmin}}{\Vert \mathbf{u}-\mathbf{z}\Vert }_2^2.$$

$${\mathcal{P}}_{{\mathcal{C}}_r}\left(\mathbf{x}\right)=g_r(f_r(\mathbf{x})).$$

$${\mathbf{x}}^{k+1}\leftarrow$$

$${\mathbf{x}}^{k+1}\leftarrow {\mathcal{P}}_{{\mathcal{C}}_r}\left({\mathbf{x}}^k+\eta {\mathrm{A}}^T\left(\mathbf{y}-\mathrm{A}{\mathbf{x}}^k\right)\right).$$

$$\Vert {\mathbf{x}}^{k+1}-\tilde{\mathbf{x}}{\Vert }_2\le 0.9\Vert {\mathbf{x}}^k-\tilde{\mathbf{x}}{\Vert }_2+2{\left(2+\sqrt{\frac{n}{m}}\right)}^2\delta +\frac{{\sigma }_z}{{\sigma }_a}\sqrt{\frac{8(1+ϵ)r}{m}},$$

$$\frac{1}{\sqrt{n}}\Vert {\mathbf{x}}^{k+1}-\tilde{\mathbf{x}}{\Vert }_2\le \frac{0.9}{\sqrt{n}}\Vert {\mathbf{x}}^k-\tilde{\mathbf{x}}{\Vert }_2+2{\left(2+\sqrt{\frac{n}{m}}\right)}^2\frac{\delta }{\sqrt{n}}+\frac{{\sigma }_z}{{\sigma }_a}\sqrt{\frac{8(1+ϵ)r}{m}}.$$

$$\frac{1}{\sqrt{n}}\Vert {\mathbf{x}}^{k+1}-\tilde{\mathbf{x}}{\Vert }_2\le \frac{0.9}{\sqrt{n}}\Vert {\mathbf{x}}^k-\tilde{\mathbf{x}}{\Vert }_2+2{\left(2+\sqrt{\frac{n}{m}}\right)}^2\frac{\delta }{\sqrt{n}}+\frac{{\sigma }_z}{{\sigma }_a}\sqrt{\frac{8(1+ϵ)r}{m}}+\frac{\xi }{\sqrt{n}}.$$

$$m>\left(\frac{16K^4(1+ϵ)}{{\mu }_0^2{\sigma }_a^4\log \mathrm{e}}\right)r,$$

$$1-2^{-4rϵ}-{\mathrm{e}}^{-\frac{m}{4}}-2^{-2rϵ},$$

$$\Vert {\mathbf{x}}^{k+1}-\tilde{\mathbf{x}}{\Vert }_2\le {\mu }_0\Vert {\mathbf{x}}^k-\tilde{\mathbf{x}}{\Vert }_2+8\left(1+\frac{3Kn}{{\sigma }_a^{2}m}\right)\delta +\frac{9K{\sigma }_z}{{\sigma }_a^2}\sqrt{\frac{r(1+ϵ)}{m}}.$$

$${\Gamma }_k^n\triangleq \{\mathbf{x}\in {\mathcal{B}}_p^n(1):\Vert \mathbf{x}{\Vert }_0\le k\}.$$

$$\frac{1}{\sqrt{n}}\Vert {\mathbf{x}}^{t+1}-\tilde{\mathbf{x}}{\Vert }_2\le \frac{0.9}{\sqrt{n}}\Vert {\mathbf{x}}^t-\tilde{\mathbf{x}}{\Vert }_2+2{\left(2+\sqrt{\frac{n}{m}}\right)}^2{n}^{-1/2-\gamma }+\frac{{\sigma }_z}{{\sigma }_a}\sqrt{\frac{8(1+ϵ)\tilde{r}}{m}},$$

$$\delta \le n^{-\gamma }.$$

$$r\le \sum _{i=0}^k\log \left(\genfrac{}{}{0ex}{}{n}{i}\right)+k(b+1)\le \log n^{k+1}+k(b+1)=(k+1)\log n+k(b+1),$$

$$\underset{t\to \infty }{\lim }\frac{1}{\sqrt{n}}\Vert {\mathbf{x}}^{t+1}-\tilde{\mathbf{x}}{\Vert }_2=O\left({\left(2+\sqrt{\frac{n}{m}}\right)}^2{n}^{-1/2-\gamma }\right).$$

$$\underset{t\to \infty }{\lim }\frac{1}{\sqrt{n}}\Vert {\mathbf{x}}^{t+1}-\tilde{\mathbf{x}}{\Vert }_2=O\left(\frac{n^{\frac{1}{2}-\gamma }}{m}\right).$$

Full text

I Introduction

The main problem of compressed sensing is to recover an unknown target signal ${\boldsymbol{x}}\in\mathbbm{R}^{n}$ from undersampled linear measurements ${\boldsymbol{y}}\in\mathbbm{R}^{m}$,

[TABLE]

where ${\boldsymbol{\rm A}}\in\mathbbm{R}^{m\times n}$ and ${\boldsymbol{z}}\in\mathbbm{R}^{m}$ denote the measurement matrix and the measurement noise, respectively. Since $m<n$, ${\boldsymbol{x}}$ cannot be recovered accurately, unless we leverage some prior information about the signal ${\boldsymbol{x}}$. Such information can mathematically be expressed by assuming that ${\boldsymbol{x}}\in{\cal S}$, where ${\cal S}\subset\mathbbm{R}^{n}$ is a known set. Intuitively, it is expected that, the “smaller” the set ${\cal S}$, the fewer number of measurements $m$ required for recovering ${\boldsymbol{x}}$. In other words, having more information about the target signal ${\boldsymbol{x}}$ limits the feasible set in a way that enables the recovery algorithm to identify a high-quality estimate of ${\boldsymbol{x}}$ using fewer measurements. In the last decade, researchers have explored several different instances of set ${\cal S}$, such as the class of sparse signals or low-rank matrices [1, 2, 3, 4, 5]. Despite the mathematical elegance of these research directions, their impact has not yet met expectations in many application areas. A major barrier is that real-world signals often exhibit far more complex structures than the simple ones that are studied in theoretical researches. In response to this shortcoming, we propose a new complimentary approach to elevate the scope of compressed sensing much beyond sparsity, low-rankness, etc. Our approach hinges upon the following simple assumption: for a given class of signals, there exists an efficient compression algorithm that is able to represent the signals in that class with few number of bits per symbol. In many application areas, such as image and video processing, thanks to the extensive research performed in the past fifty years, such compression algorithms exist and use sophisticated structures shared by signals in a class. Compressive sensing recovery algorithms that take advantage of similar complex structures potentially considerably outperform current algorithms. This raises the following question that we address in this paper: How can one effectively employ a compression algorithm in a compressed sensing recovery method?

In response to this question, we propose an iterative compression-based compressed sensing recovery algorithm, called compression-based gradient descent (C-GD). C-GD, with no extra effort, elevates the scope of the structures used in signal recovery from simple structures, such as sparsity, to the more sophisticated ones that are used in the state-of-the-art compression codes. Our simulation results confirm that combining C-GD with the state-of-the-art image compression codes such as JPEG2000 yields state-of-the-art performance in compressive imaging applications. In addition to its state-of-the-art performance in the imaging applications, C-GD comes with a theoretical framework that derives its convergence rate, analyzes its performance in the presence of noise and other non-idealities in the system, and determines the impact of different parameters of C-GD on its performance.

Using a compression code for the compressed sensing problem was first discussed in [6] and [7]. Both papers studied this problem from a theoretical standpoint, for deterministic signal models and for stationary processes, respectively. In [6], inspired by the Occam’s razor principle, the compressible signal pursuit (CSP) optimization was proposed. We will review this algorithm briefly in Section II-C. The goal of the CSP optimization as a compression-based recovery algorithm is to, among the compressible signals, find the one that gives the lowest measurement error. Theoretically, it can be proved that, asymptotically, as $n$ goes to infinity and the per-symbol reconstruction distortion approaches zero, the CSP optimization achieves the optimal sampling rate, in cases where such optimal rates are already known in the literature. While these results theoretically support the idea of designing compression-based compressed sensing algorithms, they do not provide a recipe of how to achieve this goal in practice. Solving the CSP requires an exhaustive search over all compressible signals, or in other words, all signals in the codebook, whose size is often too large. For instance if we consider the class of images that can be compressed in $1000$ bits only, then the size of the codebook will be $2^{1000}$. Hence, CSP cannot be used in most real-world applications. On the other hand, the C-GD algorithm that is proposed in this paper is inspired by projected gradient descent. Each step of the algorithm involves two operations: i) moving in the direction of the gradient of $\|{\boldsymbol{y}}-{\boldsymbol{\rm A}}{\boldsymbol{x}}\|^{2}$, and ii) projecting on the set of compressible signals. Both steps can be performed efficiently.

There are several other works in the literature that aim at elevating the scope of signal recovery beyond sparsity and low-rankness, and hence, in spirit, have a similar objective as this paper [8, 9, 10]. Here we briefly mention those results and describe how they compare with our work. To use the framework developed in [8] for a class of signals ${\cal S}$, one should construct a set of atoms such that every ${\boldsymbol{x}}\in{\cal S}$ can be represented using a few of those atoms. Finding non-standard atoms for real-world signals, such as images, is a sophisticated task and hence this approach has not found any application in compressive imaging or video recoding yet. The work of [9] is only concerned with independently and identically distributed (i.i.d.) measurement matrices. While in our theoretical results, we also only consider i.i.d. measurement matrices, as shown later in the section on simulation results, C-GD works well with partial Fourier matrices that are used for MRI and radar applications too. On the contrary, D-AMP developed in [9] fails to work with such matrices. The framework developed in [10] still aims to use sparsity. However, it assumes that many of ${n\choose k}$ subspaces of $k$-sparse vectors, do not belong to our signal space. This framework is still only capable of using simple structures and hence does not offer the state-of-the-art performance in applications of compressed sensing, such as imaging. This point is demonstrated in [9].

Finally, there are many papers in the literature that have considered the compressive imaging application and proposed heuristic algorithms to employ sophisticated structures in the signal. (See [11] and the references therein for a complete list of references.) In the simulation results section, we compare the performance of our algorithm with the state-of-the-art heuristic algorithm NLR-CS proposed in [11]. It turns out that C-GD offers the state-of-the-art performance. We should also emphasize that our framework has several other advantages over this line of work: (i) It is general and can be applied to different applications with no extra effort, (ii) it comes with a theoretical framework that shows its robustness to noise and other non-idealities in the measurement process.

The organization of the paper is as follows. Section II reviews some background information. Section III presents the main contributions of the paper, namely the C-GD algorithms and its theoretical analysis. Section IV consider the application of the C-GD algorithm to some standard classes of signals. Section V presents our simulation results, which show that the C-GD algorithm achieves state-of-the-art performance in image compressed sensing. Section VI provides the proofs of the mains results of the paper. Finally, Section VII concludes the paper.

II Background

In this section, we first review the notation used throughout the paper and also the definitions of Gaussian and sub-Gaussian random variables. Then we review the rate-distortion function of a compression code, and formally define the objective of compression-based signal recovery. Finally we review the CSP optimization, which is the first proposed method for employing compression algorithms for compressed sensing.

II-A Notations and defenitions

Scalar values are denotes by lower-case letters such as $x$. A column vector is denotes by bold letters as ${\boldsymbol{x}}$. The $i$-th element of ${\boldsymbol{x}}$ is denoted by ${x}_{i}$. Given vector ${\boldsymbol{x}}\in\mathds{R}^{n}$, $\|{\boldsymbol{x}}\|_{\infty}=\max_{i}|x_{i}|$ and $\|{\boldsymbol{x}}\|_{2}=\sqrt{\sum_{i=1}^{n}x_{i}^{2}}$ denote the $\ell_{\infty}$ norm and the $\ell_{2}$ norm of ${\boldsymbol{x}}$, respectively. A matrix is denoted by bold capital letters such ${\boldsymbol{\rm X}}$ and its $(i,j)$-th element by ${X}_{i,j}$. The transpose of ${\boldsymbol{\rm X}}$ is given by ${\boldsymbol{\rm X}}^{T}$. $\sigma_{\max}({\boldsymbol{\rm A}})$ denotes the maximum singular value of ${\boldsymbol{\rm A}}$. Calligraphic letters such as $\mathcal{D}$ and ${\cal C}$ denote sets. The size of a set ${\cal C}$ is denotes as $|{\cal C}|$. The unit sphere in $\mathds{R}^{n}$ is denoted by $S^{n-1}$, i.e., $S^{n-1}\triangleq\{{\boldsymbol{u}}\in\mathbb{R}^{n}\ |\ \|{\boldsymbol{u}}\|_{2}=1\}$. Throughout the paper $\log$ and $\ln$ refer to logarithm in base $2$ and natural logarithm, respectively. Finally, we use $O$ and $\Omega$ notation, defined below, to describe the limiting behavior of certain quantities. $f(n)=O(g(n))$ as $n\rightarrow\infty$, if and only if there exist $n_{0}$ and $c$ such that for any $n>n_{0}$, $|f(n)|\leq c|g(n)|$. Likewise, $f(n)=\Omega(g(n))$ as $n\rightarrow\infty$, if and only if there exist $n_{0}$ and $c$ such that for any $n>n_{0}$, $|f(n)|\geq c|g(n)|$.

In this paper, the main results are proved for both Gaussian and sub-Gaussian measurement matrices. In the following, we briefly review the definition of sub-Gaussian and sub-exponential random variables.

Definition 1** (Sub-Gaussian).**

A random variable $X$ is sub-Gaussian when

[TABLE]

Note that Gaussian random variables are also sub-Gaussian random variables. For $X\sim\mathcal{N}(0,\sigma^{2})$, ${\left\|{X}\right\|}_{\psi_{2}}{}=\sqrt{8\over 3}\sigma$.

Definition 2** (Sub-exponential).**

A random variable $X$ is a sub-exponential random variable when

[TABLE]

Using the above definitions, it is straightforward to show the following result.

Lemma 1**.**

Let $X$ and $Y$ be sub-Gaussian random variables. Then $XY$ is sub-exponential. Moreover, ${\left\|{XY}\right\|}_{\psi_{1}}\leq{\left\|{X}\right\|}_{\psi_{2}}{\left\|{Y}\right\|}_{\psi_{2}}.$

II-B Objective of compression-based compressed sensing

Consider a compact set ${\cal Q}\subset\mathbbm{R}^{n}$. (Throughout the paper, we focus on this deterministic signal model. However, all the results can be extended to the stochastic setting as well. Refer to [12] for further information.) A rate-$r$ lossy compression code for set ${\cal Q}$ is characterized by its encoding and decoding mappings $(f,g)$, where

[TABLE]

and

[TABLE]

Encoding and decoding mappings $(f,g)$ define a codebook $\mathcal{C}$, where

[TABLE]

Clearly, $|\mathcal{C}|\leq 2^{r}$. The performance of a code defined by $(f,g)$ is characterized by its i) rate $r$ and ii) maximum distortion $\delta$ defined as

[TABLE]

The problem of compression-based compressed sensing can be formally stated in the following way. Suppose that a family of compression codes $(f_{r},g_{r})$, parameterized with rate $r$ is given for ${\cal Q}$. For instance, JPEG or JPEG2000 compression algorithms at different rates can be considered as a family of compression algorithms for the class of natural images. The deterministic distortion-rate function of this family is given by

[TABLE]

Note that for any reasonable code, $\delta(r)$ is expected to be a monotonically non-increasing function of $r$. We also define the deterministic rate-distortion function of the compression code as

[TABLE]

Remark 1*.*

Typically a family of compression codes is defined as a sequence of compression codes that are indexed by blocklength $n$ and operate either at a fixed rate or a fixed distortion. In this paper, we are more interested in the setting where the blocklength is fixed and the rate or the distortion changes. Therefore, we consider a family of compression codes with fixed blocklength $n$ indexed by rate $r$.

Based on these definitions we can formally state the objectives of this paper as follows. Consider the problem of compressed sensing of signals in a set ${\cal Q}$. Suppose that instead of knowing the set ${\cal Q}$ explicitly, for signals in ${\cal Q}$, we have access to a family of compression algorithms with rate-distortion function $r(\delta)$. For ${\boldsymbol{x}}\in{\cal Q}$ our goal is not to compress it, but to recover it from its undersampled set of linear projections ${\boldsymbol{y}}={\boldsymbol{\rm A}}{\boldsymbol{x}}+{\boldsymbol{z}}$. The goal of compression-based compressed sensing is summarized in the following two questions:

Question 1: Can one employ a given compression code in an efficient (polynomial time) signal recovery algorithm?

Question 2: Can we characterize the number of observations such an algorithm requires to accurately recover ${\boldsymbol{x}}$, in terms of the rate-distortion performance of the code?

Note that the algorithms that are developed in response to the above questions will automatically employ the structure captured by the compression algorithms, and hence can immediately elevate the scope of compressed sensing much beyond simple structures. For instance, an MPEG4-based compressed sensing recovery algorithm will use not only the intra-frame but also the inter-frame dependencies among different pixels.

II-C Compressible signal pursuit

Consider the following simplified version of our main questions: Can one employ a compression code for signal recovery? Can we characterize the number of observations such algorithm requires for an accurate recovery of ${\boldsymbol{x}}$ in terms of the rate-distortion performance of the code? Note that the only simplification is that we have removed the constraint on the computational complexity of the recovery scheme. In response to this simplified question, [6] proposed the compressible signal pursuit (CSP) optimization that estimates signal ${\boldsymbol{x}}$ based on measurements ${\boldsymbol{y}}$ as follows. Among all the signals ${\boldsymbol{u}}\in{\cal Q}$ that satisfy the measurement constraint, i.e.  ${\boldsymbol{y}}={\boldsymbol{\rm A}}{\boldsymbol{u}}$, it searches for the one that can be compressed well by the compression code described by $(f_{r},g_{r})$. More formally, given a lossy compression code with codebook $\mathcal{C}_{r}=\{g_{r}(f_{r}({\boldsymbol{x}}):{\boldsymbol{x}}\in{\cal Q}\}$, the CSP optimization recovers ${\boldsymbol{x}}$ from its measurements ${\boldsymbol{y}}={\boldsymbol{\rm A}}{\boldsymbol{x}}$ as follows

[TABLE]

The performance of the CSP optimization is characterized in [6] and [7], under deterministic and stochastic signal models, respectively. Before we mention the theoretical results, we should emphasize that at this point CSP is based on an exhaustive search over the codebook and hence is computationally infeasible. The following result from [6] characterizes the performance of the CSP optimization in the noiseless setting, where ${\boldsymbol{z}}={\boldsymbol{0}}$.

Theorem 1** (Corollary 1 in [6]).**

Consider a family of compression codes $(f_{r},g_{r})$ for set ${\cal Q}$ with corresponding codebook ${\cal C}_{r}$ and rate-distortion function $r(\delta)$. Let ${\boldsymbol{\rm A}}\in\mathbbm{R}^{m\times n}$, where $A_{i,j}$ are i.i.d. ${\cal N}(0,1)$. For ${\boldsymbol{x}}\in{\cal Q}$ and ${\boldsymbol{y}}={\boldsymbol{\rm A}}{\boldsymbol{x}}$, let $\hat{{\boldsymbol{x}}}$ denote the solution of (3). Given $\nu>0$ and $\zeta>1$, such that ${\zeta\over\log{1\over{\rm e}\delta}}<\nu$, let

[TABLE]

Then,

[TABLE]

where $\theta=2{\rm e}^{-(1+\nu)/\zeta}$.

For a simpler interpretation of this result define the $\alpha$-dimension of a family of compression codes with rate-distortion function $r(\delta)$ as

[TABLE]

Suppose that a small value of $\delta$ (or large value of $r$) is used in CSP. Then, roughly speaking, Theorem 1 implies that CSP returns an almost accurate estimate of ${\boldsymbol{x}}$ as long as $m>\alpha$. Note that $\alpha$ is usually much smaller than $n$, and hence the number of measurements CSP requires is much smaller than the ambient dimension of the signal. (See Section IV for some classical examples.)

Remark 2*.*

In this paper, we focus entirely on the deterministic setting. The performance of CSP in the stochastic setting is studied in [12]. Using the connection between the Rényi information dimension and rate distortion dimension of a random variable [13], it has been proved in [12, 14] that for i.i.d. sources with a mixture of discrete and continuous distribution, CSP achieves the optimal sampling rate.

Remark 3*.*

The robustness of CSP to deterministic and stochastic measurement noises has also been proved in [6]. However, for the sake of brevity we do not repeat those results here and only discussed the noiseless setting in Theorem 1.

Unfortunately, the positive theoretical properties of CSP are overshadowed by the fact that implementing it requires an exhaustive search over the set of all codewords. This makes CSP computationally infeasible for all real-world applications. In the next section, we propose an efficient CS recovery algorithm that employs compression code and compare its performance with that of CSP.

III Our main contributions

III-A Compression-based gradient descent (C-GD)

As discussed in the last section, CSP is based on an exhaustive search and is computationally infeasible for real-world signals. In response to this drawback of CSP, we propose a computationally efficient and theoretically analyzable approach to approximate the solution of CSP. Towards this goal, inspired by the projected gradient descent (PGD) algorithm [15], we propose the following iterative algorithm: Start from some ${\boldsymbol{x}}^{0}\in\mathds{R}^{n}$. For $k=1,2,\ldots$,

[TABLE]

where

[TABLE]

Here index $k$ denotes the iteration number and $\eta_{k}\in\mathbbm{R}$ denotes the step size. We refer to this algorithm as compression-based gradient descent (C-GD). Each iteration of this algorithm involves performing two operations. In the first step, it moves in the direction of the negative of $\|{\boldsymbol{y}}-{\boldsymbol{\rm A}}{\boldsymbol{x}}\|_{2}^{2}$ to find solutions that are closer to the ${\boldsymbol{y}}={\boldsymbol{\rm A}}{\boldsymbol{u}}$ hyperplane. The second step, i.e., the projection step, ensures that the estimate C-GD obtains belongs to the codebook.

The first step of C-GD is straightforward and requires two matrix-vector multiplications. For the second step, ideally, applying the encoder and decoder to a signal ${\boldsymbol{x}}$ yields the closest codeword of the compression code. Hence, we make the following assumption about the compression code to ${\boldsymbol{x}}$.

Assumption 1**.**

In analyzing the performance of C-GD, we assume that the compression code $(f_{r},g_{r})$ satisfies

[TABLE]

Under Assumption 1, the projection step of C-GD can be implemented efficiently. More precisely, under this assumption, the C-GD algorithm is simplified to

[TABLE]

Hence, each step of this algorithm requires two matrix-vector multiplication and an application of the encoder and the decoder of the given compression code. In the next section, we summarize our theoretical results regarding the performance of C-GD. Note that for the notational simplicity, we present all our results under Assumption 1. However as will be discussed after Corollary 1, we can relax this assumption and still analyze the iterative algorithm proposed in (7).

III-B *Convergence Analysis of C-GD *

The objective of this section is to theoretically analyze some of the properties of C-GD. As discussed before, the measurement vector is denoted with ${\boldsymbol{y}}={\boldsymbol{\rm A}}{\boldsymbol{x}}+{\boldsymbol{z}}$, where ${\boldsymbol{x}}\in{\cal Q}$, ${\boldsymbol{\rm A}}\in\mathbbm{R}^{m\times n}$, and ${\boldsymbol{z}}$ is the noise. Furthermore, we assume that a family of compressions codes $(f_{r},g_{r})$ parameterized with the rate $r$ that is known for ${\cal Q}$. Starting with ${\boldsymbol{x}}^{0}$, C-GD uses iterations

[TABLE]

to obtain a good estimate of ${\boldsymbol{x}}$. In our theoretical analysis of C-GD, we focus on popular measurement matrices in compressed sensing area namely dense i.i.d. Gaussian and sub-Gaussian matrices. In Section V, we present our simulation results that confirm the success of C-GD for Fourier matrices as well. However, the theoretical study of this important class of matrices is left for future research. We first study the performance of C-GD for i.i.d. Gaussian measurement matrices. Then, we extend our results to i.i.d. sub-Gaussian measurement matrices.

Theorem 2**.**

Let ${\boldsymbol{\rm A}}\in\mathbbm{R}^{m\times n}$ be a random Gaussian measurement matrix with i.i.d entries ${A}_{i,j}\sim\mathcal{N}\left(0,\sigma_{a}^{2}\right)$, and ${\boldsymbol{z}}\in\mathbbm{R}^{m}$ be an i.i.d. Gaussian noise vector with $z_{i}\sim\mathcal{N}\left({\boldsymbol{0}},\sigma_{z}^{2}\right)$. Let $\eta=\frac{1}{m\sigma_{a}^{2}}$ and define $\tilde{{\boldsymbol{x}}}=\mathcal{P}_{\mathcal{C}_{r}}({\boldsymbol{x}})$, where $\mathcal{P}_{\mathcal{C}_{r}}(\cdot)$ is defined in (6). If Assumption 1 holds, then, given $\epsilon>0$, for $m\geq 80r\left(1+\epsilon\right),$ with a probability larger than $1-2^{-2\epsilon r+1}$, we have

[TABLE]

for $k=0,1,2,\ldots$.

The proof of Theorem 2 is given in Section VI-B.

A few important features of this theorem are discussed in the following remarks.

Remark 4*.*

According to Theorem 2, C-GD requires $\Omega(r(\delta))$ measurements (for small values of $\delta$) to obtain an accurate estimate. Hence, according to this theorem, even in the noiseless setting we should not let $\delta\rightarrow 0$. Otherwise, for many signal classes $r(\delta)\rightarrow\infty$ and hence C-GD will require more measurements than the ambient dimension of ${\boldsymbol{x}}_{o}$. As we will demonstrate in several examples in Section IV, one can set $\delta$ to a small dimension-dependent value, e.g. $\delta=1/n$, to ensure that C-GD can obtain a very good estimate with a few observations. Section IV studies how $\delta$ is set and connects Theorem 2 to some classical results in compressed sensing.

Remark 5*.*

According to Theorem 1, CSP requires \Omega\Big{(}\frac{r(\delta)}{\log(1/\delta)}\Big{)}. This implies that, in the noiseless setting, for fixed $n$ and $m$, the estimate of CSP improves as $\delta$ decreases. However, this seems not to be the case for C-GD. According to Theorem 2, C-GD requires more than $\Omega(r(\delta))$. Hence, as $\delta$ decreases, C-GD requires more measurements. This is not an issue in almost all the applications of compressed sensing, where $n$ is very large and hence there is not much difference between setting $\delta=1/n$ and $\delta=0$. However, from a theoretical perspective it is interesting to discover whether this mismatch is an artifact of our proof technique or it is a fundamental loss that is incurred by the reduction in the computational complexity. This question is left for future research.

We postpone the discussion of the convergence rate and the reconstruction error to Section IV, where we discuss some classical examples and compare the conclusion of this theorem with some classical results in compressed sensing. In the setup considered in Theorem 2, as $n$ increases the per measurement’s signal to noise ratio (SNR) is approaching infinity. Note that by scaling the measurement matrix, we can also obtain results for fixed SNR. The next corollary clarifies our claim.

Corollary 1**.**

Consider the setup of Theorem 2, where now ${\boldsymbol{\rm A}}\in\mathbbm{R}^{m\times n}$ is a random Gaussian measurement matrix with i.i.d entries ${A}_{i,j}\sim\mathcal{N}\left(0,\frac{\sigma_{a}^{2}}{n}\right)$ and $\eta={n\over\sigma_{a}^{2}m}$. If Assumption 1 holds, then, given $\epsilon>0$, for $m\geq 80r\left(1+\epsilon\right)$, with a probability larger than $1-2^{-2\epsilon r+1}$, for $k=0,1,2,\ldots$, we have

[TABLE]

Finally, Assumption 1 seems to play a critical role in Theorem 2 and Corollary 1. However, thanks to the linear convergence of C-GD one can relax Assumption 1 in several ways and still obtain recovery guarantees for C-GD. For the sake of brevity we only mention one such result in the paper.

Theorem 3**.**

Consider the setup of Theorem 2, where now ${\boldsymbol{\rm A}}\in\mathbbm{R}^{m\times n}$ is a random Gaussian measurement matrix with i.i.d entries ${A}_{i,j}\sim\mathcal{N}\left(0,\frac{\sigma_{a}^{2}}{n}\right)$ and $\eta={n\over\sigma_{a}^{2}m}$. Suppose that $\sup_{{\boldsymbol{x}}}\|g_{r}(f_{r}({\boldsymbol{x}}))-\mathcal{P}_{\mathcal{C}_{r}}({\boldsymbol{x}})\|_{2}\leq\xi$. Then, given $\epsilon>0$, for $m\geq 80r\left(1+\epsilon\right)$, with a probability larger than $1-2^{-2\epsilon r+1}$, for $k=0,1,2,\ldots$, we have

[TABLE]

The proof can be found in Section VI-C. Note that at every iteration the imperfect projection introduces an error. These errors accumulate as the algorithm proceeds. However, thanks to the linear convergence of the algorithm the overall error caused by the imperfect projection remains at the order of $O(\xi/\sqrt{n})$.

All our results so far have been stated for Gaussian measurement matrices. However, they can be generalized to subgaussian matrices too. To prove our claim we extend one of our results, i.e. Theorem 2, to sub-Gaussian matrices below. Our next theorem shows that Theorem 2 can be extended to sub-Gaussian measurement matrices.

Theorem 4**.**

Let ${\boldsymbol{\rm A}}\in\mathbbm{R}^{m\times n}$ be a zero-mean random sub-Gaussian measurement matrix with i.i.d entries, such that ${\left\|{{A}_{i,j}}\right\|}_{\psi_{2}}\leq K$ and ${\rm E}\left[{A}_{i,j}^{2}\right]=\sigma_{a}^{2}$. The noise vector ${\boldsymbol{z}}$ is distributed as $\mathcal{N}\left({\boldsymbol{0}},\sigma_{z}^{2}{\boldsymbol{\rm I}}_{m\times m}\right)$. Set $\eta=\frac{1}{m\sigma_{a}^{2}}$. Then, given $\epsilon>0$ and $\mu_{0}\in(0,1)$, such that $\mu_{0}\sigma_{a}^{2}\leq 2K^{2}$, for

[TABLE]

with probability at least

[TABLE]

for $k=0,1,2,\ldots$,

[TABLE]

Proof is given in Section VI-D. Since all the terms in Theorem 4 are similar to the corresponding terms in Theorem 2 we do not discuss them here. We only discuss the convergence rate. Note that the convergence rate $\mu_{0}$ has a direct impact on the number of measurements. If we want the convergence to be fast ($\mu_{0}$ to be small), we should either increase the number of measurements or decrease the rate $r$. If we decrease the rate, the two error terms $8\left(1+\frac{3Kn}{\sigma_{a}^{2}m}\right)\delta+{9K\sigma_{z}\over\sigma_{a}^{2}}\sqrt{r(1+\epsilon)\over m}$ grow.

IV Standard signal classes

In this section we discuss the corollaries of our main theorem for the following two standard signal classes that have been studied extensively in the literature: (i) sparse signals, and (ii) piecewise polynomials. For each class, we first construct a simple compression algorithm that can be efficiently implemented in practice, and then explain the implications of C-GD and its analysis for that class. These examples enable us to shed light on different aspects of C-GD, such as (i) convergence rate, (ii) number of measurements, (iii) reconstruction error in noiseless setting, and (iv) reconstruction error in the presence of noise.

IV-A Sparse signals

Let $\mathcal{B}_{p}^{n}(\rho)\triangleq\{{\boldsymbol{x}}\in\mathds{R}^{n}:\ \|{\boldsymbol{x}}\|_{p}\leq\rho\}$ represent a ball of radius $\rho$ in $\mathds{R}^{n}$. Also, let $\Gamma_{k}^{n}$ denote the set of all $k$-sparse signals in $\mathcal{B}_{p}^{n}(\rho)$, i.e.,

[TABLE]

In order to apply C-GD to this class of signals, we first need to construct a family of compression codes for such sparse bounded signals. Consider a family of compression codes for set $\Gamma_{k}^{n}$ defined as follows. For ${\boldsymbol{x}}\in\Gamma_{k}^{n}$, (i) encode the locations of its at most $k$ non-zero entries ($\approx\log{n\choose k})$ bits) and (ii) apply a uniform quantizer to the magnitudes of the non-zero entries. (Using $b$ bits for the magnitude of each entry and one bit for its sign, this step spends $(b+1)k$ bits.) Using this specific compression algorithm in the C-GD framework yields an algorithm which is very similar to iterative hard thresholding (IHT) [16]. At every iteration, after moving in the opposite direction of the gradient, the standard IHT algorithm keeps the $k$ largest elements and sets the rest to zero. The C-GD algorithm on the other hand, while having the same first step, performs the projection step in a slightly different manner. For the projection onto codewords, similar to IHT, it first finds the $k$ largest entries. Then, for each entry $x_{i}$, it first limits it between $[-1,1]$ by computing $x_{i}{\boldsymbol{1}}_{x_{i}\in(-1,1)}+{\boldsymbol{1}}_{x_{i}\geq 1}-{\boldsymbol{1}}_{x_{i}\leq-1}$. Then, it quantizes the result by $b+1$ bits. The following corollary enables us to compare our results with that of hard thresholding.

Consider ${\boldsymbol{x}}\in\Gamma_{k}^{n}$ and let ${\boldsymbol{y}}={\boldsymbol{\rm A}}{{\boldsymbol{x}}}+{\boldsymbol{z}}$, where ${A}_{i,j}\stackrel{{\scriptstyle\rm i.i.d.}}{{\sim}}\mathcal{N}\left(0,{\sigma_{a}^{2}}/{n}\right)$ and ${z}_{i}\stackrel{{\scriptstyle\rm i.i.d.}}{{\sim}}\mathcal{N}\left(0,{\sigma_{z}^{2}}\right)$. Let $\tilde{{\boldsymbol{x}}}$ denote the projection of ${\boldsymbol{x}}$ onto the codebook of the above-described code. The following corollary of Theorem 2 characterizes the convergence performance of C-GD applied to ${\boldsymbol{y}}$ when using this code.

Corollary 2**.**

Given $\gamma>0$, set the quantization level of the compression code to $b+1=\lceil\gamma\log n+{1\over 2}\log k\rceil+1$ bits. Also, set $\eta={n\over\sigma_{a}^{2}m}$. Then, given $\epsilon>0$, for $m\geq 80\tilde{r}(1+\epsilon)$, where $\tilde{r}=(1+\gamma)k\log n+{k\over 2}\log k+2k$,

[TABLE]

for $t=1,2,\ldots$, with probability larger than $1-2^{-2\epsilon\tilde{r}}$.

Proof.

Consider $u\in[-1,1]$. Quantizing $u$ by a uniform quantizer that uses $b+1$ bits yields $\hat{u}$, which satisfies $|u-\hat{u}|<2^{-b}$. Therefore, using $b+1$ bits to quantize each non-zero element of ${\boldsymbol{x}}\in\Gamma_{k}^{n}$ yields a code which achieves distortion $\delta\leq 2^{-b}\sqrt{k}$. Hence, for $b+1=\lceil\gamma\log n+{1\over 2}\log k\rceil+1$,

[TABLE]

On the other hand, the code rate $r$ can be upper-bounded as

[TABLE]

where the last inequality holds for all $n$ large enough. The rest of the proof follows directly from inserting these numbers in the statement of Theorem 2. ∎

This corollary enables us to provide further intuition on the performance of C-GD. We start with the noiseless observations and for the moment we only study the required number of measurements and the reconstruction error in the absence of noise. The number of measurements required by the C-GD algorithm is $m=\Omega(k\log n)$. For the final reconstruction error, we can use (13) and obtain

[TABLE]

This implies that the recovery error satisfies

[TABLE]

Hence, if $\gamma>0.5$, the error vanishes as the dimension grows. Regarding the number of measurements, there are two points that we would like to emphasize here:

The $\alpha$-dimension of the above-described code is $k$. Hence, CSP is able to accurately recover the signal from only $k$ measurements. However, solving CSP requires an exhaustive search over all the codebooks. On the other hand, C-GD requires $k\log n$ measurements. Therefore, it seems that the extra $\log n$ factor is the price for having an efficient recovery algorithm. 2. 2.

For large values of $n$, $n^{-\gamma}$ is very small, and hence C-GD becomes very similar to IHT. The results we have obtained for C-GD in this case are slightly weaker than those provided for IHT. First, our reconstruction is not exact even in the noiseless setting. Second, the number of measurements C-GD requires is $O(k\log n)$ compared to $O(k\log(n/k))$ required by IHT. These minor differences seem to be the price of the generality of the bounds derived for C-GD.

So far, we have studied two important quantities in Corollary 1, i.e., (i) required number of measurements, and (ii) reconstruction error in the absence of noise. The last important term, is the reconstruction error in the presence of the noise. From Corollary 2, the distortion caused due to the presence of a Gaussian noise is $O\left(\frac{\sigma_{z}}{\sigma_{a}}\sqrt{\frac{\tilde{r}}{m}}\right)$, or $O\left(\frac{\sigma_{z}}{\sigma_{a}}\sqrt{\frac{k\log n}{m}}\right)$. Note that there is no result on the performance of IHT in the presence of stochastic measurement noise. However, this noise sensitivity is comparable with the performance of algorithms that are based on convex optimization such as LASSO and Dantzig selector [17, 18].

IV-B Piecewise polynomial functions

Let ${\rm Poly}_{N}^{Q}$ denote the class of piecewise-polynomial functions $p(\cdot):[0,1]\rightarrow[0,1]$ with at most $Q$ singularities111A singularity is a point at which the function is not infinitely differentiable., where each polynomial has a maximum degree of $N$. For $p\in{\rm Poly}_{N}^{Q}$, let $(x_{1},x_{2},\ldots,x_{n})$ be the samples of $p$ at

[TABLE]

For $\ell=1,\ldots,Q$, let $\{a_{i}^{\ell}\}_{i=0}^{N_{\ell}}$ denote the set of coefficients of the $\ell^{\rm th}$ polynomial in $p$. Here, $N_{\ell}\leq N$ denotes the degree of the $\ell$-th polynomial. For notational simplicity, assume that the coefficients of each polynomial belong to $[0,1]$ interval and that $\sum_{i=0}^{N_{\ell}}a^{\ell}_{i}<1$, for every $\ell$. Define

[TABLE]

Note that this class of functions is a generalization of the class of piecewise-constant functions that are popular in many applications including imaging. To apply C-GD to this class of signals, we need to design an efficient compression code for signals in $\mathcal{P}$ and describe how to project signals on the codewords of $\cal{P}$. For the first part, consider a simple code which for any signal ${\boldsymbol{x}}\in\cal{P}$, it first describes the locations of its discontinuities and then, using a uniform quantizer that spends $b$ bits per coefficient, describes the quantized coefficients of the polynomials. For the other task, which is projecting a signal ${\boldsymbol{x}}\in\mathds{R}^{n}$, Appendix A describes how we can find the closest signal to ${\boldsymbol{x}}$ in $\mathcal{P}$, i.e.,

[TABLE]

using dynamic programing. Once that signal is found, its quantized version using the described code represents the desired projection.

Note that C-GD combined with the described compression code is an extension of IHT to piecewise-polynomial functions. At every iteration, C-GD projects its current estimate of the signal to space of piecewise-polynomial functions.

Consider ${\boldsymbol{x}}\in\mathcal{P}$ and ${\boldsymbol{y}}={\boldsymbol{\rm A}}{\boldsymbol{x}}+{\boldsymbol{z}}$, where ${A}_{i,j}\stackrel{{\scriptstyle\rm i.i.d.}}{{\sim}}\mathcal{N}\left(0,{\sigma_{a}^{2}}/{n}\right)$ and ${z}_{i}\stackrel{{\scriptstyle\rm i.i.d.}}{{\sim}}\mathcal{N}\left(0,{\sigma_{z}^{2}}\right)$. Similar to Corollary 2, the following corollary characterizes the convergence performance of C-GD combined with the described compression code, when applied to measurements ${\boldsymbol{y}}$.

Corollary 3**.**

Set the step size in C-GD as $\eta={n\over\sigma_{a}^{2}m}$ and the quantization level in the compression code as

[TABLE]

where $\gamma>0$ is given. Set $\tilde{r}=\left((\gamma+0.5)(N+1)(Q+1)+Q\right)\log n+(N+1)(Q+1)(\log(N+1)+1)+1$. Then, given $\epsilon>0$, for $m\geq 80\tilde{r}\left(1+\epsilon\right)$, for $t=0,1,2,\ldots$, we have

[TABLE]

with a probability larger than $1-2^{-2\epsilon\tilde{r}+1}$.

Proof.

To apply Theorem 2, we need to find the rate-distortion performance of the described compression code. Using $b$ bits per coefficient, the described code, in total, spends at most $r$ bits, where

[TABLE]

For $b=\lceil(\gamma+0.5)\log n+\log(N+1)\rceil$,

[TABLE]

On the other hand, using $b$ bits per coefficient, the distortion in approximating each point can be bounded as

[TABLE]

where $[a^{\ell}_{i}]_{b}$ denotes the $b$-bit quantized version of $a^{\ell}_{i}$. Therefore, the overall error is bounded as

[TABLE]

Choosing $b=\lceil(\gamma+0.5)\log n+\log(N+1)\rceil$, as prescribed by the corollary, ensures that

[TABLE]

Inserting these numbers in Theorem 2 yields the desired result.

∎

The important quantities in the above corollary are explained in the following.

Required number of measurements: If we assume that $n$ is much larger than $N$ and $Q$, then the required number of measurements is $\Omega((N+1)(Q+1)\log n)$. Note that given the degrees of freedom of piecewise polynomial functions, we do not expect to be able to recovery ${\boldsymbol{x}}\in{\cal P}$ with fewer than $(N+1)(Q+1)$ observations. 2. 2.

Reconstruction error in the absence of measurement noise: Similar to the discussion of the previous section we can argue that

[TABLE]

Hence, the error goes to zero for every $\gamma>0.5$. 3. 3.

Reconstruction error in the presence of measurement noise: In this case, the impact of Gaussian noise in the upper bound is $O\left(\frac{\sigma_{z}}{\sigma_{a}}\sqrt{\frac{(Q+1)(N+1)\log n}{m}}\right)$.

V Simulation Results and Discussion

In this section, we assess the performance of the C-GD algorithm in various settings. Furthermore, we compare our results with the state-of-the-art recovery algorithms, such as denoising-based approximate message passing (D-AMP) [9] and nonlocal low-rank regularization (NLR-CS) [11]. Throughout this section, when compression algorithm $\mathcal{X}$ is used in the platform of C-GD, the resulting algorithm is referred to $\mathcal{X}$-GD. For instance, the recovery algorithm that employs the JPEG code is called JPEG-GD.

V-A Parameters setting

Running C-GD involves specifying three free parameters: (i) step-size $\eta$, (ii) compression rate $r$, and (iii) the number of iterations we run the algorithm. The success of C-GD relies on proper tuning of the first two parameters, i.e., step-size and compression rate. In this section, we explain how we tune these parameters in practice. In Section III-B, we theoretically showed that the algorithm converges to the optimal solution for $\eta=\frac{1}{m\sigma_{a}}.$ However, this choice of step size might yield a very slow convergence in practice. Hence, in our simulations we follow an adaptive strategy for setting $\eta$. Let $\eta_{k}$ denote the step size at the $k^{\rm th}$ iteration. Then, we set $\eta_{k}$ to

[TABLE]

where, for ${\boldsymbol{u}}\in\mathds{R}^{n}$, $l\left({\boldsymbol{u}}\right)\triangleq{\left\|{{\boldsymbol{y}}-{\boldsymbol{\rm A}}{\boldsymbol{u}}}\right\|}_{2}$. In other words, $\eta_{k}$ is set such that the next estimate is moved as close as possible to the subspace $\mathcal{V}=\{{\boldsymbol{u}}\ |\ {\boldsymbol{y}}={\boldsymbol{\rm A}}{\boldsymbol{u}}\}$. Note that regardless of the value of $\eta$, the next estimate will be a codeword. Hence, intuitively speaking, the closer this codeword is to the subspace $\mathcal{V}$, a better estimate ${\boldsymbol{x}}^{k+1}$ will be. The optimization problem proposed in (18) is a simple scalar optimization problem, and we use derivative-free methods, such as Nelder-Mead method (or downhill simplex) method [19], to solve it. In our simulations we noticed that this strategy speeds up the convergence rate of C-GD.

Finding the optimal choice of $r$ is an instance of the model selection problem in statistics and machine learning. (See Chapter 7 of [20].) Hence, standard techniques such as multi-fold cross validation can be used. Note that multi-fold cross validation increases the computational complexity of our recovery algorithm. Reducing the computational complexity of such model selection techniques is left for future research.

To control the number of iterations in the C-GD method, we consider two standard stopping rules. One is to limit the maximum number of iterations, which is defined as the parameter $K_{1,\max}$ in Algorithm 1. The second stopping rule is a predefined threshold on the reduction of squared-error in each iteration, i.e., ${\left\|{{\boldsymbol{x}}^{k+1}-{\boldsymbol{x}}^{k}}\right\|}_{2}$. In our numerical simulation we set $K_{1,\max}=50$ and $\varepsilon_{T}=0.001$. The specific algorithm that is employed in our simulations is presented below in Algorithm 1.

V-B Algorithms and comparison criteria

We explore the performance of our C-GD algorithm for the compressive imaging application. We employ the standard image compression algorithms JPEG and JPEG2000 in our C-GD framework and obtain JPEG-GD and JP2K-GD recovery schemes. In our numerical simulation, we use the implementation of JPEG2000 and JPEG codecs in the Matlab-R2016b Image and Video processing package.

We compare the performance of our algorithm on six standard test images that are shown in Fig. 1. To quantitatively evaluate the quality of an estimated image, we use the peak signal-to-noise ratio (PSNR) defined as

[TABLE]

where for a noise-free grayscale image ${\boldsymbol{x}}$ and its reconstructed image $\hat{{\boldsymbol{x}}}$, the mean square error (MSE) is defined as ${\rm MSE}=\frac{1}{n}{\left\|{{\boldsymbol{x}}-\hat{{\boldsymbol{x}}}}\right\|}_{2}^{2}$. Furthermore, we define the signal-to-noise ratio (${\rm SNR}$) in the measurement vector ${\boldsymbol{y}}$ as

[TABLE]

Even though our theoretical results consider i.i.d. Gaussian matrices, we evaluate the performance of our algorithm with both i.i.d. Gaussian and partial-DCT measurement matrices, which are closer to the matrices used in radar and magnetic resonance imaging applications. We summarize the results of our simulations for Gaussian matrices and partial-DCT matrices in Sections V-C and V-D, respectively.

V-C Compressive imaging with i.i.d. measurement matrices

V-C1 Noiseless

In Table I, we compare the results of JPEG-GD and JP2K-GD with that of the state-of-the-art BM3D-AMP method in reconstructing several test images from their compressive measurements. In these simulations, the measurement matrices are i.i.d. Gaussian. Furthermore, the test images are resized to $128\times 128$ pixels. We consider two sampling ratios $\frac{m}{n}=$ $30\%$ and $50\%$. At each iteration, the step-size parameter $\eta$ is set by solving (18). To find the solution of this optimization problem we used $K_{2,\max}=25$ iterations of Nelder-Mead method algorithm. Furthermore, we used the stopping criteria discussed in Section V-A. For BM3D-AMP we also used the default setting proposed in [9].

The results of this simulation are presented in Table I. Interestingly, results in Table III indicate that JP2K-GD considerably outperforms JPEG-GD. The main reason is that JP2K codec exploits more complex structures of natural images compared with JPEG codec. Therefore, intuitively we can state that for a same amount of measurements, JP2K-GD can perform better. Also the performance of JP2K-GD is comparable with that of BM3D-AMP. When an image has more geometry (as in House), BM3D-AMP outperforms JP2K-GD. However, when an image has more irregular structures and texture such as the Dog image or the Panda image, then JP2K-GD seems to often outperform BM3D-AMP.

V-C2 Noisy

In Table II, we present the performance results of our proposed JPEG-GD and JP2K-GD, and compare them with the performance of the BM3D-AMP method for image reconstruction from noisy compressive measurements. Similar to the previous section, the measurement matrix is i.i.d.  Gaussian, the images are resized to $128\times 128$, and two sampling ratios of $30\%$ and $50\%$ are considered. Unlike before, the measurements are corrupted by i.i.d.  Gaussian noise. We consider two different values of ${\rm SNR}=10\,dB$ and ${\rm SNR}=30\,dB$. The results of our simulations are presented in Table II. As is again clear from this table our results are comparable and in most cases better than the results of the state-of-the-art BM3D-AMP.

V-D Compressive imaging with partial-DCT matrices

In many application areas, measurement matrices such as partial-DCT matrices are employed. In this section we evaluate the performance of our algorithm on partial-DCT matrices. We note that in our numerical results we observe that even though BM3D-AMP performs well for i.i.d. Gaussian measurements, its performance degrades dramatically for partial-DCT matrices. Hence, we compare the performance of our algorithm with the state-of-the-art algorithm for partial-DCT (or partial-Fourier) matrices, i.e., NLR-CS [11]. As in the previous section, we consider both noiseless and noisy measurements.

V-D1 Noiseless

In Table III, we present the performance results of our proposed JPEG-GD and JP2K-GD, and compare them with the performance of the NLR-CS method in image reconstruction from compressive samples, sampled by a random partial-DCT measurement matrix. Images in this numerical comparison are resized to $512\times 512$. We consider two sampling ratios $10\%$ and $30\%$. For $\frac{m}{n}=10\%$ (sampling rate), JP2K-GD performs comparable and in some cases, e.g., Dog and Snake, even better than the state-of-the-art NLR-CS. Increasing the sampling ratio to $\frac{m}{n}=30\%$, JP2K-GD outperform both JPEG-GD and NLR-CS methods.

Note that NLR-CS method has two main steps. In the first step, it estimates an initial image $\hat{{\boldsymbol{x}}}$ using a standard compressed sensing (CS) recovery method based on the sparsity of image coefficients in DCT/Wavelet domain. Then, in the second step, it enforces a low-rankness and group-sparsity constraint on the group of similar patches detected in the estimated image [11]. This step involves singular value decomposition of a matrix and hence is computationally expensive for large images. Furthermore, since in the second step, detection of similar patches is performed using an estimated image, exploiting structures in the second step heavily relies on the performance of the first step. For this particular reason, as observed in Section V-D2, the performance of NLR-CS method degrades significantly once noise is added to the observations. On the other hand, results in Section V-D2 show that both JPEG-GD and JP2K-GD are robust to the measurement noise.

V-D2 Noisy

In Table IV, we present the performance results of JPEG-GD and JP2K-GD, and NLR-CS methods for image reconstruction from noisy compressive measurements. Similar to the previous section, the measurement matrix is a random partial-DCT matrix. Images in this numerical comparison are resized to $512\times 512$, and we consider two sampling ratios $10\%$ and $30\%$. The measurements are corrupted by i.i.d.  Gaussian noise. We consider two different values of ${\rm SNR}=10\,dB$ and ${\rm SNR}=30\,dB$. The results of our simulations are presented in Table IV. As is again clear from this table JP2K-GD method outperforms both JPEG-GD and NLR-CS for all SNRs and sampling ratios ($m/n$). Interestingly we observe that for low SNRs, e.g., SNR = $10$ dB, even the JPEG-GD algorithm performs much better than NLR-CS.

Our next goal is to visually compare the reconstruction of JP2K-GD with that of NLR-CS. Fig. 2 shows the reconstructed images for three different sampling ratios $m/n$. The size of the test image for all scenario is $512\times 512$ and the measurement SNR is set to $30$ dB. As is clear from Fig. 2, in all cases, the reconstruction from JP2K-GD looks more appealing than the reconstruction from NLR-CS.

V-E Convergence rate evaluation

As proved in our main theorems, we expect the convergence to be linear when the measurement matrix is Gaussian. It turns out that the convergence is also fast for partial-DCT matrices. Fig. 3 depicts the normalized mean-squared-error (MSE) of image reconstruction using JP2K-GD method. We consider two different sampling ratios $\frac{m}{n}=5\%$ and $10\%$ in this test. Results in Fig. 3 and 4 show that (i) the algorithm converges very fast (often in less than 50 iterations), (ii) by increasing the number of measurements the convergence of JP2K-GD improves, and (iii) the final reconstructed image has a better PSNR. Note that all these conclusions are consistent with the results we proved for sub-Gaussian matrices.

VI Proofs

VI-A Background

In this section we briefly review some useful results that are going to be used in the proofs.

Lemma 2** (see Lemma 5.9 in [21]).**

Let $\left\{X_{i}\right\}_{i=1}^{n}$ be independent, mean zero, sub-Gaussian random variables and $\left\{a_{i}\right\}_{i=1}^{n}$ are real numbers. Then $\sum_{i=1}^{n}a_{i}X_{i}$ is also a sub-Gaussian random variable, and

[TABLE]

Theorem 5** (Bernstein Type Inequality, see e.g., [21]).**

Suppose that $\left\{X_{i}\right\}_{i=1}^{n}$ are independent, and that, for $i=1,\cdots,n$, $X_{i}$ is a sub-exponential random variable. Let $\max_{i}{\left\|{X_{i}}\right\|}_{\psi_{1}}\leq K$, for some $K>0$. Then for every $t\geq 0$ and every ${\boldsymbol{w}}=\left[w_{1},\cdots,w_{n}\right]^{T}\in\mathbbm{R}^{n\times 1}$, we have

[TABLE]

Lemma 3** (Lemma 3 from [22]).**

Consider two independent random vectors ${\boldsymbol{X}}=[X_{1},\cdots,X_{n}]^{T}\in\mathbbm{R}^{n}$ and ${\boldsymbol{Y}}=[Y_{1},\cdots,Y_{n}]^{T}\in\mathbbm{R}^{n}$. Assume that $X_{i}\stackrel{{\scriptstyle\rm i.i.d.}}{{\sim}}\mathcal{N}(0,1)$ and $Y_{i}\stackrel{{\scriptstyle\rm i.i.d.}}{{\sim}}\mathcal{N}(0,1)$. Then $\left<{\boldsymbol{X}},{\boldsymbol{Y}}\right>$ and $G{\left\|{{\boldsymbol{X}}}\right\|}_{2}$ have the same distribution, where $G\sim\mathcal{N}(0,1)$ and is independent of ${\left\|{{\boldsymbol{X}}}\right\|}_{2}$.

Lemma 4** (Lemma 2 from [22]).**

Let $G_{i}$, $i=1,2,\cdots,m$, be i.i.d. $\mathcal{N}(0,1)$. Then, for $\tau\in(0,1)$,

[TABLE]

and for $\tau>0$,

[TABLE]

Theorem 6** (see e.g., [23]).**

Let ${\boldsymbol{\rm A}}\in\mathbbm{R}^{m\times n}$ be a dense random matrix whose entries are i.i.d. zero-mean Gaussian random variables with unit variance. Then, $t>0$,

[TABLE]

If we substitute $t\leftarrow t\sqrt{m}$, then Theorem 6 results in Corollary 4.

Corollary 4**.**

Let ${\boldsymbol{\rm A}}\in\mathbbm{R}^{m\times n}$ be a random matrix whose entries are independent, zero-mean Gaussian random variables with unit variance. Then

[TABLE]

$t>0$* is arbitrary variable.*

Theorem 7**.**

Let ${\boldsymbol{\rm A}}\in\mathbbm{R}^{m\times n}$ be an i.i.d. matrix such that $A_{i,j}$ is a zero-mean sub-Gaussian random variable with ${\left\|{{A}_{i,j}}\right\|}_{\psi_{2}}\leq K$ and ${\rm E}\left[{A}_{i,j}^{2}\right]=\sigma_{a}^{2}$. Then, if $m<n$, for any $t>0$,

[TABLE]

Proof.

Let $\mathcal{N}_{\varepsilon}$ denote a maximal $\varepsilon$-separated subset of $S^{n-1}$. It is straightforward to show that [Lemma 5.2 from [23]]

[TABLE]

Consider vector ${\boldsymbol{u}}\in S^{n-1}$ that satisfies $\sigma_{\max}(A)=\max_{{\boldsymbol{u}}}\|A{\boldsymbol{u}}\|$. Let ${\boldsymbol{u}}_{0}\in\mathcal{N}_{n,\varepsilon}$ be such that ${\left\|{{\boldsymbol{u}}-{\boldsymbol{u}}_{0}}\right\|}_{2}\leq\varepsilon$. Then, by the triangle inequality, we have

[TABLE]

On the other hand, again by the triangle inequality,

[TABLE]

Combining (27) and (28) yields

[TABLE]

To finish the proof, we need to upper bound $\max_{{\boldsymbol{x}}\in\mathcal{N}_{n,\varepsilon}}\left<{\boldsymbol{\rm A}}{\boldsymbol{x}},{\boldsymbol{\rm A}}{\boldsymbol{x}}\right>$. For a fixed ${\boldsymbol{x}}\in\mathcal{N}_{\varepsilon}$, by Theorem 5,

[TABLE]

Therefore,

[TABLE]

Let $\varepsilon={1\over 4}$. Then, from (26),

[TABLE]

Substituting $t$ by $6n(1+t)K/m$ and noting that for this value of $t$, if $m<3n$, ${t\over 2K}$ is always larger than 1, it follows that

[TABLE]

Therefore, in summary, form (29), with probability larger than $1-{\rm e}^{-3nt}$,

[TABLE]

∎

Substituting $t$ by $\frac{m\sigma_{a}^{2}}{6nK}$ in Theorem 7 results in the following corollary.

Corollary 5**.**

Let ${\boldsymbol{\rm A}}\in\mathbbm{R}^{m\times n}$ be an i.i.d. matrix such $A_{i,j}$ is a zero-mean sub-Gaussian random variable with ${\left\|{{A}_{i,j}}\right\|}_{\psi_{2}}\leq K$ and ${\rm E}\left[{A}_{i,j}^{2}\right]=\sigma_{a}^{2}$. Then, for $m<n$,

[TABLE]

Lemma 5** (Lemma 5 in [14].).**

Consider ${\boldsymbol{u}},{\boldsymbol{v}}\in S^{n-1}$ and dense random Gaussian matrix ${\boldsymbol{\rm A}}\in\mathbbm{R}^{m\times n}$ with i.i.d zero mean Gaussian entries as $\mathcal{N}\left(0,\sigma_{a}^{2}\right)$. Then, for any $t>0$

[TABLE]

where $f^{*}(t)=\min_{u\in[-1,1]}\max_{s\in\left(0,\frac{1}{1-u}\right)}\left\{s(t-u)+\frac{1}{2}\ln\left[\left(1+su\right)^{2}-s^{2}\right]\right\}.$

Corollary 6**.**

Consider ${\boldsymbol{u}},{\boldsymbol{v}}\in S^{n-1}$ and dense random Gaussian matrix ${\boldsymbol{\rm A}}\in\mathbbm{R}^{m\times n}$ with i.i.d. zero mean Gaussian entries as $\mathcal{N}\left(0,\sigma_{a}^{2}\right)$. Then,

[TABLE]

Lemma 6**.**

Consider ${\boldsymbol{u}},{\boldsymbol{v}}\in S^{n-1}$ and dense matrix ${\boldsymbol{\rm A}}\in\mathbbm{R}^{m\times n}$ with i.i.d. zero-mean sub-Gaussian entries with ${\left\|{{\boldsymbol{\rm A}}_{i,j}}\right\|}_{\psi_{2}}\leq K$ and ${\rm E}\left[A_{i,j}^{2}\right]=\sigma_{a}^{2}$. Then, for any $t>0$,

[TABLE]

Proof.

Define ${\boldsymbol{y}}(u)={\boldsymbol{\rm A}}{\boldsymbol{u}}$ and ${\boldsymbol{y}}(v)={\boldsymbol{\rm A}}{\boldsymbol{v}}$. Using these definitions, $\left<{\boldsymbol{\rm A}}{\boldsymbol{u}},{\boldsymbol{\rm A}}{\boldsymbol{v}}\right>=\left<{\boldsymbol{y}}(u),{\boldsymbol{y}}(v)\right>$. Let ${\boldsymbol{\rm A}}_{i}\in\mathbbm{R}^{1\times n}$ denote the $i$-the row of matrix ${\boldsymbol{\rm A}}$. Thus, ${\boldsymbol{y}}_{i}(u)=\left<{\boldsymbol{\rm A}}_{i},{\boldsymbol{u}}\right>$ and ${\boldsymbol{y}}_{i}(v)=\left<{\boldsymbol{\rm A}}_{i},{\boldsymbol{v}}\right>$ are both sub-Gaussian random variables. Using Lemma 2, we have

[TABLE]

Furthermore, ${\rm E}\left[{\boldsymbol{y}}_{i}(u){\boldsymbol{y}}_{i}(v)\right]={\boldsymbol{u}}^{T}{\rm E}\left[{\boldsymbol{\rm A}}_{i}^{T}{\boldsymbol{\rm A}}_{i}\right]{\boldsymbol{v}}=\sigma_{a}^{2}\left<{\boldsymbol{u}},{\boldsymbol{v}}\right>.$ Note that

[TABLE]

By Lemma 1, ${\boldsymbol{y}}_{i}(u){\boldsymbol{y}}_{i}(v)$ is a sub-exponential random variable with

[TABLE]

Therefore, by applying Theorem 5 to sub-exponential random variables ${\boldsymbol{y}}_{i}(u){\boldsymbol{y}}_{i}(v)$, and setting all weights equal to $-1$, we derive

[TABLE]

∎

VI-B Proof of Theorem 2

Define

[TABLE]

Using this notation, we have ${\boldsymbol{x}}^{k+1}=\mathcal{P}_{\mathcal{C}}({\boldsymbol{s}}^{k+1})$. But since $\tilde{{\boldsymbol{x}}}=\mathcal{P}_{\mathcal{C}}({\boldsymbol{x}})$, $\tilde{{\boldsymbol{x}}}$ is also in $\mathcal{C}$. Hence, ${\left\|{{\boldsymbol{x}}^{k+1}-{\boldsymbol{s}}^{k+1}}\right\|}_{2}^{2}\leq{\left\|{\tilde{{\boldsymbol{x}}}-{\boldsymbol{s}}^{k+1}}\right\|}_{2}^{2},$ or, equivalently, ${\left\|{\left({\boldsymbol{x}}^{k+1}-\tilde{{\boldsymbol{x}}}\right)-\left({\boldsymbol{s}}^{k+1}-\tilde{{\boldsymbol{x}}}\right)}\right\|}_{2}^{2}\leq{\left\|{\tilde{{\boldsymbol{x}}}-{\boldsymbol{s}}^{k+1}}\right\|}_{2}^{2}$. By removing the common terms from both sides, we have

[TABLE]

For $k=0,1,\ldots$, define the error vector and its normalized version as

[TABLE]

and

[TABLE]

respectively. Also, given ${\boldsymbol{\theta}}^{k}\in\mathbbm{R}^{n}$, ${\boldsymbol{\theta}}^{k+1}\in\mathbbm{R}^{n}$, $\eta\in\mathbbm{R}^{+}$, and ${\boldsymbol{\rm A}}\in\mathbbm{R}^{m\times n}$, define coefficient $\mu$ as

[TABLE]

Using this definition, substituting for ${\boldsymbol{s}}^{k+1}$ from (37) and noting that ${\boldsymbol{y}}={\boldsymbol{\rm A}}{\boldsymbol{x}}+{\boldsymbol{z}}$, from (38), it follows that

[TABLE]

where $\|{\boldsymbol{\rm A}}\|_{S^{n-1}}=\sigma_{\max}\left({\boldsymbol{\rm A}}\right)$. We next find upper bounds for the three terms on the right hand side of (39):

1. (i)

Bounding $\mu\left(\underline{{\boldsymbol{\theta}}}^{k+1},\underline{{\boldsymbol{\theta}}}^{k},\eta\right)$: We show that given the parameter setting of the theorem, with high probability

[TABLE]

where

[TABLE]

To achieve this goal, we define event $\mathcal{E}_{1}$ as

[TABLE]

From Corollary 6 (or Lemma 5), given ${\boldsymbol{u}},{\boldsymbol{v}}\in\mathcal{C}^{\prime}$, we have

[TABLE]

Therefore, by the union bound,

[TABLE]

Note that $|\mathcal{C}^{\prime}|\leq|\mathcal{C}|^{2}\leq 2^{2r}$. Therefore,

[TABLE]

Therefore, for $m\geq 80r\left(1+\epsilon\right)$, where $\epsilon>0$, with probability at least $1-2^{-40r\epsilon}$, event $\mathcal{E}_{1}$ happens. 2. (ii)

Bounding $\|{\boldsymbol{\rm A}}\|_{S^{n-1}}^{2}{\left\|{{\boldsymbol{x}}-\tilde{{\boldsymbol{x}}}}\right\|}_{2}$: Define event $\mathcal{E}_{2}$ as

[TABLE]

From Corollary 4, for $t=1$ we have

[TABLE]

Also, since the compression code has supremum distortion $\delta$, ${\left\|{{\boldsymbol{x}}-\tilde{{\boldsymbol{x}}}}\right\|}_{2}\leq\delta$. Therefore, conditioned on $\mathcal{E}_{2}$, we have

[TABLE] 3. (iii)

Bounding $2\eta\left<\underline{{\boldsymbol{\theta}}}^{k+1},{\boldsymbol{\rm A}}^{T}{\boldsymbol{z}}\right>$: Note that $2\eta\left<\underline{{\boldsymbol{\theta}}}^{k+1},{\boldsymbol{\rm A}}^{T}{\boldsymbol{z}}\right>=\frac{2}{m\sigma_{a}^{2}}\left<{\boldsymbol{\rm A}}\underline{{\boldsymbol{\theta}}}^{k+1},{\boldsymbol{z}}\right>$. Let ${\boldsymbol{\rm A}}_{i}\in\mathbbm{R}^{n}$ be the $i$-th row of matrix ${\boldsymbol{\rm A}}$. Then,

[TABLE]

For any fixed $\underline{{\boldsymbol{\theta}}}^{k+1}$, $\left\{\left<{\boldsymbol{\rm A}}_{i},\underline{{\boldsymbol{\theta}}}^{k+1}\right>\right\}_{i=1}^{n}$ are i.i.d. $\mathcal{N}(0,1)$ random variables. Hence, from Lemma 3, we know that the distribution of $\left<\underline{{\boldsymbol{\theta}}}^{k+1},{\boldsymbol{\rm A}}^{T}{\boldsymbol{z}}\right>$ is the same as $\sigma_{a}{\left\|{{\boldsymbol{z}}}\right\|}_{2}\left<\underline{{\boldsymbol{\theta}}}^{k+1},{\boldsymbol{g}}\right>$, where ${\boldsymbol{g}}=\left[g_{1},\cdots,g_{n}\right]^{T}$ is independent of ${\left\|{{\boldsymbol{z}}}\right\|}_{2}$ and $g_{i}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathcal{N}(0,1)$. To bound $\sigma_{a}{\left\|{{\boldsymbol{z}}}\right\|}_{2}\left<\underline{{\boldsymbol{\theta}}}^{k+1},{\boldsymbol{g}}\right>$ we will bound $\frac{1}{\sigma_{z}^{2}}{\left\|{{\boldsymbol{z}}}\right\|}_{2}^{2}$ and $\left|\left<\underline{{\boldsymbol{\theta}}},{\boldsymbol{g}}\right>\right|^{2}$ separately. Given $\tau^{\prime}_{1}>0$ and $\tau^{\prime}_{2}>0$, define events $\mathcal{E}_{3}$ and $\mathcal{E}_{4}$ as follows

[TABLE]

and

[TABLE]

Following Lemma 4, we have

[TABLE]

and letting $m=1$ in Lemma 4, for fixed $\underline{{\boldsymbol{\theta}}}^{k+1}$, it follows that

[TABLE]

Hence, by the union bound,

[TABLE]

where the last inequality holds for $\tau^{\prime}_{2}>7$. Setting $\tau^{\prime}_{2}=4(1+\epsilon)r-1$, where $\epsilon>0$, ensures that $P\left(\mathcal{E}_{4}^{c}\right)\leq 2^{-2\epsilon r+0.5}$. Setting $\tau^{\prime}_{1}=1$, $\mathcal{P}(\mathcal{E}_{3}^{c})\leq{\rm e}^{-0.15m}$, and conditioned on $\mathcal{E}_{3}\cap\mathcal{E}_{4}$, we have

[TABLE]

Combining (49), (45), and (43) with (39) yields the desired bound on the reduction of error. Finally, note that, by the union bound,

[TABLE]

VI-C Proof of Theorem 3

The proof of this result is a simple extension of the proof of Theorem 2 presented in Section VI-B. Define

[TABLE]

Note that ${\boldsymbol{x}}^{k+1}=g_{r}(f_{r}({\boldsymbol{s}}^{k+1}))$. Hence,

[TABLE]

It is now straightforward to follow exactly the same step as the one discussed in the proof of Theorem 2 and show that with probability at least $1-2^{-4r\epsilon}-{\rm e}^{-{m\over 4}}-2^{-2r\epsilon}$

[TABLE]

Combining (51) and (52) completes the proof.

VI-D Proof of Theorem 4

Following the proof of Theorem 2, and defining ${\boldsymbol{\theta}}^{k+1}={\boldsymbol{x}}^{k+1}-\tilde{{\boldsymbol{x}}}$, for $k=0,1,2,\ldots$, it follows from (39) that

[TABLE]

where $\mu\left(\underline{{\boldsymbol{\theta}}}^{k+1},\underline{{\boldsymbol{\theta}}}^{k},\eta\right)\triangleq\left<\underline{{\boldsymbol{\theta}}}^{k+1},\underline{{\boldsymbol{\theta}}}^{k}\right>-\eta\left<{\boldsymbol{\rm A}}\underline{{\boldsymbol{\theta}}}^{k+1},{\boldsymbol{\rm A}}\underline{{\boldsymbol{\theta}}}^{k}\right>$. Let $\mathcal{C}^{\prime}$ denote the set of normalized distance vectors of the codewords in $\cal C$, which is defined in (41). Define event $\mathcal{E}_{1}$ as

[TABLE]

Similar to the proof of Theorem 2, we show that the probability of occurrence of $\mathcal{E}_{1}^{c}$ approaches [math]. Given ${\boldsymbol{u}},{\boldsymbol{v}}\in\mathcal{C}^{\prime}$, from Lemma 6, we have

[TABLE]

Therefore, by the union bound, since $|\mathcal{C}^{\prime}|\leq|\mathcal{C}|^{2}=2^{2r}$, we have

[TABLE]

Therefore, for $m>\frac{8K^{2}r(1+\epsilon)}{\mu_{0}\sigma_{a}^{2}\min\left(1,\frac{\mu_{0}\sigma_{a}^{2}}{2K^{2}}\right)\log{\rm e}}$,

[TABLE]

But, since by assumption $\mu_{0}\sigma_{a}^{2}\leq 2K^{2}$, we have $\min\left(1,\frac{\mu_{0}\sigma_{a}^{2}}{2K^{2}}\right)=\frac{\mu_{0}\sigma_{a}^{2}}{2K^{2}}$.

Define event $\mathcal{E}_{2}$ as $\mathcal{E}_{2}\triangleq\left\{\sigma_{\max}({\boldsymbol{\rm A}})\leq 2\sigma_{a}\sqrt{m+n\frac{3K}{\sigma_{a}^{2}}}\right\}$. From Corollary 5, we have $P\left(\mathcal{E}_{2}\right)\geq 1-{\rm e}^{-\frac{m\sigma_{a}^{2}}{2K}}$. Since the compression code is such that ${\left\|{{\boldsymbol{x}}-\tilde{{\boldsymbol{x}}}}\right\|}_{2}\leq\delta$, conditioned on $\mathcal{E}_{2}$, we have

[TABLE]

To complete the proof, we need to bound $2\eta\left<\underline{{\boldsymbol{\theta}}}^{k+1},{\boldsymbol{\rm A}}^{T}{\boldsymbol{z}}\right>=\frac{2}{m\sigma_{a}^{2}}\left<{\boldsymbol{\rm A}}\underline{{\boldsymbol{\theta}}}^{k+1},{\boldsymbol{z}}\right>$, which is the term related to the noise ${\boldsymbol{z}}$. Again, let ${\boldsymbol{\rm A}}_{i}\in\mathbbm{R}^{1\times n}$ denote the $i$-th row of matrix ${\boldsymbol{\rm A}}$ and, for a given ${\boldsymbol{u}}\in\mathbbm{R}^{n}$, let ${\boldsymbol{y}}_{{\boldsymbol{u}}}={\boldsymbol{\rm A}}{{\boldsymbol{u}}}.$ Hence, for $i=1,\ldots,m$,

[TABLE]

To upper bound the term corresponding to noise, given $\tau>0$, define event $\mathcal{E}_{3}$ as

[TABLE]

From Lemma 2, for ${\boldsymbol{u}}\in\mathcal{C}^{\prime}$, we know that $\left\{{\boldsymbol{y}}_{{\boldsymbol{u}}}(i)\right\}_{i=1}^{m}$ are independent sub-Gaussian random variables. Also, for ${\boldsymbol{u}}\in\mathcal{C}^{\prime}$, Lemma 2 states that

[TABLE]

where the last inequality follows because for ${\boldsymbol{u}}\in\mathcal{C}^{\prime}$, ${\left\|{{\boldsymbol{u}}}\right\|}_{2}=1$. Since every Gaussian random variable is also a sub-Gaussian random variable, ${\boldsymbol{z}}(i)$ is a sub-Gaussian random variable with ${\left\|{{\boldsymbol{z}}(i)}\right\|}_{\psi_{2}}=\sigma_{n}\sqrt{8\over 3}$. As a result, ${\left\|{{\boldsymbol{y}}_{{\boldsymbol{u}}}(i){\boldsymbol{z}}(i)}\right\|}_{\psi_{1}}\leq{\left\|{{\boldsymbol{y}}_{{\boldsymbol{u}}}(i)}\right\|}_{\psi_{2}}{\left\|{{\boldsymbol{z}}(i)}\right\|}_{\psi_{2}}\leq K\sqrt{8\over 3}\sigma_{n}$. Using Theorem 5, for ${\boldsymbol{u}}\in\mathcal{C}^{\prime}$, we have

[TABLE]

where the last line follows because ${3\over 8}>{1\over 3}$. Therefore, by the union bound, since $|\mathcal{C}^{\prime}|\leq 2^{2r}$,

[TABLE]

Choosing

[TABLE]

and given our choice of $m$, it follows that

[TABLE]

But, since $\sqrt{96\over\log{\rm e}}\leq 9$,

[TABLE]

Finally, combining (53) with the bounds derived on the three terms on the right hand side of (53) yields the desired result.

VII Conclusions

In this paper, we have studied the problem of designing efficient compression-based compressed sensing recovery algorithms. Specifically, we have proposed C-GD, an iterative robust-to-noise compression-based compressed sensing algorithm. Given measurements ${\boldsymbol{y}}={\boldsymbol{\rm A}}{\boldsymbol{x}}+{\boldsymbol{z}}$ and a compression code with codebook $\mathcal{C}$, at iteration $k$, C-GD updates its current estimate of ${\boldsymbol{x}}$, ${\boldsymbol{x}}^{k}$, by moving towards the negative of the gradient of the cost function ($f({\boldsymbol{u}})=\|{\boldsymbol{y}}-{\boldsymbol{\rm A}}{\boldsymbol{u}}\|^{2}$) and then projecting the result onto the set of codewords $\mathcal{C}$. For a given compression code, the projection step can typically be implemented by applying the compression code’s encoder and decoder. We have proved that, given enough measurements, with high probability, C-GD has linear convergence rate and is robust to additive Gaussian noise. In summary, C-GD provides a platform for using commercial compression code such as JPEG2000 or MPEG4 for compressed sensing of images and videos, respectively. In our simulation results, we have focused on compressed sensing of images and have shown that C-GD combined with the state-of-the-art compression codes yields state-of-the-art compressed sensing performance, both for i.i.d. Gaussian and partial Fourier measurement matrices.

Appendix A Finding the best piecewise polynomial approximation

Consider the following problem: given ${\boldsymbol{x}}\in\mathds{R}^{n}$, find $\tilde{{\boldsymbol{x}}}\in\mathcal{P}$, $\mathcal{P}$ defined in (14), such that

[TABLE]

In this section, we briefly describe how $\tilde{{\boldsymbol{x}}}$ can be found using dynamic programing. Note that, given singularity points $s_{1},s_{2},\ldots,s_{Q}$, one can easily find the best polynomial fit in each piece. Hence, the challenge is to find the optimal singularity points. Each singularity point $s_{i}$ is a point in the set $\{{1\over n},\ldots,{n-1\over n}\}$. Given $i_{1},i_{2}\in\{0,\ldots,n\}$, $i_{1}\leq i_{2}$, let $e(i_{1},i_{2})$ denote the minimum error achievable in approximating $(x_{i_{1}},\ldots,x_{i_{2}})$ as samples of $\sum_{j=0}^{N}a_{j}y^{j}$ at ${i_{1}\over n},\ldots,{i_{2}\over n}$, where $\sum_{j=0}^{n}a_{j}\leq 1$ and $a_{j}\in(0,1)$, $j=0,\ldots,N$. That is,

[TABLE]

Using this definition, given singularity points $s_{0}=0,s_{1},\ldots,s_{Q},s_{Q+1}=1\in\{0,{1\over n},\ldots,1\}$, the minimum achievable error in approximating ${\boldsymbol{x}}$ by signals in $\mathcal{P}$ whose singularities happen at $s_{1},\ldots,s_{Q}$ can be written as

[TABLE]

This representation suggests that the minimizer $\tilde{{\boldsymbol{x}}}$ can be found using the Viterbi algorithm. In summary, the Viterbi algorithm will operate on a Trellis diagram with $Q$ full stages corresponding to the possible $Q$ singularities and two single-state stages corresponding to the start and the end of the interval. Each intermediate stage has $n-1$ states, which correspond to the possible $n-1$ singularity points. State $s$ at stage $t\in{1,\ldots,Q}$ is connected to state $s^{\prime}$ at state $t+1$, if $s\leq s^{\prime}$. The weight of this edge is set as $e(ns,ns^{\prime})$, defined in (A.1). Otherwise, if $s^{\prime}>s$, there is no edge between the two states. Let $E_{i}(s)$ denote the minimum cost associated with state $s$ at stage $i$. Also, let $E_{0}(s_{0})=0$. The goal is to find the path on the Trellis diagram that achieves $E_{Q+1}(s_{Q+1})=E_{Q+1}(1)$. It is straightforward to show that, for $t=1,\ldots,Q$,

[TABLE]

where the minimum is taken over all states $s^{\prime}$ which are connected to $s$, i.e., $s^{\prime}<s$. This breakdown of the cost function describes the essence of how the Viterbi algorithm operates. At stage $t$, among its incoming edges, each state $s$ only keeps the edge that achieve $E_{t+1}(s)$. At the end, backtracking from the final state $s_{Q+1}=1$ at stage $Q+1$ reveals the optimal singularities.