# Compressed sensing in the presence of speckle noise

Wenda Zhou  Shirin Jalali  Arian Maleki  Wenda Zhou, Shirin Jalali, Arian
Maleki W. Zhou is with the NYU’s Center for Data Science, New York, NY
(e-mail: wz2247@nyu.edu).S. Jalali is with Nokia Bell Labs, Murray Hills, NJ
(e-mail: shirin.jalali@nokia-bell-labs.com).A. Maleki is with the Statistics
Department of Columbia University, New York, NY. (e-mail:
arian@stat.columbia.edu).

###### Abstract

The problem of recovering a structured signal from its linear measurements in
the presence of speckle noise is studied. This problem appears in many imaging
systems such as synthetic aperture radar and optical coherence tomography. The
current acquisition technology oversamples signals and converts the problem
into a denoising problem with multiplicative noise. However, this paper
explores the possibility of reducing the number of measurements below the
ambient dimension of the signal. The sophistications that appear in the study
of multiplicative noises have so far impeded theoretical analysis of such
problems. This paper aims to present the first theoretical result regarding
the recovery of signals from their undersampled measurements under the speckle
noise. It is shown that if the signal class is structured, in the sense that
the signals can be compressed efficiently, then one can obtain accurate
estimates of the signal from fewer measurements than the ambient dimension. We
demonstrate the effectiveness of the methods we propose through simulation
results.

##  I Introduction

###  I-A Problem statement

Various modern imaging methods, such as synthetic aperture radar (SAR) [1] and
optical coherence tomography (OCT) [2], rely on coherent imaging. The main
drawback of such coherence-based imaging systems is that they all suffer from
a granular noise that is typically referred to as speckle (or multiplicative)
noise. At a high level, the reason for observing such a noise is the
following. In applications of such imaging methods, surfaces reflecting the
incoming coherent waves are all rough, when considered at a resolution
comparable to the wavelength of the signal. Therefore, the reflected rays will
have different phases, which means that at detection point, reflected rays
corresponding to some pixels will add up constructively, while reflected
points corresponding to other pixels add up destructively. This phenomenon can
be modeled as a multiplicative noise.

In this paper, our goal is to develop a theoretical framework for studying the
imaging problem in the presence of speckle noise. We focus on the problem of
compressed sensing, i.e., recovering a structured signal from its
underdetermined measurements in the presence of such noise. Different speckle
noise filtering methods (for the case that the measurement matrix is identity)
have been proposed in the literature over the past couple of decades. However,
to the best of our knowledge, there has not been any theoretical analysis of
speckle noise and corresponding optimal recovery methods when fewer
measurements than the ambient dimension are available. As will be discussed
later in the paper, compared to the well-studied traditional additive noise,
which appears in other imaging systems such as magnetic resonance imaging,
multiplicative noise poses many more theoretical and practical challenges. In
this paper, we discuss and address some of these challenges.

Let ${\cal Q}$ denote a compact subset of $\mathds{R}^{n}$ that describes the
class of structured signals we are interested in, e.g., class of natural
images or more abstractly, class of bounded $k$-sparse signals or class of
bounded piece-wise-constant signals. For ${\bf x}\in\mathds{R}^{n}$, define
$X=\mathop{\rm diag}\nolimits({\bf x})=\mathop{\rm
diag}\nolimits(x_{1},\ldots,x_{n})$. In imaging systems with speckle noise,
the measurement vector ${\bf y}\in\mathds{R}^{m}$ is defined as

${\bf y}=AX{\bf w}+{\bf z}.$

Here, the multiplicative noise ${\bf w}\in\mathds{R}^{n}$ and additive noise
${\bf z}\in\mathds{R}^{m}$ are i.i.d. $\mathcal{N}(0,\sigma_{w}^{2})$ and
i.i.d. $\mathcal{N}(0,\sigma_{z}^{2})$, respectively.111Note that in real-
world systems, the speckle noise is typically complex-valued and has a
Gaussian distribution. The signal and the measurement matrix are also complex-
valued. However, for notational simplicity we have focused on real-valued
signals and noise here. Extension of the results to complex-valued signals and
noises is straightforward.  The difference between such a system and a
standard linear measurement system is clearly in the multiplicative noise
${\bf w}$ that distorts each input pixel independently. Note that in such
systems, unlike the additive noise ${\bf z}$, the dimensions of the
multiplicative noise depends on the dimensions of the input signal ${\bf x}$
and not the number of measurements $m$. Given such a measurement process,
there are various fundamental questions one can ask:

  1. 1.

Is accurate recovery feasible in the presence of speckle noise?

  2. 2.

What is the performance of maximum likelihood-based (ML-based) recovery method
that takes the source structure into account?

  3. 3.

Given a class of structured signals $\mathcal{Q}$ with a certain intrinsic
dimension, what is the minimum number of measurements $m$ that guarantees the
feasibility of accurate recovery?

In this paper, we aim to address these questions. First, inspired by maximum
likelihood estimation, we derive a compression-based ML recovery method that
employs compression codes designed for our desired class of signals
($\mathcal{Q}$) to define, capture and exploit the source structure. We
characterize the performance of the derived compression-based ML recovery
method and prove that given enough number of measurements (related to the
desired accuracy and rate-distortion performance of the compression code), it
is able to recover the signal within desired accuracy, with probability
approaching one, as the dimensions of the problem grow without bound.

###  I-B Related work

As mentioned earlier, speckle noise is an inherent problem in coherence-based
imaging systems. The current technology in such imaging systems is to collect
at least as many measurements as the ambient dimension of the signal, and
convert the problem into a denoising problem, i.e. recovering the signal ${\bf
x}$ from

${\bf y}={\bf x}\odot{\bf w}+{\bf z},$

where $\odot$ denotes element-wise multiplication, i.e., $y_{i}=x_{i}w_{i}$,
and ${\bf w}$ and ${\bf z}$ denote iid Gaussian noises.

For each type of such imaging systems, various denoising methods have been
developed over the years. One such type of imaging that is most relevant to
the measurement model studied in this paper is SAR imaging [1]. The classical
techniques for denoising speckle noise in SAR imaging have been reviewed in
[3] and [4]. In [5], the authors propose a sparsity-based total-variation
approach for smoothing and denoising speckled images. Application of non-local
means approaches to speckle noise reduction are explored in [6, 7]. In recent
years, inspired by the success of deep neural networks, especially
convolutional neural nets, in solving various inference tasks involving
images, researchers have explored application of such tools in SAR speckle
noise denoising as well [8, 9].

One of the well-known applications of compressed sensing is inverse SAR (ISAR)
imaging. Compressed sensing enables such systems to produce images of equal or
even higher quality in shorter acquisition time compared to conventional
systems. This is achieved by requiring much lower number of samples [10, 11].
Various challenges faced in applications of compressed sensing to ISAR imaging
have been explored in the literature. (Refer to [12, 13, 14, 15, 16] for a
noncomprehensive list of such works.) However, compressed sensing in the
presence of speckle noise is not addressed directly in the literature, and
unlike conventional setups, the issue of speckle noise, while still present,
is ignored. In this paper, we address this gap and study the problem of
compressed sensing in the presence of speckle noise. We take a theoretical
approach to the problem and derive a compression-based ML compressed sensing
recovery method in the presence of speckle noise. We characterize the
performance of the derived solution for the case that the additive noise
approaches zero.

###  I-C Notations and definitions

Sets are denoted by calligraphic letters. The size of a set $\mathcal{A}$ is
denoted by $|\mathcal{A}|$. The $\ell_{2}$-norm of an $m\times n$ matrix $A$
is defined as $\|A\|=\max_{{\bf x}\neq 0}\|A{\bf x}\|_{2}/{\bf x}\|_{2}$. The
Hilbert-Schmidt (or Frobenius) norm of $A$ is $\|A\|_{\rm
HS}=(\sum_{i,j}a_{i,j}^{2})^{0.5}$. Throughout the paper, $\log$ refers to
natural logarithm. For a vector ${\bf x}\in\mathds{R}^{n}$, ${\rm diag}({\bf
x})$ denotes the $n\times n$ diagonal matrix with diagonal elements determined
by ${\bf x}$. For $x\in\mathds{R}$, $b$-bit quantized version of $x$ is
defined as $[x]_{b}=2^{-b}\lfloor 2^{b}x\rfloor$.

###  I-D Paper organization

To show the effect of multiplicative noise and how using the source structure
can improve the recovery performance, Section II reviews a simple denoising
problem in the presence of multiplicative noise. Section III reviews
compressible codes and motivates using compression codes to define and enforce
structure of signals. Section IV characterizes the log likelihood function
corresponding to the problem of compressed sensing in the presence of
multiplicative noise and derives some of its properties. V introduces an
optimization recovery based on compression-based ML for compressed sensing in
the presence of speckle noise is derived. The performance of the proposed
optimization is also analyzed. Section VI discusses application of projected
gradient descent for approximating the solution of the optimization
corresponding to ML recovery. The proofs and simulation results are presented
in Section VIII and Section VII, respectively. Section IX concludes the paper.

##  II An illustrative example

To illustrate how inference in the presence of multiplicative noise can be
done, we start with a simple denoising example. Let

${\bf y}={\bf x}\odot{\bf w},$

where $\odot$ denotes element-wise multiplication, i.e., $y_{i}=x_{i}w_{i}$.
Assume that ${\bf w}\stackrel{{\scriptstyle\rm
i.i.d.}}{{\sim}}\mathcal{N}(0,1)$. Note that if we consider $x_{i}=u_{i}$,
then $y_{i}\sim\mathcal{N}(0,u_{i}^{2})$. Therefore, ignoring the constant
terms, the log likelihood function is equal to

$\log\prod_{i=1}^{n}{1\over u_{i}}\exp(-{y_{i}^{2}\over
2u_{i}^{2}})=\sum_{i=1}^{n}\Big{(}-\log|u_{i}|-{y_{i}^{2}\over
2u_{i}^{2}}\Big{)}.$

As expected, since the noise has a symmetric distribution around zero, the
log-likelihood function is symmetric with respect to $u_{i}$. In other words,
the sign of $u_{i}$ cannot be inferred from the given measurements. Hence, in
the rest of this section we assume that all the $u_{i}$s are positive. The log
likelihood can be optimized over the individual signal components $u_{i}$. The
function $-\log u_{i}-{y_{i}^{2}\over 2u_{i}^{2}}$ is a quasi-concave function
of $u_{i}$ and maximizing the likelihood function leads to $\hat{{\bf x}}_{\rm
ML}=|{\bf y}|$. Hence, the mean square error of the maximum likelihood
estimate is given by

$\operatorname{E}[\|\hat{{\bf x}}_{\rm ML}-{\bf x}\|_{2}^{2}]=2\|{\bf
x}\|_{2}^{2}(1-\sqrt{\frac{2}{\pi}}).$ (1)

Now, we would like to consider a class of structured signals and see the
effect of the structure on the MSE. Consider a very simple structure for the
signal. Let $\mathcal{Q}=\\{a{\bf 1}_{n}:\;a\in(x_{m},x_{M})\\}$, where ${\bf
1}_{n}=[1,\ldots,1]^{T}\in\mathds{R}^{n}$. That is, $\mathcal{Q}$ denotes
constant $n$-dimensional vectors with a value in $(x_{m},x_{M})$. Because of
the sign ambiguity issue mentioned earlier, we assume that $0<x_{m}<x_{M}$. To
take the known structure of the source into account, we could maximize the
likelihood function over the set of constant signals. Let ${\bf u}=\alpha{\bf
1}_{n}$. Then, the likelihood function simplifies to

$-n\log\alpha-{1\over 2\alpha^{2}}\sum_{i=1}^{n}y_{i}^{2},$

which is again a quasi-concave function in terms of $\alpha$. Maximizing the
likelihood function leads to

$\hat{{\bf x}}_{\rm ML}=\hat{\alpha}_{\rm ML}{\bf 1}_{n},\;\;\;{\rm
where}\;\;\hat{\alpha}_{\rm ML}=({1\over n}\sum_{i=1}^{n}y_{i}^{2})^{1\over
2}.$

Given ${\bf x}=a{\bf 1}_{n}$ and ${\bf y}={\bf x}\odot{\bf w}$, we have

$\displaystyle\hat{\alpha}_{\rm ML}=\Big{(}{a^{2}\over
n}\sum_{i=1}^{n}w_{i}^{2}\Big{)}^{1\over 2}=a\Big{(}{1\over
n}\sum_{i=1}^{n}w_{i}^{2}\Big{)}^{1\over 2}.$

Now, suppose that $n$ is large. Then the calculations in the appendix show
that

$\operatorname{E}[\|\hat{{\bf x}}_{\rm ML}-{\bf
x}\|_{2}^{2}]=\operatorname{E}[n(\hat{\alpha}_{\rm
ML}-a)^{2}]=a^{2}n\operatorname{E}[(({1\over
n}\sum_{i=1}^{n}w_{i}^{2})^{1\over
2})-1)^{2}]=a^{2}n(2-2\operatorname{E}({1\over
n}\sum_{i=1}^{n}w_{i}^{2})^{1\over 2})\rightarrow\frac{a^{2}}{2}.$

Therefore, in this case, ML estimation that takes the signal structure into
account asymptotically recovers the underlying signal much more accurately
than the one that does not use the information about the signal structure.
There are a few points that we would like to emphasize here:

  1. 1.

The likelihood in the multiplicative noise problem is more complicated than
the likelihood in additive noise models. While the likelihood is quasi-convex
in the denoising problem we discussed in this section, in general the
likelihood can become non-convex in the sensing problem.

  2. 2.

As expected, given that the noise is always present in such systems, the exact
recovery is never possible. As can be seen in our simple example, the best
that can be expected from such systems is that the normalized mean square
error $\frac{1}{n}{\rm MSE}=O(\frac{1}{n})$. Note that the structure of the
signal, i.e. the fact that the signal is constant, was the major help in
reducing the MSE. Otherwise, the MSE of the MLE would be proportional to the
power of the signal as shown in (1). We will clarify the notion of structure
that will be used in this paper in the next section.

##  III Structured signals

In this paper, we study the problem of compressed sensing in the presence of
speckle noise. For ${\bf x}\in\mathds{R}^{n}$, define $X=\mathop{\rm
diag}\nolimits({\bf x})=\mathop{\rm diag}\nolimits(x_{1},\ldots,x_{n})$. In
imaging systems with the speckle noise, the measurement vector ${\bf
y}\in\mathds{R}^{m}$ is defined as

${\bf y}=AX{\bf w}+{\bf z}.$

Here, ${\bf w}\in\mathds{R}^{n}$ and ${\bf z}\in\mathds{R}^{m}$ are i.i.d.
$\mathcal{N}(0,\sigma_{w}^{2})$ and i.i.d. $\mathcal{N}(0,\sigma_{z}^{2})$,
respectively. The main goal of this paper is to study the problem of
recovering the vector ${\bf x}$ from measurements ${\bf y}$, under the
assumption that the number of measurements $m<n$. To solve this undersampled
problem, the underlying signal has to be structured. Otherwise, there are
infinitely many signals that satisfy the measurement constraints, even without
any noise in the system. Furthermore, as illustrated in the denoising example
in Section II, the structure can potentially improve the performance of the
recovery algorithms.

In this paper, we define the signal structure based on a family of compression
algorithms. This method that was introduced and developed in a series of
papers [17, 18, 19, 20], has the following advantages over the popular
sparsity structure:

  1. 1.

Compression-based compressed sensing naturally expands the scope of compressed
sensing algorithms for a class of signals from sparsity in a transform domain
to that of compression codes designed for the same class. The latter is much
richer and contains sparsity as its special case.

  2. 2.

Given that the state-of-the-art compression algorithms, such as JPEG-2000, use
sophisticated structures in the data, compression-based compressed sensing
algorithms can potentially use such sophisticated structures for the signal
recovery. Hence, compression-based algorithms that treat off-the-shelf
compression codes as blackboxes that both define and enforce the source
structure achieve very competitive performance in various applications such as
compressed sensing of natural images [18] and phase retrieval [19].

  3. 3.

In cases where the minimum required sampling rate ($m/n$) is known, ideal
compression-based compressed sensing methods theoretically achieve those
bounds [20].

In the rest of this section, we clarify the way we consider signal structure
based on a family of compression codes.

Let ${\cal Q}$ denote a compact subset of $\mathds{R}^{n}$ representing a
class of structured signals. We assume that this class is structured in the
sense that a family of compression codes indexed by their rate is known for
${\cal Q}$. More formally, we have access to a family of encoder mappings
$\mathcal{E}_{r}:\mathds{R}^{n}\to\\{1,2,\ldots,2^{nr}\\}$ and decoder
mappings $\mathcal{D}_{r}:\\{1,2,...,2^{nr}\\}\to\mathds{R}^{n}$ indexed with
rate $r$. For the signal ${\bf x}\in{\cal Q}$, $\mathcal{E}_{r}({\bf x})$
denotes its compressed version (with $nr$ bits), while
$\mathcal{D}_{r}(\mathcal{E}_{r}({\bf x}))$ denotes the reconstructed version
of ${\bf x}$ based on the $nr$ bits which were available to the compression
algorithm. To see a concrete example, we can consider ${\cal Q}$ as the class
of natural images of certain size and let $\mathcal{E}_{r}$ and
$\mathcal{D}_{r}$ denote the compression and decompression algorithms of JPEG
or JPEG2000 at a given rate $r$. Given that the compression algorithm is
usually lossy in real-world applications, we define the distortion of our
compression code as

$\delta_{r}=\sup_{{\bf x}\in\mathcal{Q}}{1\over n}\|{\bf
x}-\mathcal{D}_{r}(\mathcal{E}_{r}({\bf x}))\|^{2}.$

Naturally, it is expected that the distortion $\delta_{r}$ decreases as the
rate $r$ increases. The rate-distortion function of the given family of lossy
compression codes is defined as

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

Finally, any given encoder and decoder mappings $(\mathcal{E},\mathcal{D})$
define a codebook

$\mathcal{C}=\\{\mathcal{D}(\mathcal{E}({\bf x})):\;{\bf x}\in\mathcal{Q}\\}$

of size as most $2^{nr}$. The codebook represents all the fixed points of the
mapping $\mathcal{D}_{r}(\mathcal{E}_{r}(\cdot))$. Intuitively, according to
our compression code, the elements of the codebook are the simplest signals
because they can be exactly reconstructed from $2^{nr}$ bits.

Note that one can construct lossy compression codes for any class of signals
${\cal Q}$, for instance by simply quantizing the elements of every ${\bf
x}\in\mathcal{Q}$. However, if the class ${\cal Q}$ is structured, then for a
given compression rate $r$, there exist algorithms that achieve lower
distortion compared to those that perform element-wise quantization. The
following example further clarifies this point. Let $\mathcal{Q}$ denote
$[0,1]^{n}$. As is clear, the signals in this class do not have a particular
structure. Hence, a standard compression algorithm for this class cannot do
better than standard quantization. Suppose that we quantize each element of
${\bf x}$ at resolution $\epsilon$. It is straightforward to see that the rate
and distortion of this compression code become

$\displaystyle\delta$ $\displaystyle=$ $\displaystyle\epsilon^{2},$
$\displaystyle r$ $\displaystyle=$
$\displaystyle\frac{1}{n}\log\left(\frac{1}{\epsilon}\right)^{n}=\log\left(\frac{1}{\epsilon}\right).$
(2)

Hence, $r_{Q}(\delta)=0.5\log\left(\frac{1}{\delta}\right)$, where the
subscript $\mathcal{Q}$ only denotes that this rate is calculated for the set
$\mathcal{Q}$. (The base of all the logarithms is $2$.)
$r_{Q}(\delta)=0.5\log\left(\frac{1}{\delta}\right)$ is our baseline as it
shows the rate-distortion function of a compression algorithm on an
unstructured class of signals. To compare this rate distortion function with
an achievable rate-distortion function for a class of structured signals,
consider the set

$\tilde{\mathcal{Q}}=\\{{\bf x}\in\mathcal{Q}\ :\
x_{k+1}=x_{k+2}=\ldots=x_{n}=0\\}.$

It is straightforward to see that the same coder that we used for $Q$ offers
the following rate-distortion performance on $\tilde{\mathcal{Q}}$

$r_{\tilde{\mathcal{Q}}}(\delta)=\frac{k}{2n}\log\left(\frac{k}{n\delta}\right)$

Note that if the sparsity level $k$ is much smaller than the ambient dimension
$n$, then $\frac{r_{\tilde{\mathcal{Q}}}(\delta)}{r_{\mathcal{Q}}(\delta)}\ll
1$. In the rest of the paper, whenever we mention that a class of signals
$\mathcal{Q}$ is structured, it is assumed that
$r_{\mathcal{Q}}(\delta)\ll\log\left(\frac{1}{\delta}\right)$.

##  IV Characterizing and understanding the log likelihood function

Our goal is to recover ${\bf x}\in\mathds{R}^{n}$ from measurements

${\bf y}=AX{\bf w}+{\bf z},$ (3)

where, ${\bf w}\in\mathds{R}^{n}$ and ${\bf z}\in\mathds{R}^{m}$ are i.i.d.
$\mathcal{N}(0,\sigma_{w}^{2})$ and $\mathcal{N}(0,\sigma_{z}^{2})$,
respectively. Given the complexity of the model, in this section, we ignore
the structure of ${\bf x}$, and only study the log-likelihood function of (3).
In the next section, we will use the results we derive in this section to
obtain results for the problem of structured signal recovery.

To compute the likelihood $p({\bf y}|X;A)$, ignoring the terms that do not
depend on $X$, we have

$\displaystyle p({\bf y}|X;A)$ $\displaystyle=\int p({\bf y}|{\bf
w},X;A)p({\bf w})d{\bf w}$ $\displaystyle\propto\int\exp(-{1\over
2\sigma_{z}^{2}}\|{\bf y}-AX{\bf w}\|^{2}-{1\over 2\sigma_{w}^{2}}\|{\bf
w}\|^{2})d{\bf w}$ $\displaystyle\propto\int\exp\Big{(}-{1\over 2}{\bf
w}^{T}({1\over\sigma_{w}^{2}}I_{n}+{1\over\sigma_{z}^{2}}XA^{T}AX){\bf
w}+{1\over\sigma_{z}^{2}}{\bf y}^{T}AX{\bf w}\Big{)}d{\bf w}$
$\displaystyle\propto|\Sigma|^{1\over 2}\exp({1\over
2}\mu^{T}\Sigma^{-1}\mu),$ (4)

where $\Sigma$ and $\mu$ satisfy

$\Sigma^{-1}={1\over\sigma_{w}^{2}}I_{n}+{1\over\sigma_{z}^{2}}XA^{T}AX,$

and

$\mu^{T}\Sigma^{-1}={1\over\sigma_{z}^{2}}{\bf y}^{T}AX.$

Therefore, again ignoring the terms not depending on $X$, the log-likelihood
function can be written as

$\displaystyle 2\ell(X)$ $\displaystyle=\log\det\Sigma+\mu^{T}\Sigma^{-1}\mu$
$\displaystyle=-\log\det({1\over\sigma_{w}^{2}}I_{n}+{1\over\sigma_{z}^{2}}XA^{T}AX)+{1\over\sigma_{z}^{4}}{\bf
y}^{T}AX({1\over\sigma_{w}^{2}}I_{n}+{1\over\sigma_{z}^{2}}XA^{T}AX)^{-1}XA^{T}{\bf
y}.$ (5)

Since the main focus of this paper is on the speckle noise, we next consider
the case where the power of additive noise converges to zero ($\sigma_{z}\to
0$) and simplify the likelihood function accordingly. Depending on the number
of measurements, the limit as $\sigma_{z}\rightarrow 0$ changes. Hence, we
calculate the limit in two separate cases:

  * $\bullet$

Case I: $AA^{T}$ is invertible ($m<n$). The following theorem simplifies the
log-likelihood function for the small additive noise. Note that here we have
assumed that the number of measurements $m$ is less than the ambient dimension
$n$.

######  Theorem IV.1.

Assume that $AA^{T}$ is invertible. Then, as $\sigma_{z}\to 0$, $\ell(X)$
converges to

$\displaystyle\ell(X)$ $\displaystyle=-{1\over 2}\log\det(AX^{2}A^{T})-{1\over
2\sigma_{w}^{2}}{\bf y}^{T}(AX^{2}A^{T})^{-1}{\bf y}.$ (6)

The proof of this theorem can be found in Section VIII-B. Note that the
dependence of the log-likelihood on the signal ${\bf x}$ is not as simple as
the log-likelihood of the additive noise. This likelihood comes out of matrix
manipulations and algebra, and it is not clear why the likelihood has this
form. To get a better intuition, consider the following hypothetical
acquisition model where the measurements have the following form:

${\bf y}=AX_{o}\mathbf{1}_{n}+AX_{o}{\bf w}+{\bf z},$ (7)

The only difference here compared to the original problem is that now we have
an extra term of the form $AX_{o}\mathbf{1}_{n}$ in the measurement. We should
emphasize again that we are not aware of any real-world acquisition system for
which the model in (7) is accurate. However, as will be discussed later, this
model enables us to better understand some of the terms that appear in (6).
Using the model in (7), and following the same strategy as the one we used in
Theorem IV.1 for $\sigma_{z}\rightarrow 0$, it is straightforward to show that
the corresponding likelihood function converges to

$\displaystyle 2\ell(X)=-\log\det(AX^{2}A^{T})-\frac{1}{\sigma_{w}^{2}}({\bf
y}-A\bm{x})^{T}(AX^{2}A^{T})^{-1}({\bf y}-A\bm{x}).$ (8)

To gain some intuition behind this function, we now rederive it using some
simplifying assumptions. Given that $\sigma_{z}\rightarrow 0$, the measurement
model simplifies to ${\bf y}=AX_{o}+AX_{o}{\bf w}$. Interpreting ${\bf
u}=AX_{o}{\bf w}$ as an additive noise and assuming that (even though not
accurate) the noise is Gaussian and independent of with $AX_{o}$. It is
straightforward to see that under these assumptions, ${\bf
u}\sim\mathcal{N}(0,\sigma_{w}^{2}AX_{o}^{2}A^{T})$. Now, if we write the log-
likelihood for the model

${\bf y}=AX_{o}\mathbf{1}_{n}+{\bf u},$

with ${\bf u}\sim\mathcal{N}(0,\sigma_{w}^{2}AX^{2}A^{T})$, we obtain (8).
Note that the covariance matrix of ${\bf u}$ also depends on the signal and
should not be ignored in the likelihood. While the assumptions we have made in
our heuristic argument are not accurate, they explain why for instance we
should expect $(AX^{2}A^{T})^{-1}$ to be multiplied by our measurements. Also,
it explain where the logarithm term in our expressions comes from.

  * $\bullet$

Case II: $A^{T}A$ is invertible ($m>n$). While the main focus of this paper is
on the underdetermined settings where $m<n$, for the sake of completeness, for
completeness, we also derive the log-likelihood for the $m>n$ case. For this
case, note that multiplying both sides of ${\bf y}=AX_{o}{\bf w}$ with
$A^{T}$, it follows that $X_{o}{\bf w}=(A^{T}A)^{-1}A^{T}{\bf y}$. In other
words, in this case, we need to solve a denoising problem in the presence of a
multiplicative noise. But given $y=xw$, where
$z\sim\mathcal{N}(0,\sigma_{w}^{2})$ the likelihood of $y$ given $x$ is
${1\over\sqrt{2\pi x^{2}\sigma_{w}^{2}}}{\rm e}^{-{y^{2}\over
2x^{2}\sigma_{w}^{2}}}$. Therefore, in this case the likelihood function of
${\bf x}\in\mathds{R}^{n}$ can be derived in the following way.

######  Theorem IV.2.

Assume that $A^{T}A$ is invertible. Let ${\bf b}=(A^{T}A)^{-1}A^{T}{\bf y}$.
Then, as $\sigma_{z}\to 0$, $\ell({\bf x})$ converges to

$\displaystyle\ell({\bf x})$ $\displaystyle=-\sum_{i=1}^{n}\Big{(}{1\over
2}\log x_{i}^{2}+{b_{i}^{2}\over 2\sigma_{w}^{2}x_{i}^{2}}\Big{)}.$ (9)

While we derived this formula by converting the inverse problem to a denoising
problem, one can start with (5) and simplify it to reach the same conclusion.

##  V Compressed sensing in the presence of speckle noise

###  V-A Our proposed recovery optimization problem

Consider recovering ${\bf x}_{o}\in\mathcal{Q}\subset\mathds{R}^{n}$ from
measurements ${\bf y}=AX_{o}{\bf w}$, where $X_{o}=\mathop{\rm
diag}\nolimits({\bf x}_{o})$ and $m<n$. Also, assume that ${\bf
x}_{o}\in\mathcal{Q}\subset\mathds{R}^{n}$, where $\mathcal{Q}$ is a set of
structured signals for which we have a family of compression codes indexed
with rate $r$ with corresponding codebook $\mathcal{C}_{r}$. Then, employing
Theorem IV.1, when $m<n$, we use the following optimization problem to obtain
an estimate of ${\bf x}_{o}$:

$\displaystyle\hat{{\bf
x}}=\operatorname*{arg\,min}_{\begin{subarray}{c}X=\mathop{\rm
diag}\nolimits({\bf x}),\\\ {\bf
x}\in\mathcal{C}_{r}\end{subarray}}\Big{[}\log\det(AX^{2}A^{T})+{1\over\sigma_{w}^{2}}{\bf
y}^{T}(AX^{2}A^{T})^{-1}{\bf y}\Big{]}.$ (10)

In other words, instead of minimizing the negative log-likelihood over all
signals, we only focus on the signals that have simple representation
according to our compression code, i.e. they can be represented exactly with
$2^{nr}$ bits. Given that we have access to a family of compression codes, the
rate $r$ can be considered as a free parameter that the user can tune to
obtain better performance. We will discuss this parameter in more details
later.

Note that this optimization problem (10) is hard not only because the fact
that the negative log-likelihood is non-convex, but also because the set on
which we would like to solve our optimization problem is discrete and very
large. We will later discuss practical approaches to approximate the solution
of this optimization problem.

###  V-B Main theoretical result

As described in Section $\bullet$ ‣ IV, when the number of measurements
exceeds the ambient dimension of the problem, i.e., $m\geq n$, the problem
reduces to a denoising problem. Therefore, we focus on the case of $m<n$.
Given $A\in\mathds{R}^{m\times n}$, define mapping
$\Sigma:\mathds{R}^{n}\to\mathds{R}^{m\times m}$ as

$\displaystyle\Sigma({\bf x})=(AX^{2}A^{T})^{-1},$ (11)

if the inverse is well-defined. (Here $X=\mathop{\rm diag}\nolimits({{\bf
x}})$.) Also, define function $\mathds{R}^{m\times m}\to\mathds{R}$, as

$f(\Sigma)=-\log\det\Sigma+{1\over\sigma_{w}^{2}}{\rm Tr}(\Sigma{\bf y}{\bf
y}^{T}).$

Let ${\bf y}=AX_{o}{\bf w}$, where $X_{o}=\mathop{\rm diag}\nolimits({{\bf
x}_{o}})$. Consider a rate-$r$ distortion-$\delta$ compression code with
codebook $\mathcal{C}_{r}$. As argued in Section $\bullet$ ‣ IV, an ML-based
recovery algorithm that takes advantage of the given compression code solves
the following optimization problem:

$\displaystyle\hat{{\bf
x}}_{o}=\operatorname*{arg\,min}_{\begin{subarray}{c}X=\mathop{\rm
diag}\nolimits({{\bf x}}),\\\ {\bf
x}\in\mathcal{C}_{r}\end{subarray}}f(\Sigma({\bf x})).$ (12)

The main theoretical result of this paper is the following theorem, which
shows that given sufficient number of measurements which is characterized in
terms of the properties of the compression code, the described ML-based
recovery method is able to recover ${\bf x}$ from the measurements. To avoid
sign ambiguity issue we discussed in Section II we assume that the elements of
${\bf x}_{o}$ are all positive.

######  Theorem V.1.

Let ${\bf y}=AX_{o}{\bf w}$, where $X_{o}=\mathop{\rm diag}\nolimits({{\bf
x}_{o}})$. Let $\hat{{\bf x}}_{o}$ denote the solution of (12). Assume that
$m<n/4$ and define $\gamma\triangleq({1+2\sqrt{m/n}\over 1-2\sqrt{m/n}})^{2}$,
$\alpha\triangleq{x_{\max}\over x_{\min}}$, $\rho_{1}\triangleq
4\sqrt{2}\alpha^{8}\gamma^{5}(1+2\alpha^{2}\gamma)^{2}(1+2\sqrt{m/n})^{2}$ and
$\rho_{2}\triangleq(1+2\alpha^{2}\gamma)^{2}\gamma^{7}\alpha^{14}$. Then, for
any $\epsilon>0$,

$\displaystyle{1\over n}\|{\bf x}_{o}-\hat{{\bf
x}}_{o}\|_{2}^{2}\leq\rho_{1}\sqrt{(1+\epsilon)nr\over
m}+\rho_{2}x_{\max}^{2}\delta,$ (13)

with a probability higher than

$1-n{\rm e}^{-0.09m}-n{\rm e}^{-0.84m}-2^{-nr\epsilon+1}-2{\rm e}^{-{m\over
2}}.$

Note that $r$ and $\delta$ denote the rate and distortion of the compression
algorithm used in (12).

We would like to make a few remarks about this result.

######  Remark 1.

Unlike the classical compressed sensing problem, here the exact recovery is
not possible even though there is no additive noise in the measurements. This
is due to two issues:

  1. 1.

The existence of multiplicative noise. Given the multiplicative noise that is
present in the system, it should be clear that the exact recovery is not
possible.

  2. 2.

The distortion in our compression algorithm: Given that our search space is
only on the codewords, we should not expect to be able to get to the exact
solution, because the exact solution may not be even a codeword. We expect
this effect to diminish as $\delta\rightarrow 0$. As is clear the term
$\rho_{2}x_{\max}^{2}\delta\rightarrow 0$ as $\delta\rightarrow 0$. However,
note that as we let $\delta\rightarrow 0$, we most probably are letting
$r(\delta)$ grow to infinity, which in turn will blow up the first term in the
MSE, i.e. $\rho_{1}\sqrt{(1+\epsilon)nr\over m}$. Hence, the best choice of
$\delta$ depends on the rate-distortion function. In practice, one may use a
cross-validation technique for finding the best choice of $\delta$.

######  Remark 2.

The condition $m<n/4$ does not allow us to let $m$ go to infinity without
bound in our theorem. Note that as we discussed in Section IV, in order to
derive the optimization problem (12) we assumed that $m<n$. Hence, it is
expected that solving (12) does not generate any reliable result for $m>n$.
Furthermore, as $m$ gets close to $n$, the matrices $AX^{2}A^{T}$ start to
have eigenvalues close to zero for many signals $X^{2}$. This seems to have
some adverse algorithmic and statistical effect on the recovery of ${\bf
x}_{o}$. While in practice, we would always like to employ these systems in
highly underdetermined regime, still a better understanding of the problem
when $n/4<m<n$ will shed more light on the landscape and shape of the
likelihood function in multiplicative noise systems.

Before we proceed to the discussion of the algorithmic issues of solving (12),
we would like to clarify the statement of this theorem through a well-known
example. In this example we assume that the family of compression codes on
$\mathcal{Q}$ satisfy
$r(\delta)\leq\frac{k}{n}\log\frac{1}{\delta}+\frac{k}{n}\log(n)$. For
instance, it is straightforward to construct such a family compression codes
for the class of $k$-sparse signals in unit sphere, or class of piecewise
constant signals with $k$ jumps.

######  Corollary 1.

Suppose that the family of compression algorithms satisfy
$r(\delta)\leq\frac{k}{n}\log\frac{1}{\delta}+\frac{k}{n}\log n$. For any
$\epsilon>0$, if we use the compression code with distortion
$\delta=\frac{1}{n}$, then we have (the rest of the notations are the same as
the ones introduced in Theorem V.1)

$\displaystyle{1\over n}\|{\bf x}_{o}-\hat{{\bf
x}}_{o}\|_{2}^{2}\leq\rho_{1}\sqrt{2(1+\epsilon)k\log n\over
m}+\rho_{2}x_{\max}^{2}\frac{1}{n},$ (14)

with a probability higher than

$1-n{\rm e}^{-0.09m}-n{\rm e}^{-0.84m}-2^{-2\epsilon k\log n+1}-2{\rm
e}^{-{m\over 2}}.$

Note that given the fact that $m<\frac{n}{4}$ the dominant term in the right
hand side of (14) is the first term. The term $\frac{k\log n}{m}$ is the term
that appears in the compressed sensing problem with additive noise as well
[21]. Note however that, in the additive noise case, instead of having
$\sqrt{k\log n\over m}$, we have ${k\log n\over m}$. Whether this is an
artifact of our proof technique or a fundamental difficulty of the
multiplicative noise problem is an important open problem for our future
research.

##  VI Recovery algorithms

###  VI-A Roadmap

In Section IV, we showed that a compression-based ML recovery method recovers
signal ${\bf x}_{o}\in\mathcal{Q}$ from measurements ${\bf y}=AX{\bf w}_{o}$,
where $X_{o}=\mathop{\rm diag}\nolimits({{\bf x}_{o}})=\mathop{\rm
diag}\nolimits({x_{o,1}},\ldots,{x_{o,n}})$, by solving the optimization
described in (10). However, (10) is a challenging optimization as the cost
function is a high-dimensional non-convex function of ${\bf x}$ over an
exponentially large discrete set. Therefore, obtaining the solution of (10) is
not straightforward.

As an initial step in solving this optimization problem, we propose two
different algorithms which attempt to solve (10): (1) a projected gradient
descent algorithm, and (2) a gradient-free multi-level optimization algorithm.
We describe the derivations of the algorithms in this section, and the
simulation results in section VII.

###  VI-B Projected Gradient Descent Algorithm

####  VI-B1 General algorithm

Projected gradient descent (PGD) algorithm and its close relative proximal
gradient descent [22] are among the most popular algorithms in the fields of
signal processing, machine learning, and optimization. In each step, PGD takes
a small step in the direction of the gradient of the cost function (ignoring
all the constraints on the solution) and then project the current estimate on
the constraint set. These algorithms are known to converge to a global
minimizer of convex optimization problems (if the step size is picked
according to certain rules). PGD algorithms have also been applied in practice
to non-convex optimization problems. While they have exhibit a good
performance for non-convex problems as well, there is still no general theory
that can explain the good performance of PGD for a large class of non-convex
problems. However, researchers have proved the success of such algorithms in
particular instances of non-convex optimization problems [23, 18, 19]. In this
paper, we derive the PGD algorithm for optimization problem (10). We also
mention some details about how we implement the algorithm, and evaluate its
performance through simulations. We leave the theoretical evaluation of this
algorithm for a future research. As we discussed before, we are interested in
solving the following optimization problem with PGD:

$\displaystyle\hat{{\bf
x}}=\operatorname*{arg\,min}_{\begin{subarray}{c}X=\mathop{\rm
diag}\nolimits({\bf x}),\\\ {\bf
x}\in\mathcal{C}_{r}\end{subarray}}\Big{[}\log\det(AX^{2}A^{T})+{1\over\sigma_{w}^{2}}{\bf
y}^{T}(AX^{2}A^{T})^{-1}{\bf y}\Big{]}.$

The PGD algorithm has two steps: (i) moving in the direction of the gradient
of $\log\det(AX^{2}A^{T})+{1\over\sigma_{w}^{2}}{\bf
y}^{T}(AX^{2}A^{T})^{-1}{\bf y}$, and (ii) projecting onto the constraint set
$\mathcal{C}_{r}$. We discuss each step and some other details below.

  * $\bullet$

Gradient calculation: To derive the PGD algorithm, we first need to compute
the gradient of the cost function $f(\Sigma({\bf
x}))=\log\det(AX^{2}A^{T})+{1\over\sigma_{w}^{2}}{\bf
y}^{T}(AX^{2}A^{T})^{-1}{\bf y}$ with respect to ${\bf x}$. Note that because
$X$ is a diagonal matrix, $AX^{2}A^{T}=\sum_{i=1}^{n}x_{i}^{2}{\bf a}_{i}{\bf
a}_{i}^{T}$, where ${\bf a}_{i}$ is the $i$-th column of matrix $A$. Define
${\bf u}={\bf x}^{2}$, i.e., $u_{i}=x_{i}^{2}$. We first compute the partial
derivative of $\underline{f}({\bf u})=\log\det(\sum_{i=1}^{n}u_{i}{\bf
a}_{i}{\bf a}_{i}^{T})+{1\over\sigma_{w}^{2}}{\bf
y}^{T}(\sum_{i=1}^{n}u_{i}{\bf a}_{i}{\bf a}_{i}^{T})^{-1}{\bf
y}=f(\Sigma({\bf x}))$ with respect to $u_{i}$. Let ${\bf e}_{i}$ denote the
unit vector in direction $i$. Let $B=\sum_{j=1}^{n}u_{j}{\bf a}_{j}{\bf
a}_{j}^{T}$. Then,

$\displaystyle\underline{f}({\bf u}+\delta u_{i}{\bf
e}_{i})-\underline{f}({\bf u})$ $\displaystyle=\log\det(B+\delta_{i}{\bf
a}_{i}{\bf a}_{i}^{T})+{1\over\sigma_{w}^{2}}{\bf y}^{T}(B+\delta_{i}{\bf
a}_{i}{\bf a}_{i}^{T})^{-1}{\bf y}-\log\det B+{1\over\sigma_{w}^{2}}{\bf
y}^{T}B^{-1}{\bf y}$ $\displaystyle=\log\det(I_{m}+\delta_{i}B^{-{1\over
2}}{\bf a}_{i}{\bf a}_{i}^{T}B^{-{1\over 2}})+{1\over\sigma_{w}^{2}}{\bf
y}^{T}((B+\delta_{i}{\bf a}_{i}{\bf a}_{i}^{T})^{-1}-B^{-1}){\bf y}.$ (15)

Using the Woodbury matrix identity,

$(B+\delta_{i}{\bf a}_{i}{\bf a}_{i}^{T})^{-1}=B^{-1}-{\delta_{i}B^{-1}{\bf
a}_{i}{\bf a}_{i}^{T}B^{-1}\over 1+\delta_{i}{\bf a}_{i}^{T}B^{-1}{\bf
a}_{j}}.$

Also, given that all but one of the eigenvalues of
$I_{m}+\delta_{i}B^{-{1\over 2}}{\bf a}_{i}{\bf a}_{i}^{T}B^{-{1\over 2}}$ are
equal to one, we have $\log\det(I_{m}+\delta_{i}B^{-{1\over 2}}{\bf a}_{i}{\bf
a}_{i}^{T}B^{-{1\over 2}})=\log(1+\delta_{i}{\bf a}_{i}^{T}B^{-1}{\bf
a}_{i})$. Using these identities, it follows that ${\partial\underline{f}({\bf
u})\over\partial u_{i}}={\bf a}_{i}^{T}B^{-1}{\bf a}_{i}-({\bf
a}_{i}^{T}B^{-1}{\bf y})^{2}$, and

${\partial f(\Sigma({\bf x}))\over\partial x_{i}}=2x_{i}({\bf
a}_{i}^{T}B^{-1}{\bf a}_{i}-({\bf a}_{i}^{T}B^{-1}{\bf y})^{2}),$ (16)

where $B=AX^{2}A^{T}$. Equation (16) gives us the gradients that are necessary
for the PGD algorithm.

  * $\bullet$

Projection, proper step-size, and initialization: As discussed before, PGD
algorithm requires the projection onto $\mathcal{C}_{r}$ defined as

$\pi_{\mathcal{C}_{r}}({\bf u})=\operatorname*{arg\,min}_{{\bf
u}^{\prime}\in\mathcal{C}}\lVert{\bf u}^{\prime}-{\bf u}\rVert.$

Note that, because $\mathcal{C}_{r}$ is often a very large non-convex set,
this projection can be computationally demanding. While in some cases, such as
the one we will mention in the next section, the projection can be calculated
efficiently and accurately using for instance dynamic programming, in many
other examples this is not the case. One successful approximation for
$\pi_{\mathcal{C}_{r}}(u)$ that has shown promising results in other
applications (e.g. [18, 19]) is

$\hat{\pi}_{\mathcal{C}_{r}}({\bf u})=\mathcal{D}_{r}(\mathcal{E}_{r}({\bf
u}))$

where $\mathcal{E}_{r}$ and $\mathcal{D}_{r}$ are the encoder and decoder of
the compression algorithms respectively. There are two reasons for using this
approximation for the projection operator: (i) All the state-of-the-art
compression algorithms have computationally efficient decoders and encoders.
Hence, the approximate projection will be fast as well. (ii) All state-of-the-
art compression algorithms try to make sure that they are projecting each data
point onto the closest codeword or a codeword that is in the vicinity of the
closest codeword. Hence, this approximation is expected to be quite accurate.

Based on our discussion so far, our PGD algorithm proceeds according to the
following iteration:

${\bf x}_{t}=\pi_{\mathcal{C}_{r}}({\bf x}_{t-1}-\mu_{t}\bm{g}_{t}),$ (17)

where ${\bf x}_{t}$ is the estimate of ${\bf x}_{o}$ at iteration $t$,
$\bm{g}_{t}$ is the gradient of
$\log\det(AX^{2}A^{T})+{1\over\sigma_{w}^{2}}{\bf y}^{T}(AX^{2}A^{T})^{-1}{\bf
y}$ based on (16), and $\mu_{t}$ is the step-size at iteration $t$. There are
two remaining ingredients in our algorithm: (i) the choice of the step-size
$\mu_{t}$, and (ii) the choice of initialization. To set the step size we use
the line search to find a value of $\mu_{t}$ that makes
$f(\Sigma(\pi_{\mathcal{C}}({\bf x}_{t-1}-\mu\bm{g}_{t})))$ smaller than
$f(\Sigma({\bf x}_{t-1}))$. Additionally, we initialize the algorithm with
${\bf x}_{0}$ being a constant vector. Based on all our discussions, the PGD-
based algorithm operates as described in Algorithm 1. As discussed before,
whenever possible we will use the exact choice of the projection function.
Otherwise, we will use the approximation presented in (17).

Result: $\hat{{\bf x}}_{T}$

initialize ${\bf x}_{0}\in\mathds{R}^{n}$

$X_{0}=\mathop{\rm diag}\nolimits({\bf x}_{0})$,

$B_{0}=AX_{0}^{2}A^{T}$;

for  _$t=1,\ldots,T$_ do

for  _$i=1,\ldots,n$_ do

$s_{t,i}=x_{t-1,i}-\mu x_{t-1,i}({\bf a}_{i}^{T}B_{t-1}^{-1}{\bf a}_{i}-({\bf
a}_{i}^{T}B_{t-1}^{-1}{\bf y})^{2})$

end for

${\bf x}_{t}=\pi_{\mathcal{C}_{r}}({\bf s}_{t})$

$X_{t}=\mathop{\rm diag}\nolimits({\bf x}_{t})$

$B_{t}=AX_{t}^{2}A^{T}$

end for

Algorithm 1 PGD-based recovery from under-determined measurements corrupted by
multiplicative noise

####  VI-B2 Piecewise constant functions

To provide a concrete example for the discussions of the previous section, we
mention a popular class of functions in imaging systems, i.e. the class of
piecewise constant functions. Refer to [24] for applications of this model
beyond imaging. Generalizations of our discussion to the class of piecewise
polynomial functions are straightforward. However, to keep the notation
simple, we focus on piecewise constant functions here.

Let $\mathcal{Q}^{J,n}(x_{\min},x_{\max})$ denote the set of piecewise
constant signals in $\mathds{R}^{n}$ with the maximum of $J$ jumps and values
bounded between $x_{\min}$ and $x_{\max}$, where $0<x_{\min}\leq x_{\max}$. To
apply Alg. 1, we need to have a compression code for this class of signals. In
our simulations, we use a compression code that operates as follows.

  * $\bullet$

Encoder $\mathcal{E}_{n}$: Consider ${\bf
x}\in\mathcal{Q}^{J,n}(x_{\min},x_{\max})$. Let $i_{1},\ldots,i_{j}$ denote
the location of the jumps in ${\bf x}$. That is, $x_{i_{k}}\neq x_{i_{k}+1}$.
(Note that $j\leq J$.) The encoder describes the number of jumps and the jump
locations (using $(j+1)\lceil\log_{2}J\rceil$ bits). Then, it describes the
$b$-bit quantized versions of the $j+1$ values corresponding to the $j+1$
constant intervals (overall using
$(j+1)(\lceil\log_{2}(x_{\max}-x_{\min})\rceil+b)$ bits).

  * $\bullet$

Decoder $\mathcal{D}_{n}$: The decoder receives the number of jumps, locations
of the jumps and the $b$-bit quantized values and reconstructs the signal
accordingly.

Let $\mathcal{C}_{r}$ denote the codebook corresponding to this compression
code. Note that both the rate and the distortion depend on $J$ and $b$. For
this class of functions (and also the class of piecewise polynomial
functions), and for the codewords we have constructed, we can implement the
exact projection function $\pi_{\mathcal{C}_{r}}({\bf u})$ efficiently using
dynamic programming. Consider ${\bf c}\in\mathcal{C}_{r}$, with $j$ jumps at
$i_{1}<\ldots<i_{j}$ and values $a_{1},\ldots,a_{j+1}\in\mathcal{X}_{m}$,
where

$\mathcal{X}_{b}\triangleq\\{[x]_{b}:\;x\in(x_{\min},x_{\max})\\}.$

Let $i_{0}\triangleq 0$ and $i_{j+1}\triangleq n$. Then,

$\displaystyle\|{\bf x}-{\bf
c}\|^{2}=\sum_{k=0}^{j}\sum_{i=i_{k}+1}^{i_{k+1}}(a_{k+1}-x_{i})^{2}.$ (18)

Fixing the jump locations, it is straightforward to find the minimizer of
(18), as

$\min_{(a_{1},\ldots,a_{j+1})\in\mathcal{X}_{b}^{j+1}}\sum_{k=0}^{j}\sum_{i=i_{k}+1}^{i_{k+1}}(a_{k+1}-x_{i})^{2}=\sum_{k=0}^{j}\min_{a_{k+1}\in\mathcal{X}_{b}}\sum_{i=i_{k}+1}^{i_{k+1}}(a_{k+1}-x_{i})^{2}.$

Therefore, the complexity of minimizing (18) lies in finding jumps’ locations.
But, the expression of (18) suggests that we can employ the Viterbi (dynamic
programming) algorithm [25] to find the jumps’ locations. Consider a Trellis
diagram with $J+1$ layers and $n$ nodes at each layer. In this graph, layer
$k$ corresponds to the $k$-th jump and node $i$ in layer $k$ correspond to
having the $k$-th jump at $i$. Using this construction, node $i$ in layer $k$
only connects to nodes $i+1,\ldots,n$ at layer $k+1$. The weight assigned to
the edge connection node $i_{1}$ in layer $k$ to node $i_{2}$ in layer $k+1$
is equal to

$\min_{a\in\mathcal{X}_{b}}\sum_{i=i_{1}+1}^{i_{2}}(a-x_{i})^{2}.$

Applying the Viterbi algorithm on the described Trellis diagram yields the
solution of (18).

###  VI-C Multilevel formulation

In addition to the projected gradient descent algorithm proposed in the
previous section, in some cases, such as the piecewise polynomial functions,
we can reformulate our problem as a multilevel optimization problem. Namely,
we split the optimization into a discrete component (optimization of the
structure), and a continuous (or near continuous) component (optimization of
the values given the structure). To avoid introducing new notations, we
describe this method on the concrete example of piecewise constant functions
that we described in the previous section.

Given integers $1=d_{1},d_{2},\dotsc,d_{k+1}=n+1$ and positive real values
$\theta_{1},\dotsc,\theta_{k}$, we may define the following piecewise-constant
signal ${\bf x}(\bm{\theta},{\bf d})$:

${\bf x}(\bm{\theta},{\bf d})_{i}=\theta_{l}\quad\text{ if }d_{l}\leq
i<d_{l+1}.$

We may thus view the recovery problem as:

$\operatorname*{arg\,min}_{{\bf
d}\in\mathbb{N}^{k+1}}\min_{\bm{\theta}\in\mathbb{R}^{k}}f(\Sigma({\bf
x}(\bm{\theta},{\bf d}))),$

Now, note that the inner optimization problem, which we denote as:

$h({\bf d})=\min_{\bm{\theta}\in\mathbb{R}^{k}}f(\Sigma({\bf
x}(\bm{\theta},{\bf d})))$

can be solved using traditional gradient-descent or quasi-Newton type solvers,
such as L-BFGS. Our simulation results show that the non-convexity of the cost
does not present any major challenges to these methods. Hence, suppose that
once ${\bf d}$ is given the calculation of the optimal solution is
straightforward. We then attempt to optimize $h({\bf d})$ directly as a
function of ${\bf d}$. However, we note here that $h$ is a multivariate
function of discrete quantities, and gradient information (or similar) is not
easily computable for $h$. Instead, we propose to use a gradient-free method
to optimize $h({\bf d})$. We leverage recent advances in hyper-parameter
optimization methods, and make use of _optuna_ [26], a general purpose
gradient-free optimizer, to optimize $h$. Note that using alternating
minimization on ${\bf d}$ and $\bm{\theta}$, we can hope to obtain an
approximate minimizer of the $f(\Sigma({\bf x}(\bm{\theta},{\bf d})))$.

##  VII Simulation results

As a preliminary investigation, we analyze the performance of the proposed
methods on a simulated example. To evaluate the performance of our algorithms,
we consider the problem of recovering a piecewise-constant signal from its
undersampled measurements. We study the performance of the following
algorithms that are based on the two algorithms we discussed in the last
section:

  1. 1.

A vanilla implementation of Algorithm 1, where the step size $\mu$ is selected
by line search at every iteration. We refer to this implementation as _pgd_.

  2. 2.

We noticed a dependence of the optimization quality on the choice of initial
value $x_{0}$. We consider an algorithm which performs Algorithm 1 for a set
of initializations (namely all constant signals with a given magnitude), and
selects the output among those which maximizes the likelihood of the data. We
refer to this implementation as _pgd + init_.

  3. 3.

We also consider the multi-level formulation described in section VI-C. We
refer to this method as _multilevel_.

  4. 4.

Additionally, we present a variant of this multi-level approach, where
approximate break locations are first estimated using a PGD method, then
neighborhoods of those breaks are searched using a gradient-free optimizer. We
refer to this implementation as _pgd + multilevel_.

Figure 1: The piecewise constant signal used in the simulations.

In our simulations, we consider a signal with $5$ pieces in dimension
$n=1000$. This signal is shown in Figure 1. We evaluate the performance of our
algorithms for $m\in\\{200,250,300,350,400,450,500\\}$ measurements. For each
$m$, we generate 50 random measurement matrices and noise vectors and measure
the reconstruction error of each of the four algorithms we mentioned above. To
measure the amount of reconstruction error of vector $\hat{{\bf x}}$ we use
the PSNR in decibels defined as

$10\log_{10}\left(\frac{\|{\bf x}\|^{2}_{\infty}}{\|\hat{{\bf x}}-{\bf
x}\|_{2}^{2}}\right)$

A few instances of reconstructions are shown in Figures 2 and 3. A summary of
our simulation results is presented in Figure 4. We can draw the following
conclusions based on these simulations:

Figure 2: Signal ${\bf x}_{o}$ shown in Figure 1 and its reconstruction from
its undersampled measurements with multiplicative noise using the four
algorithms described in this section. The number of measurements is set to
$m=250$. Figure 3: Signal ${\bf x}_{o}$ shown in Figure 1 and its
reconstruction from its undersampled measurements with multiplicative noise
using the four algorithms described in this section. The number of
measurements is set to $m=400$.

  1. 1.

As is clear from Figure 4, pgd+multilevel algorithm overall offers a better
performance than the other algorithms. However, we should note that multilevel
algorithms are often more computationally demanding (see Table I), and face
scalability issues to larger structures (e.g. larger number of constant
pieces), or more complex structures (e.g. those that may be found in images).
Finding algorithms that are more efficient than pgd-based methods and are
still scalable for more sophisticated signals and structures is an important
direction for our future research.

  2. 2.

As can be seen in Figure 4(b), in general the true signal has higher log-
likelihood than the solutions our algorithms are converging to. This is an
indication of the fact that our algorithms have not been able to find the
global minimizer of the negative log-likelihood. This again poses an open
algorithmic problem. Can we find algorithms that achieve better solutions?
Again this is an important direction for our future research.

(a) Statistical performance (b) Optimization performance

Figure 4: Summary of our simulation results on a piecewise constant signal with $5$ pieces. Figure (a) exhibits the PSNR of the reconstruction of the four algorithms based on $50$ simulations as a function of the number of measurements $m$. In each simulation a fresh measurement matrix and random noise vector are drawn. Note that the bars represent the ninety percent confidence intervals. Figure (b) shows the negative log-likelihood of the solution of our algorithms. For comparison purposes we have also included the negative log-likelihood of ${\bf x}_{o}$ the true signal. The label for the negative log-likelihood of the true signal is data in Figure (b).  $m$ | Method | Time (s) | Evaluations
---|---|---|---
200 | PGD | 0.6 (0.2) | 89 (25)
| PGD + init | 8.8 (1.7) | 1325 (223)
| Multilevel | 33.6 (1.2) | 4492 (55)
| PGD + multilevel | 33.5 (1.6) | 4454 (77)
400 | PGD | 1.1 (0.3) | 65 (16)
| PGD + init | 20.6 (3.0) | 1173 (161)
| Multilevel | 99.4 (5.7) | 4332 (35)
| PGD + multilevel | 94.3 (15.7) | 4302 (36)
Table I: Time and number of likelihood or gradient evaluations performed by
the proposed methods: average and standard deviations are calculated based on
$50$ Monte Carlo simulations. Note that the signal is piecewise constant with
$5$ pieces in dimension $n=1000$. The number of measurements are mentioned in
the first column.

##  VIII Proofs

###  VIII-A Preliminaries

Before stating the proofs, we present some lemmas, some new and some known,
that are going to be used in the proofs on the main results.

######  Lemma VIII.1.

[27] Let the elements of an $m\times n$ (m<n) matrix $A$ be drawn
independently from $\mathcal{N}(0,1)$. Then for any $t>0$ we have

$\operatorname{P}(\sqrt{n}-\sqrt{m}-t\leq\sigma_{\min}(A)\leq\sigma_{\max}(A)\leq\sqrt{n}+\sqrt{m}+t)\geq
1-2{\rm e}^{-\frac{t^{2}}{2}}.$

######  Lemma VIII.2 (Concentration of $\chi^{2}$ [28]).

Let $Z_{1},Z_{2},\ldots,Z_{n}$ denote a sequence of independent
$\mathcal{N}(0,1)$ random variables. Then, for any $t\in(0,1)$, we have

$\operatorname{P}(\sum_{i=1}^{n}Z_{i}^{2}\leq m(1-t))={\rm
e}^{\frac{m}{2}(t+\log(1-t))}.$

Also, for any $t>0$,

$\operatorname{P}(\sum_{i=1}^{n}Z_{i}^{2}\geq m(1+t))={\rm
e}^{-\frac{m}{2}(t-\log(1+t))}.$

Define

$\|X\|_{\psi_{2}}=\inf\\{t>0:\mathbb{E}(\exp(X^{2}/t^{2}))\leq 2\\}.$ (19)

######  Theorem VIII.3 (Hanson-Wright inequality).

Let ${\bf X}=(X_{1},...,X_{n})$ be a random vector with independent components
with $\operatorname{E}[X_{i}]=0$ and $\|X_{i}\|_{\Psi_{2}}\leq K$. Let A be an
$n\times n$ matrix. Then, for $t>0$,

$\operatorname{P}\Big{(}|{\bf X}^{T}A{\bf X}-\operatorname{E}[{\bf X}^{T}A{\bf
X}]|>t\Big{)}\leq 2\exp\Big{(}-c\min({t^{2}\over K^{4}\|A\|_{\rm
HS}^{2}},{t\over K^{2}\|A\|})\Big{)}.$

######  Lemma VIII.4.

Let $\lambda_{\max}$ denote the maximum eigenvalue of $\Sigma_{o}^{-{1\over
2}}\Delta\Sigma\Sigma_{o}^{-{1\over 2}}$. Then,

$\displaystyle{1\over 2(1+\lambda_{\max})^{2}}{\rm
Tr}(\Sigma_{o}^{-1}\Delta\Sigma\Sigma_{o}^{-1}\Delta\Sigma)\leq\bar{f}(\Sigma)-\bar{f}(\Sigma_{o})$
$\displaystyle\leq{1\over 2}{\rm
Tr}(\Sigma_{o}^{-1}\Delta\Sigma\Sigma_{o}^{-1}\Delta\Sigma),$ (20)

where $\Delta\Sigma=\Sigma-\Sigma_{o}$.

###### Proof of Lemma VIII.4.

Note that

$\displaystyle\bar{f}(\Sigma)-\bar{f}(\Sigma_{o})$
$\displaystyle=-\log\det\Sigma+{\rm Tr}(\Sigma
AX_{o}^{2}A^{T})+\log\det\Sigma_{o}-{\rm Tr}(\Sigma_{o}AX_{o}^{2}A^{T}).$ (21)

But,

$\displaystyle\log\det\Sigma=\log\det(\Sigma_{o}+\Delta\Sigma)=\log\det\Sigma_{o}^{1\over
2}(I+\Sigma_{o}^{-{1\over 2}}\Delta\Sigma\Sigma_{o}^{-{1\over
2}})\Sigma_{o}^{1\over 2}.$ (22)

Therefore,

$\displaystyle\log\det\Sigma-\log\det\Sigma_{o}$
$\displaystyle=\log\det(\Sigma_{o}+\Delta\Sigma)=\log\det(I+\Sigma_{o}^{-{1\over
2}}\Delta\Sigma\Sigma_{o}^{-{1\over 2}})$
$\displaystyle=\sum_{i=1}^{n}\log(1+\lambda_{i}(\Sigma_{o}^{-{1\over
2}}\Delta\Sigma\Sigma_{o}^{-{1\over 2}})),$ (23)

where $\lambda_{i}=\lambda_{i}(\Sigma_{o}^{-{1\over
2}}\Delta\Sigma\Sigma_{o}^{-{1\over 2}})$, $i=1,\ldots,n$. Moreover, by the
mean value theorem, for $x>0$, $\log(1+x)=x-{1\over 2(1+a)^{2}}x^{2}$, for
some $a\in(0,x)$. Therefore, since $\lambda_{i}\geq 0$, for all $i$, we have

$\displaystyle\log\det\Sigma-\log\det\Sigma_{o}$
$\displaystyle\leq\sum_{i=1}^{n}\lambda_{i}-{1\over
2(1+\lambda_{\max})^{2}}\sum_{i=1}^{n}\lambda_{i}^{2}$ $\displaystyle={\rm
Tr}(\Sigma_{o}^{-1}\Delta\Sigma)-{1\over 2(1+\lambda_{\max})^{2}}{\rm
Tr}(\Sigma_{o}^{-1}\Delta\Sigma\Sigma_{o}^{-1}\Delta\Sigma).$ (24)

Combining this bound with (21), the desired result follows

$\displaystyle\bar{f}(\Sigma)-\bar{f}(\Sigma_{o})$ $\displaystyle\geq-{\rm
Tr}(\Sigma_{o}^{-1}\Delta\Sigma)+{1\over 2(1+\lambda_{\max})^{2}}{\rm
Tr}(\Sigma_{o}^{-1}\Delta\Sigma\Sigma_{o}^{-1}\Delta\Sigma)+{\rm Tr}(\Sigma
AX_{o}^{2}A^{T})-{\rm Tr}(\Sigma_{o}AX_{o}^{2}A^{T})$ $\displaystyle={1\over
2(1+\lambda_{\max})^{2}}{\rm
Tr}(\Sigma_{o}^{-1}\Delta\Sigma\Sigma_{o}^{-1}\Delta\Sigma).$ (25)

Note that to obtain an upper bound for $\bar{f}(\Sigma)-\bar{f}(\Sigma_{o})$
we should replace $\lambda_{\max}$ with zero in (VIII-A).

∎

######  Lemma VIII.5.

Given $\Sigma=(AX^{2}A^{T})^{-1}$, let $\delta
f({\Sigma})={f}({\Sigma})-\bar{f}({\Sigma})$. Then, for $t>0$,

$\operatorname{P}(|\delta f(\Sigma)|\geq t|A)\leq
2\exp\Big{(}-{ct\min(1,{t\over 2m})\over 2x_{\max}^{4}\|A^{T}\Sigma
A\|^{2}}\Big{)}.$

Also, $\|A^{T}\Sigma
A\|^{2}\leq{\lambda^{2}_{\max}(A^{T}A)\over\lambda^{2}_{\min}(AA^{T})x^{4}_{\min}}$.

###### Proof of Lemma VIII.5.

By definition,

$\delta f({\Sigma})={f}({\Sigma})-\bar{f}({\Sigma})={1\over\sigma_{w}^{2}}{\bf
w}^{T}X_{o}A^{T}\Sigma AX_{o}{\bf w}-{\rm Tr}(\Sigma AX_{o}^{2}A^{T})$

Define matrix $B\in\mathds{R}^{n\times n}$ as $B=X_{o}A^{T}\Sigma AX_{o}$.
Then, be the Hanson-Wright inequality (Theorem VIII.3), we have

$\displaystyle\operatorname{P}(|{1\over\sigma_{w}^{2}}{\bf w}^{T}B{\bf w}-{\rm
Tr}(\Sigma AX_{o}^{2}A^{T})|>t)\leq 2\exp\Big{(}-c\min({t^{2}\over 4\|B\|_{\rm
HS}^{2}},{t\over 2\|B\|^{2}})\Big{)}.$ (26)

But

$\displaystyle\|B\|_{\rm HS}^{2}={\rm
Tr}(B^{2})=\sum_{i=1}^{m}\lambda_{i}^{2}(B)\leq
m\lambda^{2}_{\max}(B)=m\|B\|^{2}.$ (27)

On the other hand, $\|B\|=\|X_{o}A^{T}\Sigma AX_{o}\|\leq
x^{2}_{\max}\|A^{T}\Sigma A\|$. Moreover,

$\displaystyle\|A^{T}\Sigma A\|^{2}=\max_{{\bf u}\in\mathds{R}^{n}}{{\bf
u}^{T}A^{T}\Sigma AA^{T}\Sigma A{\bf u}\over\|{\bf
u}\|^{2}}\leq\lambda_{\max}(A^{T}A)\lambda_{\max}(AA^{T})\lambda_{\max}^{2}(\Sigma)$
(28)

But $\Sigma=(AX^{2}A^{T})^{-1}$ and $X={\rm diag}({\bf x})$. Therefore,
$\lambda_{\max}(\Sigma)=(\lambda_{\min}(AX^{2}A^{T}))^{-1}\leq(\lambda_{\min}(AA^{T})x^{2}_{\min})^{-1}$
and

$\displaystyle\|A^{T}\Sigma
A\|^{2}\leq{\lambda_{\max}(AA^{T})\lambda_{\max}(A^{T}A)\over\lambda^{2}_{\min}(AA^{T})x^{4}_{\min}}.$
(29)

∎

######  Lemma VIII.6.

Consider $m\times m$ matrix defined as $\Delta\Sigma=\Sigma-\Sigma_{o}$, where
$\Sigma=(AX^{2}A^{T})^{-1}$ and $\Sigma=(AX_{o}^{2}A^{T})^{-1}$. Then,

${x^{4}_{\min}\lambda^{3}_{\min}(AA^{T})\over
x^{8}_{\max}\lambda^{4}_{\max}(AA^{T})}\|A(X^{2}-X_{o}^{2})\|^{2}_{\rm
HS}\leq{\rm
Tr}(\Sigma_{o}^{-1}\Delta\Sigma\Sigma_{o}^{-1}\Delta\Sigma)\leq{x^{4}_{\max}\lambda^{3}_{\max}(AA^{T})\over
x_{\min}^{8}\lambda^{4}_{\min}(AA^{T})}\|A(X^{2}-X_{o}^{2})\|^{2}_{\rm HS}$

###### Proof of Lemma VIII.6.

Using the vectorizing equality,

$\displaystyle{\rm
Tr}(\Sigma_{o}^{-1}\Delta\Sigma\Sigma_{o}^{-1}\Delta\Sigma)$
$\displaystyle\stackrel{{\scriptstyle\rm(a)}}{{=}}({\rm
Vec}(\Delta\Sigma))^{T}(\Sigma_{o}^{-1}\otimes\Sigma_{o}^{-1}){\rm
Vec}(\Delta\Sigma)$
$\displaystyle\stackrel{{\scriptstyle\rm(b)}}{{\geq}}\|{\rm
Vec}(\Delta\Sigma)\|_{2}^{2}\lambda_{\min}(\Sigma_{o}^{-1}\otimes\Sigma_{o}^{-1}))$
$\displaystyle\stackrel{{\scriptstyle\rm(c)}}{{=}}\|\Delta\Sigma\|^{2}_{\rm
HS}\lambda^{2}_{\min}(\Sigma_{o}^{-1})$
$\displaystyle=\|\Delta\Sigma\|^{2}_{\rm
HS}\lambda^{2}_{\min}(AX_{o}^{2}A^{T}),$ (30)

where $\otimes$ denotes the Kronecker product operation. Steps (a), (b) and
(c) follow from the results on vectorization and Kronecker product. (Refer for
instance to Chapter 13 of [29].) But

$\displaystyle\lambda_{\min}(AX_{o}^{2}A^{T})$ $\displaystyle=\min_{{\bf
u}\neq 0}{\|X_{o}A^{T}{\bf u}\|^{2}\over\|A^{T}{\bf u}\|^{2}}{\|A^{T}{\bf
u}\|^{2}\over\|{\bf u}\|^{2}}\geq x^{2}_{\min}\lambda_{\min}(AA^{T}),$ (31)

which yields the desired result.

On the other hand,
$\Sigma_{o}-\Sigma=\Sigma(\Sigma^{-1}-\Sigma_{o}^{-1})\Sigma_{o}=\Sigma
A(X^{2}-X_{o}^{2})A^{T}\Sigma_{o}$. Moreover, $\|BC\|_{\rm
HS}\geq\sigma_{\min}(B)\|\|C\|_{\rm HS}$ and $\|BC\|_{\rm
HS}\geq\sigma_{\min}(C)\|B\|_{\rm HS}$. Therefore,

$\displaystyle\|\Sigma_{o}-\Sigma\|_{\rm HS}$
$\displaystyle\geq\lambda_{\min}(\Sigma)\sigma_{\min}(A^{T}\Sigma_{o})\|A(X^{2}-X_{o}^{2})\|_{\rm
HS}$ $\displaystyle\geq{\lambda^{1\over
2}_{\min}(AA^{T})\over\lambda_{\max}(AX^{2}A^{T})\lambda_{\max}(AX^{2}_{o}A^{T})}\|A(X^{2}-X_{o}^{2})\|_{\rm
HS}$ $\displaystyle\geq{\lambda^{1\over 2}_{\min}(AA^{T})\over
x_{\max}^{4}\lambda^{2}_{\max}(AA^{T})}\|A(X^{2}-X_{o}^{2})\|_{\rm HS}.$ (32)

Combining (30), (31) and (34) yields the desired lower bound. Similarly, to
obtain an upper bound for ${\rm
Tr}(\Sigma_{o}^{-1}\Delta\Sigma\Sigma_{o}^{-1}\Delta\Sigma)$, note that

$\displaystyle{\rm
Tr}(\Sigma_{o}^{-1}\Delta\Sigma\Sigma_{o}^{-1}\Delta\Sigma)$
$\displaystyle{=}({\rm
Vec}(\Delta\Sigma))^{T}(\Sigma_{o}^{-1}\otimes\Sigma_{o}^{-1}){\rm
Vec}(\Delta\Sigma)$ $\displaystyle{\leq}\|{\rm
Vec}(\Delta\Sigma)\|_{2}^{2}\lambda_{\max}(\Sigma_{o}^{-1}\otimes\Sigma_{o}^{-1}))$
$\displaystyle=\|\Delta\Sigma\|^{2}_{\rm
HS}\lambda^{2}_{\max}(\Sigma_{o}^{-1})$
$\displaystyle=\|\Delta\Sigma\|^{2}_{\rm
HS}\lambda^{2}_{\max}(AX_{o}^{2}A^{T})$
$\displaystyle\leq\|\Delta\Sigma\|^{2}_{\rm
HS}x^{2}_{\max}\lambda_{\max}(AA^{T}).$ (33)

Furthermore, using similar techniques as those used in deriving (34), we have

$\displaystyle\|\Sigma_{o}-\Sigma\|_{\rm HS}$
$\displaystyle\leq\lambda_{\max}(\Sigma)\sigma_{\max}(A^{T}\Sigma_{o})\|A(X^{2}-X_{o}^{2})\|_{\rm
HS}$ $\displaystyle\leq{\lambda^{1\over
2}_{\max}(AA^{T})\over\lambda_{\min}(AX^{2}A^{T})\lambda_{\min}(AX^{2}_{o}A^{T})}\|A(X^{2}-X_{o}^{2})\|_{\rm
HS}$ $\displaystyle\leq{\lambda^{1\over 2}_{\max}(AA^{T})\over
x_{\min}^{4}\lambda^{2}_{\min}(AA^{T})}\|A(X^{2}-X_{o}^{2})\|_{\rm HS}.$ (34)

∎

######  Lemma VIII.7.

Let the elements of $m\times n$ matrix $A$ be drawn i.i.d. $\mathcal{N}(0,1)$.
For any given $\bf{d}\in\mathbb{R}^{n}$, define $D=\mathop{\rm
diag}\nolimits(\bf{d})$. Then, for any $t\in(0,1)$,

$\operatorname{P}(\|AD\|^{2}_{HS}\leq m(1-t)\|{\bf
d}\|_{2}^{2},\;\forall\;{\bf d})\leq n\exp(\frac{m}{2}(t+\log(1-t))),$

and, for any $t>0$,

$\operatorname{P}(\|AD\|^{2}_{HS}\geq m(1+t)\|{\bf
d}\|_{2}^{2},\;\forall\;{\bf d})\leq n\exp(-\frac{m}{2}(t-\log(1+t))).$

###### Proof.

Let $A=[{\bf a}_{1};\ldots,{\bf a}_{n}]$ with ${\bf a}_{i}\in\mathds{R}^{m}$.
Then, by definition,

$\|AD\|^{2}_{HS}=\sum_{i=1}^{n}d_{i}^{2}\|{\bf a}_{i}\|_{2}^{2}.$ (35)

Define event $\mathcal{E}$ as the event where $\min_{i}\|{\bf
a}_{i}\|_{2}^{2}>m(1-t)$. Using Lemma VIII.2 and applying the union bound, it
follows that $\operatorname{P}(\mathcal{E}^{c})\geq n{\rm
e}^{\frac{m}{2}(t+\log(1-t))}$. It is clear that under $\mathcal{E}$, we have

$\|AD\|^{2}_{HS}=\sum_{i=1}^{n}d_{i}^{2}\|A_{i}\|_{2}^{2}\geq m(1-t)\|{\bf
d}\|_{2}^{2}.$

∎

###  VIII-B Proof of Theorem IV.1

Using the matrix inversion lemma, we have

$\displaystyle(I_{n}+{\sigma_{w}^{2}\over\sigma_{z}^{2}}XA^{T}AX)^{-1}=I_{n}-{\sigma_{w}^{2}\over\sigma_{z}^{2}}XA^{T}\Big{(}I_{n}+{\sigma_{w}^{2}\over\sigma_{z}^{2}}AX^{2}A^{T}\Big{)}^{-1}AX.$
(36)

Inserting this in (5), it follows that

$\displaystyle 2\ell(X)$
$\displaystyle=-\log\det({1\over\sigma_{w}^{2}}I_{n}+{1\over\sigma_{z}^{2}}XA^{T}AX)+{\sigma_{w}^{2}\over\sigma_{z}^{4}}{\bf
y}^{T}AX^{2}A^{T}{\bf y}$
$\displaystyle\;\;\;\;-{\sigma_{w}^{4}\over\sigma_{z}^{6}}{\bf
y}^{T}AX^{2}A^{T}\Big{(}I_{n}+{\sigma_{w}^{2}\over\sigma_{z}^{2}}AX^{2}A^{T}\Big{)}^{-1}AX^{2}A^{T}{\bf
y}.$ (37)

But,

$\displaystyle{\sigma_{w}^{4}\over\sigma_{z}^{6}}{\bf
y}^{T}AX^{2}A^{T}\Big{(}I_{n}+{\sigma_{w}^{2}\over\sigma_{z}^{2}}AX^{2}A^{T}\Big{)}^{-1}AX^{2}A^{T}{\bf
y}$ $\displaystyle={\sigma_{w}^{2}\over\sigma_{z}^{4}}{\bf
y}^{T}(I_{n}+{\sigma_{w}^{2}\over\sigma_{z}^{2}}AX^{2}A^{T}-I_{n})\Big{(}I_{n}+{\sigma_{w}^{2}\over\sigma_{z}^{2}}AX^{2}A^{T}\Big{)}^{-1}AX^{2}A^{T}{\bf
y}$ $\displaystyle={\sigma_{w}^{2}\over\sigma_{z}^{4}}{\bf
y}^{T}AX^{2}A^{T}{\bf y}-{\sigma_{w}^{2}\over\sigma_{z}^{4}}{\bf
y}^{T}\Big{(}I_{n}+{\sigma_{w}^{2}\over\sigma_{z}^{2}}AX^{2}A^{T}\Big{)}^{-1}AX^{2}A^{T}{\bf
y}.$ (38)

Therefore, cancelling the common terms, we have

$\displaystyle 2\ell(X)$
$\displaystyle=-\log\det({1\over\sigma_{w}^{2}}I_{n}+{1\over\sigma_{z}^{2}}XA^{T}AX)+{\sigma_{w}^{2}\over\sigma_{z}^{4}}{\bf
y}^{T}\Big{(}I_{n}+{\sigma_{w}^{2}\over\sigma_{z}^{2}}AX^{2}A^{T}\Big{)}^{-1}AX^{2}A^{T}{\bf
y}$
$\displaystyle=-\log\det({1\over\sigma_{w}^{2}}I_{n}+{1\over\sigma_{z}^{2}}XA^{T}AX)+{1\over\sigma_{z}^{2}}{\bf
y}^{T}\Big{(}I_{n}+{\sigma_{w}^{2}\over\sigma_{z}^{2}}AX^{2}A^{T}\Big{)}^{-1}\Big{(}{\sigma_{w}^{2}\over\sigma_{z}^{2}}AX^{2}A^{T}+I_{n}-I_{n}\Big{)}{\bf
y}$
$\displaystyle=-\log\det({1\over\sigma_{w}^{2}}I_{n}+{1\over\sigma_{z}^{2}}XA^{T}AX)+{{\bf
y}^{T}{\bf y}\over\sigma_{z}^{2}}-{1\over\sigma_{z}^{2}}{\bf
y}^{T}\Big{(}I_{n}+{\sigma_{w}^{2}\over\sigma_{z}^{2}}AX^{2}A^{T}\Big{)}^{-1}{\bf
y}.$ (39)

Assuming that $\sigma_{z}$ is converging to zero, ignoring the terms not
depending on $X$, we have

$\displaystyle-2\ell(X)$
$\displaystyle=\log\det({1\over\sigma_{w}^{2}}I_{n}+{1\over\sigma_{z}^{2}}XA^{T}AX)+{1\over\sigma_{w}^{2}}{\bf
y}^{T}(AX^{2}A^{T})^{-1}{\bf y}.$ (40)

Let $\lambda_{1},\ldots,\lambda_{m}$ denote the non-zero eigenvalues of
$XA^{T}AX$. Then, ${1\over\sigma_{w}^{2}}I_{n}+{1\over\sigma_{z}^{2}}XA^{T}AX$
has $n-m$ eigenvalues equal to ${1\over\sigma_{w}^{2}}$ and the rest of its
eigenvalues are ${1\over\sigma_{w}^{2}}+{1\over\sigma_{z}^{2}}\lambda_{i}$,
$i=1,\ldots,m$. Therefore,

$\displaystyle\log\det({1\over\sigma_{w}^{2}}I_{n}+{1\over\sigma_{z}^{2}}XA^{T}AX)$
$\displaystyle=-2(n-m)\log\sigma_{w}+\sum_{i=1}^{m}\log({1\over\sigma_{w}^{2}}+{1\over\sigma_{z}^{2}}\lambda_{i})$
$\displaystyle=-2(n-m)\log\sigma_{w}-2m\log\sigma_{z}+\sum\log\lambda_{i}+O(\sigma_{z}^{2})$
$\displaystyle=\sum_{i=1}^{m}\log\lambda_{i}-2(n-m)\log\sigma_{w}-2m\log\sigma_{z}+O(\sigma_{z}^{2}).$
(41)

Note that if $XA^{T}AX{\bf u}=\lambda_{i}{\bf u}$, then $AX^{2}A^{T}(AX{\bf
u})=\lambda_{i}(AX{\bf u})$. Therefore, $\lambda_{i}$, $i=1,\ldots,m$, are the
eigenvalues of $AX^{2}A^{T}$ as well. Hence, in summary, again by ignoring the
terms that do not depend on $X$, by a slight abuse of notation, as
$\sigma_{z}\to 0$, we have

$\displaystyle-2\ell(X)$
$\displaystyle=\log\det(AX^{2}A^{T})+{1\over\sigma_{w}^{2}}{\bf
y}^{T}(AX^{2}A^{T})^{-1}{\bf y}.$ (42)

###  VIII-C Proof of Theorem V.1

Since ${\bf y}=AX_{o}{\bf w}$, function $f$ can be written as
$f(\Sigma)=-\log\det\Sigma+{1\over\sigma_{w}^{2}}{\rm Tr}(\Sigma AX_{o}{\bf
w}{\bf w}^{T}X_{o}A^{T})$. Taking the expected value of $f$ with respect to
the noise ${\bf w}$, we derive

$\displaystyle\bar{f}(\Sigma)=-\log\det\Sigma+{\rm Tr}(\Sigma
AX_{o}^{2}A^{T}).$ (43)

As a function of $\Sigma$, $\bar{f}$ is a convex function that achieves its
minimum at $\Sigma^{-1}=AX_{o}^{2}A^{T}$ or at $X$ satisfying
$AX^{2}A^{T}=AX_{o}^{2}A^{T}$. Define

$\Sigma_{o}\triangleq(AX_{o}^{2}A^{T})^{-1}.$

Given $\hat{{\bf x}}_{o}$ (the minimizer of $f$), let
$\hat{\Sigma}_{o}=\Sigma(\hat{{\bf x}}_{o})$. Moreover, define $\tilde{{\bf
x}}_{o}$ as the closest reconstruction signal in $\mathcal{C}$ to ${\bf
x}_{o}$, i.e.,

$\tilde{{\bf x}}_{o}=\operatorname*{arg\,min}_{{\bf x}\in\mathcal{C}}\|{\bf
x}_{o}-{\bf x}\|.$

Let $\tilde{\Sigma}_{o}\triangleq(A{\tilde{X}}_{o}^{2}A^{T})^{-1}$, where
${\tilde{X}}_{o}={\rm diag}({\tilde{{\bf x}}_{o}})$. Since $\hat{{\bf x}}_{o}$
is the minimizer of (12),

$\displaystyle f(\hat{\Sigma}_{o})\leq f(\tilde{\Sigma}_{o}).$ (44)

Define $\Delta\Sigma$ as $\hat{\Sigma}_{o}-\Sigma_{o}$ and let $\lambda_{m}$
denote the maximum eigenvalue of $\Sigma_{o}^{-{1\over
2}}\Delta\Sigma\Sigma_{o}^{-{1\over 2}}$. Then, as shown in Lemma VIII.4,
$\bar{f}(\hat{\Sigma}_{o})-\bar{f}(\Sigma_{o})$ can be lower bounded as
follows:

$\displaystyle\bar{f}(\hat{\Sigma}_{o})-\bar{f}(\Sigma_{o})$
$\displaystyle\geq{1\over 2(1+\lambda_{m})^{2}}{\rm
Tr}(\Sigma_{o}^{-1}\Delta\Sigma\Sigma_{o}^{-1}\Delta\Sigma),$ (45)

where $\Delta\Sigma=\hat{\Sigma}_{o}-\Sigma_{o}$. Let $\delta
f(\hat{\Sigma}_{o})\triangleq{f}(\hat{\Sigma}_{o})-\bar{f}(\hat{\Sigma}_{o})$
and $\delta
f(\tilde{\Sigma}_{o})\triangleq{f}(\tilde{\Sigma}_{o})-\bar{f}(\tilde{\Sigma}_{o})$.
Then,

$\displaystyle\bar{f}(\hat{\Sigma}_{o})-\bar{f}(\Sigma_{o})$
$\displaystyle={f}(\hat{\Sigma}_{o})-\delta
f(\hat{\Sigma}_{o})-(\bar{f}(\Sigma_{o})-{f}(\tilde{\Sigma}_{o})+{f}(\tilde{\Sigma}_{o})-\bar{f}(\tilde{\Sigma}_{o})+\bar{f}(\tilde{\Sigma}_{o}))$
$\displaystyle={f}(\hat{\Sigma}_{o})-{f}(\tilde{\Sigma}_{o})-\delta
f(\hat{\Sigma}_{o})+\delta
f(\tilde{\Sigma}_{o})-\bar{f}(\Sigma_{o})+\bar{f}(\tilde{\Sigma}_{o}).$ (46)

Therefore, combining (44) and (45), it follows that

${1\over 2(1+\lambda_{m})^{2}}{\rm
Tr}(\Sigma_{o}^{-1}\Delta\Sigma\Sigma_{o}^{-1}\Delta\Sigma)\leq|\delta
f(\tilde{\Sigma}_{o})|+|\delta
f(\hat{\Sigma}_{o})|+|\bar{f}(\Sigma_{o})-\bar{f}(\tilde{\Sigma}_{o})|.$

Also, applying Lemma VIII.6, it follows that

${1\over 2(1+\lambda_{m})^{2}}{x^{4}_{\min}\lambda^{3}_{\min}(AA^{T})\over
x^{8}_{\max}\lambda^{4}_{\max}(AA^{T})}\|A(\hat{X}_{o}^{2}-X_{o}^{2})\|^{2}_{\rm
HS}\leq|\delta f(\tilde{\Sigma}_{o})|+|\delta
f(\hat{\Sigma}_{o})|+|\bar{f}(\Sigma_{o})-\bar{f}(\tilde{\Sigma}_{o})|.$

Given $t_{1}\in(0,1)$, $t_{2}>0$ and $t_{3}>0$, define events
$\mathcal{E}_{1}$, $\mathcal{E}_{2}$, $\mathcal{E}_{3}$ and $\mathcal{E}_{4}$
as

$\mathcal{E}_{1}=\\{\|A(\hat{X}_{o}^{2}-X_{o}^{2})\|^{2}_{\rm HS}\geq
m(1-t_{1})\|{\bf x}^{2}_{o}-\hat{{\bf x}}^{2}_{o}\|^{2}\\},$
$\mathcal{E}_{2}=\\{|\delta f(\tilde{\Sigma}_{o})|\leq
t_{2}\\},\;\;\;\mathcal{E}_{3}=\\{|\delta f(\hat{\Sigma}_{o})|\leq t_{2}\\},$
$\mathcal{E}_{4}=\\{\sqrt{n}-2\sqrt{m}\leq\sigma_{\min}(A)\leq\sigma_{\max}(A)\leq\sqrt{n}+2\sqrt{m}\\},$

and

$\mathcal{E}_{5}=\\{\|A(\tilde{X}_{o}^{2}-X^{2}_{o})\|^{2}_{\rm HS}\leq
m(1+t_{3})\|{\bf x}_{o}^{2}-{\tilde{{\bf x}}^{2}_{o}}\|^{2}\\},$

respectively. Conditioned on
$\mathcal{E}_{1}\cap\mathcal{E}_{2}\cap\mathcal{E}_{3}\cap\mathcal{E}_{4}$, we
have

${m\over 2(1+\lambda_{m})^{2}}{x^{4}_{\min}(\sqrt{n}-2\sqrt{m})^{6}\over
x^{8}_{\max}(\sqrt{n}+2\sqrt{m})^{8}}(1-t_{1})\|{\bf x}^{2}-\hat{{\bf
x}}_{o}^{2}\|_{2}^{2}\leq
2t_{2}+|\bar{f}(\Sigma_{o})-\bar{f}(\tilde{\Sigma}_{o})|.$ (47)

But, from Lemma VIII.4, we have

$\displaystyle|\bar{f}(\Sigma_{o})-\bar{f}(\tilde{\Sigma}_{o})|$
$\displaystyle\leq{1\over 2}{\rm
Tr}(\Sigma_{o}^{-1}(\tilde{\Sigma}_{o}-\Sigma_{o})\Sigma_{o}^{-1}(\tilde{\Sigma}_{o}-\Sigma_{o}))$
$\displaystyle\leq{x^{4}_{\max}\lambda^{3}_{\max}(AA^{T})\over
2{x_{\min}^{8}}\lambda^{4}_{\min}(AA^{T})}\|A(\tilde{X}^{2}_{o}-X^{2}_{o})\|^{2}_{\rm
HS},$ (48)

where the second inequality follows from Lemma VIII.6. Then, conditioned on
$\mathcal{E}_{4}\cap\mathcal{E}_{5}$, it follows that

$\displaystyle|\bar{f}(\Sigma_{o})-\bar{f}(\tilde{\Sigma}_{o})|$
$\displaystyle\leq{x^{4}_{\max}(\sqrt{n}+2\sqrt{m})^{6}\over
2{x_{\min}^{8}}(\sqrt{n}-2\sqrt{m})^{8}}m(1+t_{3})\|{\bf
x}_{o}^{2}-{\tilde{{\bf x}}^{2}_{o}}\|^{2}.$ (49)

Define

$c_{1}={x^{4}_{\min}(1-2\sqrt{m/n})^{6}\over
x^{8}_{\max}(1+2\sqrt{m/n})^{8}},$

and

$c_{2}={x^{4}_{\max}(1+2\sqrt{m/n})^{6}\over{x_{\min}^{8}}(1-2\sqrt{m/n})^{8}}.$

Using these definitions and combining (48) and (49), it follows that

$\displaystyle{m\over 2(1+\lambda_{m})^{2}}{c_{1}\over n}(1-t_{1})\|{\bf
x}_{o}^{2}-\hat{{\bf x}}^{2}_{o}\|_{2}^{2}\leq 2t_{2}+{c_{2}\over
2n}m(1+t_{3})\|{\bf x}^{2}_{o}-{\tilde{{\bf x}}^{2}_{o}}\|^{2}.$ (50)

Recall that $\lambda_{m}$ is defined as the maximum eigenvalue of
$\Sigma_{o}^{-{1\over 2}}\Delta\Sigma\Sigma_{o}^{-{1\over 2}}$. On the other
hand, $\lambda_{m}=\|\Sigma_{o}^{-{1\over 2}}\Delta\Sigma\Sigma_{o}^{-{1\over
2}}\|\leq\|\Delta\Sigma\|\|\Sigma_{o}^{-{1\over
2}}\|^{2}=\|\Delta\Sigma\|\|\Sigma_{o}^{-1}\|$. But,
$\|\Sigma_{o}^{-1}\|=\|AX_{o}^{2}A^{T}\|\leq
x^{2}_{\max}\lambda_{\max}(AA^{T})$. Similarly,
$\|\Delta\Sigma\|=\|\hat{\Sigma}_{o}-{\Sigma}_{o}\|\leq\|\hat{\Sigma}_{o}\|-\|{\Sigma}_{o}\|\leq{1\over\|AX_{o}^{2}A^{T}\|}+{1\over\|A\hat{X}_{o}^{2}A^{T}\|}\leq{2\over
x^{2}_{\min}\lambda_{\min}(AA^{T})}$. So overall,
$\lambda_{m}\leq{2x^{2}_{\max}\lambda_{\max}(AA^{T})\over
x^{2}_{\min}\lambda_{\min}(AA^{T})}$, and conditioned on $\mathcal{E}_{4}$, we
have

$\displaystyle\lambda_{m}\leq{2x^{2}_{\max}(1+2\sqrt{m/n})^{2}\over
x^{2}_{\min}(1-2\sqrt{m/n})^{2}}.$ (51)

To finish the proof we bound the probability of
$\cap_{i=1}^{5}\mathcal{E}_{i}$. Lemma VIII.7 can be used to bound
$\operatorname{P}(\mathcal{E}_{1}^{c})+\operatorname{P}(\mathcal{E}_{5}^{c})$
as

$\operatorname{P}(\mathcal{E}_{1}^{c})\leq
n\exp(\frac{m}{2}(t_{1}+\log(1-t_{1}))),$

and

$\operatorname{P}(\mathcal{E}_{5}^{c})\leq
n\exp(-\frac{m}{2}(t_{3}-\log(1+t_{3}))).$

Setting $t_{1}=0.5$ and $t_{3}=1$, it follows that
$\operatorname{P}(\mathcal{E}_{1}^{c})\leq n{\rm e}^{-0.09m}$ and
$\operatorname{P}(\mathcal{E}_{5}^{c})\leq n{\rm e}^{-0.84m}$. From Lemma
VIII.1, $\operatorname{P}(\mathcal{E}_{4}^{c})\leq 2\exp(-\frac{m}{2})$.
Setting $t_{3}=1$, Lemma VIII.7 implies that Finally, from Lemma VIII.5
combined with the union bound (since $|\mathcal{C}|\leq 2^{nr}$) implies that

$\operatorname{P}((\mathcal{E}_{2}\cap\mathcal{E}_{3})^{c})\leq
2^{nr+1}\exp\Big{(}-{ct_{2}\min(1,{t_{2}\over
2m})\lambda^{2}_{\min}(AA^{T})x^{4}_{\min}\over
2x_{\max}^{4}\lambda^{2}_{\max}(A^{T}A)}\Big{)}.$

Let

$c_{3}\triangleq{cx^{4}_{\min}(1-2\sqrt{m/n})^{4}\over
2x_{\max}^{4}(1+2\sqrt{m/n})^{4}}.$

Then, for $t_{2}<2m$,

$\displaystyle\operatorname{P}((\mathcal{E}_{2}\cap\mathcal{E}_{3})^{c}\cup\mathcal{E}_{4})$
$\displaystyle\leq 2^{nr+1}\exp(-{c_{3}t_{2}^{2}/2m})$
$\displaystyle=\exp((nr+1)\log 2-{c_{3}t_{2}^{2}/2m}).$ (52)

Let $t_{2}=\sqrt{(2\log 2)mnr(1+\epsilon)/c_{3}}$. Then,

$\displaystyle\operatorname{P}((\mathcal{E}_{2}\cap\mathcal{E}_{3})^{c}\cup\mathcal{E}_{4})$
$\displaystyle\leq\exp((nr+1)\log 2-(\log 2)nr(1+\epsilon))$
$\displaystyle=2^{-nr\epsilon+1}.$ (53)

Let $c_{4}=c_{3}/(8\log 2)$. Since by assumption ${1\over n}\|{\bf
x}_{o}-{\tilde{{\bf x}}_{o}}\|^{2}\leq\delta$, we have ${1\over n}\|{\bf
x}_{o}^{2}-\tilde{{\bf x}}_{o}^{2}\|^{2}\leq 4x_{\max}^{2}\delta$. Also,
$\|\hat{{\bf x}}_{o}^{2}-{\bf x}_{o}^{2}\|^{2}\geq 4x_{\min}^{2}\|\hat{{\bf
x}}_{o}-{\bf x}_{o}\|^{2}$. Therefore, from (50) it follows that

$\displaystyle{c_{1}x^{2}_{\min}\over(1+\lambda_{m})^{2}}{1\over n}\|{\bf
x}_{o}-\hat{{\bf x}}_{o}\|_{2}^{2}\leq\sqrt{(1+\epsilon)nr\over
c_{4}m}+4c_{2}x^{2}_{\max}\delta.$ (54)

##  IX Conclusions

In this paper, we have studied compressed sensing recovery of structured
signals in the presence of speckle noise. In a compressed sensing system where
every input pixel is distorted independently by a multiplicative Gaussian
noise, we have derived a ML-based recovery method. We have used lossy
compression codes to model the structures of sources. We have shown that given
sufficient number of measurements the ML-based recovery method is able to
recover a signal from its under-sampled measurements, even in the presence of
speckle noise. To the best of our knowledge, this is the first theoretical
result on estimation in the presence of speckle noise. The ML-based
optimization is computationally intractable and cannot be implemented. We have
proposed employing projected gradient descent to approximate its solution. Our
simulations results show the effectiveness of the proposed method.

First note that $W=\sum_{i=1}^{n}w_{i}^{2}$ has a $\chi^{2}$ distribution $n$
degrees of freedom. Hence, its distribution is given by

$f_{W}(w)=\frac{w^{\frac{n}{2}-1}\exp(-\frac{n}{2})}{2^{\frac{n}{2}}\Gamma(\frac{n}{2})}.$

Hence,
$\operatorname{E}(W^{\frac{1}{2}})=\frac{\Gamma(\frac{n+1}{2})\sqrt{2}}{\Gamma(\frac{n}{2})}$.
Using the Stirling’s formula for the $\Gamma$ function, we have

$\displaystyle\log(\Gamma(\frac{n+1}{2}))-\log(\Gamma(\frac{n}{2}))$
$\displaystyle=$
$\displaystyle\frac{n+1}{2}\log(\frac{n+1}{2})-\frac{n+1}{2}+\frac{1}{2}\log(\frac{4\pi}{n+1})-\frac{n}{2}\log(\frac{n}{2})+\frac{n}{2}-\frac{1}{2}\log\frac{4\pi}{n}+O(\frac{1}{n^{2}})$
$\displaystyle=$
$\displaystyle\frac{1}{2}\log(\frac{n+1}{2})+\frac{n}{2}\log(\frac{n+1}{n})-\frac{1}{2}-\frac{1}{2}\log(\frac{n+1}{n})+O(\frac{1}{n^{2}})$
$\displaystyle=$
$\displaystyle\frac{1}{2}\log\frac{n}{2}+\frac{1}{2n}+\frac{n}{2}(\frac{1}{n}-\frac{1}{2n^{2}})-\frac{1}{2}-\frac{1}{2n}+O(\frac{1}{n^{2}}).$
$\displaystyle=$
$\displaystyle\frac{1}{2}\log\frac{n}{2}-\frac{1}{4n}+O(\frac{1}{n^{2}}).$

Hence,

$\displaystyle\operatorname{E}(W^{\frac{1}{2}})=\sqrt{n}(1-\frac{1}{4n})+o(\frac{1}{\sqrt{n}}).$

Therefore,

$\displaystyle\operatorname{E}\|\hat{{\bf x}}_{\rm ML}-{\bf
x}\|_{2}^{2}=a^{2}n\operatorname{E}(({1\over
n}\sum_{i=1}^{n}w_{i}^{2})^{1\over
2})-1)^{2}=a^{2}n(2-2\operatorname{E}({1\over
n}\sum_{i=1}^{n}w_{i}^{2})^{1\over 2})=\frac{a^{2}}{2}.$

## References

  * [1] A. Moreira, P. Prats-Iraola, M. Younis, G. Krieger, I. Hajnsek, and K. P. Papathanassiou.  A tutorial on synthetic aperture radar.  IEEE Geo. and Rem. Sen. Mag., 1(1):6–43, 2013.
  * [2] D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, et al.  Optical coherence tomography.  Science, 254(5035):1178–1181, 1991.
  * [3] F. Argenti, A. Lapini, T. Bianchi, and L. Alparone.  A tutorial on speckle reduction in synthetic aperture radar images.  IEEE Geo. and Rem. Sen. Mag., 1(3):6–35, Sep. 2013.
  * [4] R. Touzi.  A review of speckle filtering in the context of estimation theory.  IEEE Trans. on Geo. and Rem. Sen., 40(11):2392–2404, 2002.
  * [5] C. Ozcan, B. Sen, and F. Nar.  Sparsity-driven despeckling for SAR images.  IEEE Geo. and Rem. Sen. Letters, 13(1):115–119, 2016.
  * [6] C. Deledalle, L. Denis, G. Poggi, F. Tupin, and L. Verdoliva.  Exploiting patch similarity for SARs image processing: The nonlocal paradigm.  IEEE Sig. Proc. Mag., 31(4):69–78, 2014.
  * [7] G. Di Martino, A. Di Simone, A. Iodice, and D. Riccio.  Scattering-based nonlocal means SAR despeckling.  IEEE Trans. on Geo. and Rem. Sen., 54(6):3574–3588, 2016.
  * [8] G. Chierchia, D. Cozzolino, G. Poggi, and L. Verdoliva.  SAR image despeckling through convolutional neural networks.  In 2017 IEEE Int. Geo. and Rem. Sen. Symp. (IGARSS), pages 5438–5441, 2017.
  * [9] P. Wang, H. Zhang, and V. M. Patel.  SAR image despeckling using a convolutional neural network.  IEEE Sig. Proc. Letters, 24(12):1763–1767, 2017.
  * [10] Yeo-Sun Yoon and Moeness G. Amin.  Compressed sensing technique for high-resolution radar imaging.  In Ivan Kadar, editor, Signal Processing, Sensor Fusion, and Target Recognition XVII, volume 6968, pages 506 – 515. International Society for Optics and Photonics, SPIE, 2008.
  * [11] V. M. Patel, G. R. Easley, D. M. Healy, and R. Chellappa.  Compressed sensing for synthetic aperture radar imaging.  In IEEE Int. Conf. on Image Proc. (ICIP), pages 2141–2144, 2009\.
  * [12] N. O. Onhon and M. Çetin.  Joint sparsity-driven inversion and model error correction for radar imaging.  In 2010 IEEE Int. Conf. on Acou. Speech and Sig. Pro., pages 1206–1209, 2010.
  * [13] S. Demirci and C. Ozdemir.  Compressed sensing-based imaging of millimeter-wave ISAR data.  Mic. and Opt. Tech. Letters, 55(12):2967–2972, 2013.
  * [14] P. Cheng and J. Zhao.  Generalised Pareto distribution-based Bayesian compressed sensing inverse synthetic aperture radar imaging.  IET Radar, Son. & Nav., 12(5):549–556, 2018.
  * [15] D. Bi, Y. Xie, L. Ma, X. Li, X. Yang, and Y. R. Zheng.  Multifrequency compressed sensing for 2-d near-field synthetic aperture radar image reconstruction.  IEEE Trans. on Inst. and Mea., 66(4):777–791, 2017.
  * [16] M. Cetin, I. Stojanović, N. O. Onhon, K. Varshney, S. Samadi, W. C. Karl, and A. S. Willsky.  Sparsity-driven synthetic aperture radar imaging: Reconstruction, autofocusing, moving targets, and compressed sensing.  IEEE Sig. Proc. Mag., 31(4):27–40, 2014.
  * [17] S. Jalali and A. Maleki.  From compression to compressed sensing.  Appl. Comp. Harmonic Anal. (ACHA), 40(2):352–385, 2016.
  * [18] S. Beygi, S. Jalali, A. Maleki, and U. Mitra.  An efficient algorithm for compression-based compressed sensing.  Information and Inference: A Journal of the IMA, 8(2):343–375, 2019\.
  * [19] M. Bakhshizadeh, A. Maleki, and S. Jalali.  Using black-box compression algorithms for phase retrieval.  IEEE Trans. Inform. Theory, 66(12):7978–8001, 2020.
  * [20] F. E. Rezagah, S. Jalali, E. Erkip, and H. V. Poor.  Compression-based compressed sensing.  IEEE Trans. Inform. Theory, 63(10):6735–6752, Oct. 2017.
  * [21] P. J. Bickel, Y. Ritov, and A. B.s Tsybakov.  Simultaneous analysis of Lasso and Dantzig selector.  The Annals of statistics, 37(4):1705–1732, 2009.
  * [22] Patrick L Combettes and Valérie R Wajs.  Signal recovery by proximal forward-backward splitting.  Multiscale Modeling $\&$ Simulation, 4(4):1168–1200, 2005.
  * [23] T. Blumensath and M. E. Davies.  Iterative hard thresholding for compressed sensing.  Appl. Comp. Harmonic Anal. (ACHA), 27(3):265–274, 2009.
  * [24] R. Tibshirani, M. Saunders, S. Rosset, J. Zhu, and K. Knight.  Sparsity and smoothness via the fused Lasso.  J. of the Royal Stat. Soc.: Series B (Stat. Meth.), 67(1):91–108, 2005.
  * [25] A. Viterbi.  Error bounds for convolutional codes and an asymptotically optimum decoding algorithm.  IEEE Trans. Inform. Theory, 13(2):260 – 269, apr 1967.
  * [26] T. Akiba, S. Sano, T. Yanase, T. Ohta, and M. Koyama.  Optuna: A next-generation hyperparameter optimization framework.  In Proc. of the 25rd ACM SIGKDD Int. Conf. on Know. Dis. and Data Min., 2019.
  * [27] M. Rudelson and R. Vershynin.  Non-asymptotic theory of random matrices: extreme singular values.  In Proc. of the Int. Cong. of Math. 2010 (ICM 2010), pages 1576–1602. World Scientific, 2010.
  * [28] S. Jalali, A. Maleki, and R. G. Baraniuk.  Minimum complexity pursuit for universal compressed sensing.  IEEE Trans. Inform. Theory, 60(4):2253–2268, Apr. 2014.
  * [29] A. J. Laub.  Matrix analysis for scientists and engineers, volume 91.  Siam, 2005.