Simultaneous Sparse Recovery and Blind Demodulation

TL;DR

This paper introduces a novel method for simultaneous sparse signal recovery and blind demodulation, leveraging atomic norm minimization, with theoretical guarantees and practical validation in noisy and noiseless scenarios.

Contribution

It proposes a new approach using lifting and atomic norm minimization for joint sparse recovery and blind demodulation with theoretical sample complexity bounds.

Findings

Near optimal sample complexity bounds for perfect recovery.

Effective recovery in noisy conditions demonstrated.

Theoretical analysis supported by numerical simulations.

Abstract

The task of finding a sparse signal decomposition in an overcomplete dictionary is made more complicated when the signal undergoes an unknown modulation (or convolution in the complementary Fourier domain). Such simultaneous sparse recovery and blind demodulation problems appear in many applications including medical imaging, super resolution, self-calibration, etc. In this paper, we consider a more general sparse recovery and blind demodulation problem in which each atom comprising the signal undergoes a distinct modulation process. Under the assumption that the modulating waveforms live in a known common subspace, we employ the lifting technique and recast this problem as the recovery of a column-wise sparse matrix from structured linear measurements. In this framework, we accomplish sparse recovery and blind demodulation simultaneously by minimizing the induced atomic norm, which in this problem corresponds to the block $\ell_1$ norm minimization. For perfect recovery in the noiseless case, we derive near optimal sample complexity bounds for Gaussian and random Fourier overcomplete dictionaries. We also provide bounds on recovering the column-wise sparse matrix in the noisy case. Numerical simulations illustrate and support our theoretical results.

Equations

$$\mathbf{y}=\sum _{j=1}^M{c}_j{\mathbf{D}}_j{\mathbf{a}}_j\in {\mathbf{C}}^N$$

$${\mathbf{D}}_j=\mathrm{diag}(\mathbf{B}{\mathbf{h}}_j),$$

$$\mathbf{y}(n)$$

$$=⟨\sum _{j=1}^M{c}_j{\mathbf{h}}_j{\bar{\mathbf{a}}}_j^H,{\mathbf{b}}_n^{\prime }{\mathbf{e}}_n^H⟩=⟨\mathbf{G},{\mathbf{b}}_n^{\prime }{\mathbf{e}}_n^H⟩,$$

$$\underset{\mathbf{G}\in {\mathbf{C}}^{K\times N}}{\mathrm{minimize}}||\mathbf{G}|{|}_{\mathcal{A}}subject\; to\mathbf{y}={\mathcal{L}}^{\prime }(\mathbf{G}).$$

$$\underset{\mathbf{X}\in {\mathbf{C}}^{K\times M}}{\mathrm{minimize}}||\mathbf{X}|{|}_{2,1}subject\; to\mathbf{y}=\mathcal{L}(\mathbf{X})$$

$$\mathbf{y}(n)$$

$$||\mathbf{G}|{|}_{\mathcal{A}}$$

$$\displaystyle=\inf\Bigg{\{}\sum_{j=1}^{M}|c_{j}|:\mathbf{G}=\sum_{j=1}^{M}c_{j}\bm{h}_{j}\bar{\bm{a}}_{j}^{H},||\bm{h}_{j}||_{2}=1\Bigg{\}}$$

$$\displaystyle=\inf\Bigg{\{}\sum_{j=1}^{M}||\bm{x}_{j}||_{2}:\mathbf{G}=\sum_{j=1}^{M}\bm{x}_{j}\bar{\bm{a}}_{j}^{H}\Bigg{\}}$$

$$\displaystyle=\inf\big{\{}||\mathbf{X}||_{2,1}:\mathbf{G}=\mathbf{X}\bar{\mathbf{A}}^{H}\big{\}}.$$

$$\mathbf{y}(n)$$

$$c_j=||{\hat{\mathbf{x}}}_j|{|}_2,{\mathbf{h}}_j=\frac{{\hat{\mathbf{x}}}_j}{||{\hat{\mathbf{x}}}_j|{|}_2},\text{ and }{\mathbf{D}}_j=\text{diag}(\mathbf{B}{\mathbf{h}}_j)$$

$$\mathbf{\Phi }=[{\varphi }_{1,1}\cdots {\varphi }_{K,1}\cdots {\varphi }_{1,M}\cdots {\varphi }_{K,M}]$$

$${\mathbf{\Phi }}^H=[{\varphi }_1^{\prime }{\varphi }_2^{\prime }\cdots {\varphi }_N^{\prime }]\in {\mathbf{C}}^{KM\times N}$$

$$\mathbf{y}=\sum _{j=1}^M{c}_j{\mathbf{D}}_j{\mathbf{a}}_j+\mathbf{n}$$

$$\underset{\mathbf{X}\in {\mathbf{C}}^{K\times M}}{\mathrm{minimize}}||\mathbf{X}|{|}_{2,1}subject\; to||\mathbf{y}-\mathcal{L}(\mathbf{X})|{|}_2\le \eta .$$

$$\mathbf{y}=\mathbf{DA}(\theta )\mathbf{c}+\mathbf{n}\in {\mathbf{C}}^{N\times 1}$$

$$\mathbf{y}=\sum _{j=1}^M{c}_j{\mathbf{D}}_j\mathbf{a}({\theta }_j)+\mathbf{n}\in {\mathbf{C}}^{N\times 1}.$$

$$\mathbf{y}=Sample\left[\sum _{j=1}^M{c}_j({\mathbf{B}}^{\prime }{\mathbf{h}}_j)⊛{\mathbf{e}}_j+{\mathbf{n}}^{\prime }\right]\in {\mathbf{R}}^{N\times 1}$$

$$\displaystyle\bm{y}=Sample\Bigg{\{}IDFT\left[\sum_{j=1}^{M}c_{j}\mathbf{D}_{j}\bm{a}_{j}+\bm{n}\right]\Bigg{\}}\in\mathbf{R}^{N\times 1},$$

$$\frac{N}{{\log }^2(N)}\ge C_{\alpha }{\mu }_{max}^{2}KJ(\log (M-J)+\log (N)).$$

$$N\ge$$

$$(\log (M-J)+\log (K+1)+\log (N))$$

$${\mu }_{max}=\underset{i,j}{\max }\sqrt{N}|{\mathbf{B}}_{ij}|.$$

$$||\hat{\mathbf{X}}-{\mathbf{X}}_0|{|}_F\le \left(C_1+C_2\sqrt{J}\right)\eta$$

$$\frac{N}{{\log }^2(N)}\ge$$

$$(\log (M-J)+\log (MK)+\log (N))$$

$$||\hat{\mathbf{X}}-{\mathbf{X}}_0|{|}_F\le \left(C_1+C_2\sqrt{PJ}\right)\eta$$

$$N\ge$$

$$(\log (M-J)+\log (MK)+\log (N))$$

$$-|⟨{\mathbf{H}}_T,\mathrm{sign}({\mathbf{X}}_0)⟩|+||{\mathbf{H}}_{T^C}|{|}_{2,1}>0$$

$$||{\mathbf{X}}_0+\mathbf{H}|{|}_{2,1}=||{\mathbf{X}}_{0,T}+{\mathbf{H}}_T|{|}_{2,1}+||{\mathbf{H}}_{T^C}|{|}_{2,1}$$

$$\ge |⟨{\mathbf{X}}_{0,T}+{\mathbf{H}}_T,\mathrm{sign}({\mathbf{X}}_{0,T})⟩|+||{\mathbf{H}}_{T^C}|{|}_{2,1}$$

Full text

I Introduction

I-A Overview

111Parts of the results in this paper were presented at the 44th International Conference on Acoustics, Speech, and Signal Processing (ICASSP), 2019 [1].

In classical sparse recovery and compressive sensing problems, a system observes $\bm{y}=\mathbf{D}\mathbf{A}\bm{c}\in\mathbf{C}^{N}$ where $\mathbf{D}$, $\mathbf{A}$, and $\bm{c}$ are the sensing matrix, dictionary matrix, and sparse signal coefficient vector, respectively. The goal is to recover the sparse vector $\bm{c}$ from the measurements $\bm{y}$. Usually $\mathbf{D}$ and $\mathbf{A}$ are known, but the whole system is under-determined. This model arises naturally in a wide range of applications such as medical imaging [2], seismic imaging [3], video coding [4], and network traffic monitoring [5].

In the special case where $\mathbf{D}$ is diagonal and contains a carrier signal or the Fourier coefficients of a known source signal along in its diagonal entries, $\bm{y}$ can be viewed as a modulated version of the signal $\mathbf{A}\bm{c}$ [6] or the Fourier transform of the convolution between two source signals [7]. Recovering $\bm{c}$ can thus be viewed as a demodulation (or deconvolution) problem. Unfortunately, in problems like super resolution [8] and self-calibration [9], the modulation matrix $\mathbf{D}$ is unknown a priori, as it incorporates the unknown point spread functions or calibration parameters. Recovering $\mathbf{D}$ and $\bm{c}$ jointly is a simultaneous sparse recovery and blind demodulation problem.

In this paper, we consider a more general sparse recovery and blind demodulation problem in which each atom comprising the signal undergoes a distinct modulation process. Under the assumption that the modulating waveforms live in a known common subspace, we employ the lifting technique and recast this problem as the recovery of a column-wise sparse matrix from structured linear measurements. In this framework, we recover the sparse coefficient vector $\bm{c}$ and all of the modulating waveforms simultaneously by minimizing the induced atomic norm [10, 11], which in this problem corresponds to the block $\ell_{1}$ norm minimization and we also refer to it as the $\ell_{2,1}$ norm minimization.

I-B Setup and Notation

To better illustrate our main contributions and compare to related work, we first define our signal model and the corresponding atomic norm minimization problem.

Throughout this paper, we use bold uppercase, $\mathbf{X}$, bold lowercase, $\bm{x}$, and non-bold letters, $x$, to represent matrices, vectors, and scalars. We use $\bar{\cdot}$, $\cdot^{H}$ and $\cdot^{T}$ to denote respectively complex conjugate, matrix Hermitian, and matrix transpose. The symbol $C$ denotes a constant. $\mathbf{X}_{T}$ ($\bm{x}_{T}$, resp.) is a matrix (vector, resp.) that zeros out the columns (entries, resp.) not in $T$. We call $T$ the support of the matrix $\mathbf{X}$ (and vector $\bm{x}$), and we use $\tilde{\mathbf{X}}$ to denote the sub-matrix after removing the zero rows or columns in $\mathbf{X}$. $\operatorname*{sign}(\bm{x})=\bm{x}/||\bm{x}||_{2}$ when $||\bm{x}||_{2}\neq 0$ and $\mathbf{0}$ otherwise. $\operatorname*{sign}(\mathbf{X})=[\operatorname*{sign}(\bm{x}_{1}),\cdots,\operatorname*{sign}(\bm{x}_{M})]$. We use $||\cdot||$ to indicate the spectral norm, which returns the maximum singular value of a matrix. The $\ell_{2,1}$ norm of a matrix $\mathbf{X}=[\bm{x}_{1}\quad\cdots\quad\bm{x}_{M}]$, denoted by $||\mathbf{X}||_{2,1}$, is defined to be $\sum_{j=1}^{M}||\bm{x}_{j}||_{2}$. The inner product between vectors and matrices are defined as $\langle\bm{x},\bm{y}\rangle=\bm{y}^{H}\bm{x}$ and $\langle\mathbf{X},\mathbf{Y}\rangle=\operatorname*{Tr\ }\left(\mathbf{Y}^{H}\mathbf{X}\right)$ respectively.

I-C Problem Formulation

In this paper, we study a generalized sparse recovery and blind demodulation problem in which the coefficient vector is unknown and each atom (column) of the dictionary undergoes an unknown modulation process. Specifically, we assume the system receives a composite signal

[TABLE]

where $c_{j}\in\mathbf{C}$ is an unknown scalar, $\mathbf{D}_{j}\in\mathbf{C}^{N\times N}$ is an unknown diagonal modulation matrix, and $\bm{a}_{j}\in\mathbf{C}^{N}$ is the $j$-th atom from a known dictionary $\mathbf{A}=\begin{bmatrix}\bm{a}_{1}&\bm{a}_{2}&\cdots&\bm{a}_{M}\end{bmatrix}\in\mathbf{C}^{N\times M}$ with $N<M$. Our goal is to recover both $c_{j}$ and $\mathbf{D}_{j}$ for all $j$ from the measurement $\bm{y}$.

To make this problem well-posed, among the $M$ over-complete atoms, we assume only $J<M$ of them contribute to the observed signal; that is, at most $J$ coefficients $c_{j}$ are nonzero. We furthermore assume that each modulation matrix obeys a subspace constraint:

[TABLE]

where $\mathbf{B}\in\mathbf{C}^{N\times K}$ ($N>K$) is a known basis for the $K$-dimensional subspace of possible modulating waveforms, and $\bm{h}_{j}\in\mathbf{C}^{K}$ is an unknown coefficient vector. Similar subspace assumptions have been made in deconvolution and demixing papers [12, 13]. With this assumption, recovering $c_{j}$ and $\mathbf{D}_{j}$ equals to recovering $c_{j}$ and $\bm{h}_{j}$. Since $c_{j}\mathbf{D}_{j}\bm{a}_{j}=c_{j}\operatorname*{diag}(\mathbf{B}\bm{h}_{j})\bm{a}_{j}=(kc_{j})\operatorname*{diag}(\mathbf{B}\left(\frac{1}{k}\bm{h}_{j}\right))\bm{a}_{j}$ for any $k\neq 0$, without loss of generality, we assume $\bm{h}_{j}$ has unit norm and $c_{j}\geq 0$ with its complex phase and sign absorbed by $\bm{h}_{j}$.

Define $\mathbf{B}^{H}=[\bm{b}^{\prime}_{1}\quad\bm{b}^{\prime}_{2}\quad\cdots\quad\bm{b}^{\prime}_{N}]\in\mathbf{C}^{K\times N}$ and note that the $n$-th entry of the observed signal can be expressed as

[TABLE]

where $\mathbf{G}=\sum_{j=1}^{M}c_{j}\bm{h}_{j}\bar{\bm{a}}_{j}^{H}$, and $\bm{e}_{n}$ is the $n$-th column of the $N\times N$ identity matrix. From (I.3), we see that the measurement vector $\bm{y}$ depends linearly on the matrix $\mathbf{G}$ which encodes all of the unknown parameters of interest. We denote this linear sensing process as $\bm{y}=\mathcal{L}^{\prime}(\mathbf{G})$ and recast the recovery problem as that of recovering $\mathbf{G}$ (and its components) from the linear measurements.

The unknown matrix $\mathbf{G}$ can be viewed as a linear combination of $J$ rank-$1$ matrices from the atomic set $\mathcal{A}:=\{\bm{h}\bar{\bm{a}}^{H}:\bar{\bm{a}}\in\{\bar{\bm{a}}_{1},...,\bar{\bm{a}}_{M}\},||\bm{h}||_{2}=1\}$ and thus we propose to recover $\mathbf{G}$ using the corresponding atomic norm minimization:

[TABLE]

The atomic norm appearing in (I.4) is defined as $||\mathbf{G}||_{\mathcal{A}}:=\inf\{\sum_{k}|\tilde{c}_{k}|:\mathbf{G}=\sum_{k}\tilde{c}_{k}\mathbf{J}_{k},\mathbf{J}_{k}\in\mathcal{A}\}$. Moreover, the following result establishes its equivalence with the $\ell_{2,1}$ norm.

Proposition 1**.**

The atomic norm optimization problem (I.4) can be equivalently expressed as the following $\ell_{2,1}$ norm optimization problem

[TABLE]

where $\mathbf{X}=[c_{1}\bm{h}_{1}\quad c_{2}\bm{h}_{2}\quad\cdots\quad c_{M}\bm{h}_{M}]\in\mathbf{C}^{K\times M}$ and $\mathcal{L}$ represents the following linear sensing process

[TABLE]

in which $\bm{b}^{\prime}_{n}$ and $\bm{a}_{n}^{\prime}$ are the $n$-th column of $\mathbf{B}^{H}$ and $\mathbf{A}^{T}$.

Proof.

We first note that the atomic norm can be equivalently expressed as $||\mathbf{G}||_{\mathcal{A}}=\inf\{\sum_{j=1}^{M}|c_{j}|:\mathbf{G}=\sum_{j=1}^{M}c_{j}\bm{h}_{j}\bar{\bm{a}}_{j}^{H},||\bm{h}_{j}||_{2}=1\}$. To see this, consider any decomposition of $\mathbf{G}$ of the form $\mathbf{G}=\sum_{k}\tilde{c}_{k}\mathbf{J}_{k}$ with $\mathbf{J}_{k}\in\mathcal{A}$. Define $\mathcal{N}_{j}=\{k:\mathbf{J}_{k}=\tilde{\bm{h}}_{k}\bar{\bm{a}}_{j}^{H}\}$ and write $\mathbf{G}=\sum_{j=1}^{M}(\sum_{k\in\mathcal{N}_{j}}\tilde{c}_{k}\tilde{\bm{h}}_{k})\bar{\bm{a}}_{j}^{H}$. This is equivalent to writing $\mathbf{G}=\sum_{j=1}^{M}c_{j}\bm{h}_{j}\bar{\bm{a}}_{j}^{H}$ where $\bm{h}_{j}=\frac{\sum_{k\in\mathcal{N}_{j}}\tilde{c}_{k}\tilde{\bm{h}}_{k}}{||\sum_{k\in\mathcal{N}_{j}}\tilde{c}_{k}\tilde{\bm{h}}_{k}||_{2}}$ and $c_{j}=||\sum_{k\in\mathcal{N}_{j}}\tilde{c}_{k}\tilde{\bm{h}}_{k}||_{2}$. Finally, note that $|c_{j}|\leq\sum_{k\in\mathcal{N}_{j}}|\tilde{c}_{k}|$.

Next, to establish the equivalence with the $\ell_{2,1}$ norm, for any $c_{j}$ and $\bm{h}_{j}$ with $||\bm{h}_{j}||_{2}=1$, define $\mathbf{x}_{j}=c_{j}\mathbf{h}_{j}$ and $\mathbf{X}=[\mathbf{x}_{1}\quad\mathbf{x}_{2}\quad\cdots\quad\mathbf{x}_{M}]$. Then

[TABLE]

Finally, to establish the equivalence of the linear sensing process, (I.3) indicates that for $\mathbf{G}=\mathbf{X}\bar{\mathbf{A}}^{H}$,

[TABLE]

∎

The above optimization focuses on recovering the structured matrix $\mathbf{X}$ from linear measurements. Once the optimization is solved, the unknown parameters can be easily extracted from the solution $\hat{\mathbf{X}}$ as follows:

[TABLE]

for $\hat{\bm{x}}_{j}\neq 0$ and $1\leq j\leq M$.

The adjoint of the linear operator $\mathcal{L}$ is $\mathcal{L}^{*}(\bm{y})=\sum_{l=1}^{N}y_{l}\bm{b}^{\prime}_{l}\bm{a}^{\prime H}_{l}$. The linear operator $\mathcal{L}$ also has a matrix-vector multiplication form. Note that $\mathcal{L}(\mathbf{X})=\mathbf{\Phi}\cdot\text{vec}(\mathbf{X})$, where $\mathbf{\Phi}\in\mathbf{C}^{N\times KM}$ is

[TABLE]

in which $\mathbf{\phi}_{i,j}=\text{diag}(\bm{b}_{i})\bm{a}_{j}\in\mathbf{C}^{N\times 1}$ and $\bm{b}_{i}$ is the $i$-th column of $\mathbf{B}$. Furthermore,

[TABLE]

where $\mathbf{\phi}^{\prime}_{i}=\bar{\bm{a}}^{\prime}_{i}\otimes\bm{b}_{i}^{\prime}\in\mathbf{C}^{KM\times 1}$.

Finally, we note that the observed signal could be contaminated with noise. In this case, our measurement model becomes

[TABLE]

for some unknown noise vector $\bm{n}\in\mathbf{C}^{N\times 1}$ which we suppose satisfies $||\bm{n}||_{2}\leq\eta$. In this case, we can write $\bm{y}=\mathcal{L}(\mathbf{X}_{0})+\bm{n}$, where $\mathbf{X}_{0}$ is the ground truth solution. As an alternative to equality-constrained $\ell_{2,1}$ norm minimization (I.5), we then consider the following relaxation:

[TABLE]

I-D Applications of The Proposed Signal Model

The proposed signal model encompasses a wide range of applications. We briefly introduce some of them as follows.

I-D1 Direction of arrival estimation for antenna array

We first consider the direction of arrival (DOA) estimation problem in antenna array. Assume we have a linear array antenna consisting of $N$ elements, and we want to estimate the DOAs of several sources from a snapshot of the received signal. In addition, we consider the narrowband scenario and confine the the array and the far-field sources to a common plane as described in [14]. In this case, the DOA is determined by the azimuth angle, $\theta$, of the source, which ranges from [math] to $180$ degrees. Mathematically, after discretizing the azimuth angle into $M$ grids, the observervation of the array can be represented as [15]

[TABLE]

where $\mathbf{D}\in\mathbf{C}^{N\times N}$ is the diagonal matrix capturing the unknown calibration of the array elements [9]. Particularly, the calibration issue may arise from gain discrepancies caused by the change of temperatures and humidity of the environment [9]. Namely, the channel is not ideal. One can simulate different scenarios and collect many possible calibration vectors. By applying the singular value decomposition (SVD) on the matrix formed by those calibration vectors, we can then extract the subspace matrix, $\mathbf{B}$, with desired dimensions to approximate the calibration using $\mathbf{D}=\operatorname*{diag}(\mathbf{B}\bm{h})$ where $\bm{h}$ is the unknown coefficient vector. $\mathbf{A}(\theta)\in\mathbf{C}^{N\times M}$ is the known array manifold matrix whose columns $\bm{a}(\theta_{j})$ for $j\in\{1,2,\cdots,M\}$ are the steering vectors. For uniformly spaced linear array antenna (ULA), $\bm{a}(\theta_{j})=[1,e^{i\frac{2\pi d}{\lambda}\cos(\theta_{j})},\cdots,e^{i(N-1)\frac{2\pi d}{\lambda}\cos(\theta_{j})}]$ where $d$ is the distance between array elements and $\lambda$ is the radar operating wavelength [16]. Moreover, the entries of $\bm{c}$ indicate the strength of the impinging signals and if there exists $J(<M)$ sources, only $J$ entries of $\bm{c}$ are nonzero. $\bm{n}$ consists of the discretization error, approximation error, and additive noise.

Furthermore, let us consider a more severe while realistic situation, where the calibration is sensitive to the direction of arrival which implies that the channel responses from different angles are slightly different. So that the calibration matrix, $\mathbf{D}$, are different for different $\theta_{j}$. In this case, we can write

[TABLE]

I-D2 Super-resolution for single molecule imaging

Another application is the single molecule imaging [17] via stochastic optical reconstruction microscopy (STORM) [18]. In this application, the cellular structure of the object of interest is dyed with fluorophores, and STORM divides the imaging process into thousands of cycles. Within each cycle or observation, only a portion of the fluorophores are activated and imaged. Therefore, a typical observation is a low-resolution frame with its activated fluorophores convolved with the non-stationary point spread functions of the microscope, which can be represented as

[TABLE]

where $\bm{y}\in\mathbf{R}^{N\times 1}$ is a vectorized, imaged frame downsampled from its super-resolution image with $M(>N)$ pixels, $c_{j}$ represents the intensity of the activated fluorophores, and $\mathbf{B}^{\prime}$ is the subspace that the point spread functions live in. $\bm{e}_{j}\in\mathbf{R}^{M}$, which indicates the location of the activated flurophores, is the $j$-th column of the identity matrix and $\bm{n}^{\prime}$ denotes the noise. Moreover, $\bm{y}$ can also be represented equivalently as

[TABLE]

where $IDFT[\cdot]$ is the inverse discrete Fourier transform (DFT) operator, $\mathbf{D}_{j}=\operatorname*{diag}(\mathbf{B}\bm{h}_{j})$ with $\mathbf{B}=DFT[\mathbf{B}^{\prime}]$, and $\bm{a}_{j}$s are the DFT of spikes containing the location information. $\bm{n}=DFT[\bm{n}^{\prime}]$. The goal of this application is to recover the super-resolution image from its low-resolution frame $\bm{y}$, or mathematically, locating the nonzero $c_{j}$.

Other applications that fit into the model investigated in this work include frequency estimation with damping that appears in nuclear magnetic resonance spectroscopy [19] with damping signals approximately living in a common subspace [8] and the CDMA system with spreading sequence sensitive channel as described in Section 6.4 of [9].

I-E Main Contributions

Our contributions are twofold. First, we employ $\ell_{2,1}$ norm minimization to achieve sparse recovery and blind demodulation simultaneously given the generalized signal model from equation (I.1). Second, for perfect recovery of all parameters in the noiseless case, we derive near optimal sample complexity bounds for the cases where $\mathbf{A}$ is a random Gaussian and a random subsampled Fourier dictionary. Both of bounds require the number of measurements $N$ to be proportional to the number of degrees of freedom, $O(JK)$, up to log factors. We also provide bounds on recovering the column-wise sparse matrix in the noisy case; these bounds show that the recovery error scales linearly with respect to the strength of the noise.

I-F Related Work

The $\ell_{2,1}$ norm has been widely used to promote sparse recovery in multiple measurement vector (MMV) problems [20, 21]. The MMV problem involves a collection of sparse signal vectors that are stacked as the rows of a matrix $\mathbf{X}$. These signals have a common sparsity pattern, which results in a column-wise sparse structure for $\mathbf{X}$. As in our setup, the $\ell_{2,1}$ norm is used to recover $\mathbf{X}$ from linear measurements of the form $\bm{y}=\mathbf{\Phi}_{MMV}\cdot\text{vec}(\mathbf{X}^{T})$. However, $\mathbf{\Phi}_{MMV}$ has a block diagonal structure where all diagonal sub-matrices are the same which is the dictionary matrix. This is different from the structure of the linear measurements in our problem; see for example (I.8).

Our work is also closely related to certain recent works in model-based deconvolution, self-calibration, and demixing. When all $\mathbf{D}_{j}$ in (I.1) are the same, our signal model coincides with the self-calibration problem in [9], although that work employs $\ell_{1}$ norm minimization rather than $\ell_{2,1}$ norm minimization to recover $\mathbf{X}$. A more recent paper [22] does apply the $\ell_{2,1}$ norm for the self-calibration problem but again assumes a common modulation matrix $\mathbf{D}$. The paper [12] generalizes the work of [9] and considers a blind deconvolution and demixing problem which can be interpreted as the self-calibration scenario with multiple sensors whose calibration parameters might be different. However, the signal model in that paper is not directly comparable to our model, and the recovery approach studied in that paper involves nuclear norm minimization and requires knowledge of the number of sensors. A blind sparse spike deconvolution is studied in [13], wherein the dictionary consists of sampled complex sinusoids over a continuous frequency range and all atoms undergo the same modulation. Inspired by [13], [8] generalizes the model to the case of different modulating waveforms. Like [13], however, [8] also considers a sampled sinusoid dictionary over a continuous frequency range, and it employs a random sign assumption on the coefficient vectors $\bm{h}_{j}$ which makes it difficult to derive recovery guarantees with noisy measurements. More works considering a common modulation process can be found in [7, 23, 24].

Our work can be viewed as a generalization of the self-calibration [9] and blind deconvolution problems [7]. Moreover, our analysis is quite different from the works considering the continuous sinusoid dictionary [13, 8], since the tools in those papers are specialized to the continuous sinusoids dictionary and we consider discrete Gaussian and random Fourier dictionaries in both noiseless and noisy settings.

The rest of the paper is organized as follows. In Section II, we present our main theorems regarding perfect parameter recovery in the noiseless setting and matrix denoising in the noisy setting. Sections III and IV contain the detailed proofs of the main theorems. Several numerical simulations are provided in Section V to illustrate the critical scaling relationships, and we conclude in Section VI.

II Main Results

We present our main theorems in this section. In each of the noisless and noisy cases, we consider two models for the dictionary matrix $\mathbf{A}$. In the first model, $\mathbf{A}\in\mathbf{R}^{N\times M}$ is a real-valued random Gaussian matrix, with each entry sampled independently from the standard normal distribution. In the second model, $\mathbf{A}\in\mathbf{C}^{N\times M}$ is a complex-valued random Fourier matrix, with each of its $N(<M)$ rows chosen uniformly with replacement from the $M\times M$ discrete Fourier transform matrix $\mathbf{F}$ where $\mathbf{F}^{H}\mathbf{F}=M\mathbf{I}_{M}$. Our first theorem concerns perfect parameter recovery in the noiseless setting.

Theorem II.1**.**

(Noiseless case) Consider the measurement model in equation (I.1), assume that at most $J(<M)$ coefficients $c_{j}$ are nonzero, and furthermore assume that the nonzero coefficients $c_{j}$ are real-valued and positive. Suppose that each modulation matrix $\mathbf{D}_{j}$ satisfies the subspace constraint (I.2), where $\mathbf{B}^{H}\mathbf{B}=\mathbf{I}_{K}$ and each $\bm{h}_{j}$ has unit norm 222Theorem II.1 actually works for $\bm{h}_{j}$ with arbitrary norms as long as the relative scale between $c_{j}$ and $\bm{h}_{j}$ is known. .

Then the solution $\hat{\mathbf{X}}$ to problem (I.5) is the ground truth solution $\mathbf{X}_{0}$—which means that $c_{j}$, $\bm{h}_{j}$, and $\mathbf{D}_{j}$ can all be successfully recovered for each $j$ using (I.7)—with probability at least $1-O(N^{-\alpha+1})$

• •

if $\mathbf{A}\in\mathbf{R}^{N\times M}$ is a random Gaussian matrix and

[TABLE]

• •

if $\mathbf{A}\in\mathbf{C}^{N\times M}$ is a random Fourier matrix and

[TABLE]

where $\gamma=\sqrt{2M\log(2KM)+2M+1}$.

In both cases, $C_{\alpha}$ is a constant defined for $\alpha>1$ and the coherence parameter

[TABLE]

We note that both of the sample complexity bounds in Theorem II.1 require the number of measurements $N$ to be proportional to the number of degrees of freedom, $O(KJ)$, up to log factors. We also note that the sample complexity bounds scale with the square of the coherence parameter $\mu_{max}=\max_{i,j}\sqrt{N}|\mathbf{B}_{ij}|$. Under the assumption $\mathbf{B}^{H}\mathbf{B}=\mathbf{I}_{K}$ which requires the columns of $\mathbf{B}$ to be orthonormal, $\mu_{max}\in[1,\sqrt{N}]$. Specifically, given the system parameters with large enough $N$, (II.1) is satisfied when $1\leq\mu_{max}\leq\sqrt{\frac{N}{C_{\alpha}\log^{2}(N)KJ(\log(M-J)+\log(N))}}$. The valid range of $\mu_{max}$ for (II.2) and the noisy case can be easily derived in the same manner. And $\mu_{max}$ is minimized when the energy of each column of $\mathbf{B}$ is not concentrated on a few entries but spread across the whole column.

Our second theorem provides bounds on recovering the column-wise sparse matrix in the noisy case; these bounds show that the recovery error scales linearly with respect to the strength of the noise.

Theorem II.2**.**

(Noisy case) Consider the measurement model in equation (I.10), assume that at most $J(<M)$ coefficients $c_{j}$ are nonzero, and furthermore assume that the norm of the noise is bounded, $||\bm{n}||_{2}\leq\eta$. Suppose also that each modulation matrix $\mathbf{D}_{j}$ satisfies the subspace constraint (I.2), where $\mathbf{B}^{H}\mathbf{B}=\mathbf{I}_{K}$.

Then with probability at least $1-O(N^{-\alpha+1})$, the solution $\hat{\mathbf{X}}$ to problem (I.11) satisfies

• •

if $\mathbf{A}\in\mathbf{R}^{N\times M}$ is a random Gaussian matrix,

[TABLE]

when

[TABLE]

where $C$ is a constant.

• •

if $\mathbf{A}\in\mathbf{C}^{N\times M}$ is a random Fourier matrix,

[TABLE]

when

[TABLE]

where $\gamma=\sqrt{2M\log(2KM)+2M+1}$ and $P\geq\log(4\sqrt{2J}\gamma)/\log 2$.

In both cases, $C_{\alpha}$ is defined for $\alpha>1$. $C_{1}$ and $C_{2}$ are constant.

Although Theorem II.2 focuses exclusively on bounding the recovery error of the matrix $\mathbf{X}_{0}$, one can also attempt to estimate the parameters $c_{j}$, $\bm{h}_{j}$, and $\mathbf{D}_{j}$ from $\hat{\mathbf{X}}$ using (I.7). And according to Theorem II.2, for any $\hat{\bm{x}}_{j}=\hat{c}_{j}\hat{\bm{h}}_{j}$ and $\bm{x}_{0,j}=c_{0,j}\bm{h}_{0,j}$ where $\hat{\bm{x}}_{j}$ and $\bm{x}_{0,j}$ are the $j$-th columns of the solution $\hat{\mathbf{X}}$ and the ground truth $\mathbf{X}_{0}$ respectively, we would have $||\hat{c}_{j}\hat{\mathbf{D}}_{j}-c_{0,j}\mathbf{D}_{0,j}||_{F}=||\hat{c}_{j}\hat{\bm{h}}_{j}-c_{0,j}\bm{h}_{0,j}||_{2}\leq\left(C_{1}+C_{2}\sqrt{J}\right)\eta$ with random Gaussian dictionary and $||\hat{c}_{j}\hat{\mathbf{D}}_{j}-c_{0,j}\mathbf{D}_{0,j}||_{F}=||\hat{c}_{j}\hat{\bm{h}}_{j}-c_{0,j}\bm{h}_{0,j}||_{2}\leq\left(C_{1}+C_{2}\sqrt{PJ}\right)\eta$ for random Fourier dictionary. In addition, as results on structured matrix recovery from (possibly noisy) linear measurements, we believe that Theorems II.1 and II.2 may be of independent interest outside of the sparse recovery and blind demodulation problem.

III Proof of Theorem II.1

To begin our proof of the main theorem in the noiseless case, we first derive sufficient conditions for exact recovery.

III-A Sufficient Conditions for Exact Recovery

Sufficient conditions for exact recovery are the null space property and an alternative sufficient condition derived from the null space property. Similar sufficient conditions with complete proofs are available for minimization problems using other types of norms [9, 25, 26, 27]. However, since we cannot find sufficient conditions that suit our purpose and in order to be self-contained, we provide a short proof for the ones specific to the $\ell_{2,1}$ norm minimization problem in this section.

Proposition 2**.**

(The null space property) The matrix $\mathbf{X}_{0}=[c_{1}\bm{h}_{1}\quad c_{2}\bm{h}_{2}\quad....\quad c_{M}\bm{h}_{M}]\in\mathbf{C}^{K\times M}$ with support $T$ is the unique solution to the inverse problem (I.5) if

[TABLE]

for any $\mathbf{H}\neq\textbf{0}$ in the nullspace of $\mathcal{L}$.

Proof.

Let $\hat{\mathbf{X}}=\mathbf{X}_{0}+\mathbf{H}$ be a solution to problem (I.5), with $\mathcal{L}(\mathbf{H})=\mathbf{0}$. To prove $\mathbf{X}_{0}$ is the unique solution, it is sufficient to show that $||\hat{\mathbf{X}}||_{2,1}>||\mathbf{X}_{0}||_{2,1}$ if $\mathbf{H}\neq\textbf{0}$. We start by observing that

[TABLE]

where $\operatorname*{sign}(\mathbf{X}_{0,T})=\operatorname*{sign}(\mathbf{X}_{0})$ and the first inequality comes from the fact that

[TABLE]

Therefore, as long as $-|\langle\mathbf{H}_{T},\operatorname*{sign}(\mathbf{X}_{0})\rangle|+||\mathbf{H}_{T^{C}}||_{2,1}>0$ for any $\mathbf{H}\neq\textbf{0}$ in the nullspace of $\mathcal{L}$, $\mathbf{X}_{0}$ is the unique solution. ∎

Proposition 3**.**

The matrix $\mathbf{X}_{0}\in\mathbf{C}^{K\times M}$ with support $T$ is the unique solution to the inverse problem (I.5) if there exists $\gamma>0$ and a matrix $\mathbf{Y}$ in the range space of $\mathcal{L}^{*}$ such that

[TABLE]

and the operator $\mathcal{L}$ satisfies ($\mathcal{L}_{T}(\mathbf{X})=\{\bm{b}^{\prime H}_{n}\mathbf{X}\bm{a}^{\prime}_{n,T}\}_{n=1}^{N}$)

[TABLE]

Proof.

Proposition 2 shows that to establish uniqueness, it is sufficient to prove that $-|\langle\mathbf{H}_{T},\operatorname*{sign}(\mathbf{X}_{0})\rangle|+||\mathbf{H}_{T^{C}}||_{2,1}>0$ for any $\mathbf{H}\neq\textbf{0}$ in the nullspace of $\mathcal{L}$. Note that

[TABLE]

since $\langle\mathbf{H}_{T},\mathbf{Y}_{T}\rangle=-\langle\mathbf{H}_{T^{C}},\mathbf{Y}_{T^{C}}\rangle$. By applying the Hölder inequality, we get a stronger condition

[TABLE]

Since $||\mathcal{L}_{T}^{*}\mathcal{L}_{T}-\mathbf{I}_{T}||\leq\frac{1}{2}$ and $||\mathcal{L}||\leq\gamma$, we have $||\mathcal{L}(\mathbf{H}_{T})||_{F}\geq\frac{1}{\sqrt{2}}||\mathbf{H}_{T}||_{F}$, $||\mathcal{L}(\mathbf{H}_{T^{C}})||_{F}\leq\gamma||\mathbf{H}_{T^{C}}||_{F}$ and

[TABLE]

Plugging (III.3) into the stronger condition above yields

[TABLE]

Therefore, if $||\mathbf{Y}_{T}-\operatorname*{sign}(\mathbf{X}_{0,T})||_{F}\leq\frac{1}{4\sqrt{2}\gamma}$, $||\mathbf{Y}_{T^{C}}||_{2,\infty}\leq\frac{1}{2}$, and $\mathbf{H}_{T^{C}}\neq\textbf{0}$, the left hand side is positive. On the other hand, if $\mathbf{H}_{T^{C}}=\textbf{0}$, from (III.3), $\mathbf{H}_{T}=\textbf{0}$ and $\mathbf{H}=\textbf{0}$. ∎

III-B Bounding The Isometry Constant and Operator Norm

In this section, we bound the isometry constant and operator norm $\gamma$ appearing in (III.2) based on the randomness in the matrix $\mathbf{A}$. The isometry bound for the linear operator $\mathcal{L}$ can be found in Lemma 4.3 in [9].

Lemma III.1**.**

[9, Lemma 4.3]** (Isometry) For the linear operator $\mathcal{L}$ defined in (I.5) with $\mathbf{B}^{H}\mathbf{B}=\mathbf{I}_{K}$ and $\delta>0$,

[TABLE]

with probability at least $1-N^{-\alpha+1}$ where $\mathbf{I}_{T}$ is the identity operator on the support $T$ such that $\mathbf{I}_{T}(\mathbf{X})=\mathbf{X}_{T}$,

• •

if $\mathbf{A}$ is a random Gaussian matrix and $N\geq C_{\alpha}\mu_{max}^{2}KJ\max\{\log(N)/\delta^{2},\log^{2}(N)/\delta\}$.

• •

if $\mathbf{A}$ is a random Fourier matrix and $N\geq C_{\alpha}\mu_{max}^{2}KJ\log(N)/\delta^{2}$.

Here $C_{\alpha}$ is a constant that grows linearly with $\alpha>1$.

$\mathcal{L}_{T}(\mathbf{X})=\mathbf{\Phi}_{T}\cdot\text{vec}(\mathbf{X})$ and $\mathbf{\Phi}_{T}$ can be viewed as $\mathbf{\Phi}$ constructed using $\mathbf{A}_{T}$, whose $i$-th column is zero if $i\in\mathbf{T}^{C}$, following (I.8). Therefore, $\mathbf{\Phi}_{T}\in\mathbf{C}^{N\times KM}$ has many zero columns and removing those zero columns results in $\tilde{\mathbf{\Phi}}_{T}\in\mathbf{C}^{N\times KJ}$. If $||\mathbf{\Phi}^{H}_{T}\mathbf{\Phi}_{T}-\mathbf{I}_{T}||=||\tilde{\mathbf{\Phi}}^{H}_{T}\tilde{\mathbf{\Phi}}_{T}-\tilde{\mathbf{I}}_{T}||\leq\delta<1$, $\tilde{\mathbf{\Phi}}^{H}_{T}\tilde{\mathbf{\Phi}}_{T}$ is invertible and $||(\tilde{\mathbf{\Phi}}^{H}_{T}\tilde{\mathbf{\Phi}}_{T})^{-1}||\leq(1-\delta)^{-1}$ according to Lemma A.12 in [26]. This property will be applied in (III.5) and Theorem III.1. To bound the operator norm of $\mathcal{L}$, we use results from [7] and [9].

Lemma III.2**.**

[7, 9]** For the linear operator $\mathcal{L}$ defined in (I.5) with $\mathbf{B}^{H}\mathbf{B}=\mathbf{I}_{K}$ and $\alpha\geq 1$,

• •

if $\mathbf{A}$ is a random Gaussian matrix,

[TABLE]

with probability at least $1-N^{-\alpha}$.

• •

if $\mathbf{A}$ is a random Fourier matrix,

[TABLE]

with probability at least $1-N^{-\alpha}$ when $N\geq\alpha\mu_{max}^{2}K\log(N)$.

III-C Constructing The Dual Certificate for The Gaussian Case

In the case where $\mathbf{A}$ is a random Gaussian matrix, we construct a certificate matrix $\mathbf{Y}$ that satisfies the conditions in Proposition 3. When $||\mathbf{\Phi}^{H}_{T}\mathbf{\Phi}_{T}-\mathbf{I}_{T}||\leq\frac{1}{2}$, we can set

[TABLE]

where

[TABLE]

By construction, $\mathbf{Y}_{T}=\operatorname*{sign}(\mathbf{X}_{0,T})$, and we need only to verify that $||\mathbf{Y}_{T^{C}}||_{2,\infty}\leq 1/2$.

Theorem III.1**.**

If $||\mathbf{\Phi}^{H}_{T}\mathbf{\Phi}_{T}-\mathbf{I}_{T}||\leq\frac{1}{2}$, there exists $\mathbf{Y}$ in the range space of $\mathcal{L}^{*}$ such that

[TABLE]

with probability at least $1-(M-J)e^{-\alpha}$ when $N\geq 40\alpha\mu_{max}^{2}KJ$ for $\alpha\geq\log(M-J)$.

Proof.

To simplify the notation, without loss of generality, we assume the support of $\mathbf{X}_{0}$ is the first $J$ columns. Let $\mathbf{Y}$ be the dual certificate matrix defined in (III.4). After removing the columns of $\mathbf{Y}$ on support $T$, we obtain $\text{vec}(\tilde{\mathbf{Y}}_{T^{C}})\in\mathbf{C}^{K(M-J)\times 1}$ which takes the form

[TABLE]

The columns of $\tilde{\mathbf{\Phi}}_{T^{C}}$ are independent of $\bm{p}$ since $\bm{p}$ is constructed with $\bm{a}_{i}$ $(i\in T)$. Equivalently,

[TABLE]

Thus $||\mathbf{Y}_{T^{C},j}||_{2}=||\mathbf{P}\bm{a}_{j}||_{2}$ $(j>J)$ where $\bm{a}_{j}$ is real and

[TABLE]

We set $\mathbf{\Sigma}=\mathbf{P}^{H}\mathbf{P}\in\mathbf{C}^{N\times N}$ and have

[TABLE]

since each row of $\mathbf{P}$ can be bounded by

[TABLE]

since we assume $||\mathbf{\Phi}^{H}_{T}\mathbf{\Phi}_{T}-\mathbf{I}_{T}||\leq\frac{1}{2}$ which implies $||(\tilde{\mathbf{\Phi}}_{T}^{H}\tilde{\mathbf{\Phi}}_{T})^{-1}||\leq 2$. By generalizing Proposition 1 in [28] to our case, we have

[TABLE]

In addition, because $\mathbf{\Sigma}$ is positive semi-definite and all its eigenvalues are non-negative, $\operatorname*{Tr\ }(\mathbf{\Sigma}^{2})=\sum_{i=1}^{N}\lambda_{i}^{2}\leq(\sum_{i=1}^{N}\lambda_{i})^{2}=\operatorname*{Tr\ }(\mathbf{\Sigma})^{2}$ where $\lambda_{i}$ is the $i$-th eigenvalue of $\mathbf{\Sigma}$. $||\mathbf{\Sigma}||=\sigma_{max}\leq\sum_{i=1}^{N}\lambda_{i}=\operatorname*{Tr\ }(\mathbf{\Sigma})$ where $\sigma_{max}$ is the maximum singular value of $\mathbf{\Sigma}$. Therefore, for $\alpha>1$, we obtain

[TABLE]

If we pick $N\geq 40\alpha\mu_{max}^{2}KJ$, $||\mathbf{P}\bm{a}_{j}||_{2}>1/2$ with probability at most $e^{-\alpha}$. Taking the union over all $(M-J)$ non-zero columns of $\mathbf{Y}_{T^{C}}$ gives

[TABLE]

Therefore, $||\mathbf{Y}_{T^{C}}||_{2,\infty}\leq 1/2$ with probability at least $1-(M-J)e^{-\alpha}$ when $N\geq 40\alpha\mu_{max}^{2}KJ$. To make the probability meaningful, $\alpha$ should be greater than $\log(M-J)$.

∎

III-D Proof of Theorem II.1 for Random Gaussian Dictionary

In this section, we assemble the pieces to complete the proof of Theorem II.1 in the Gaussian case. To do so, we ensure that all sufficient conditions in Proposition 3 are met. First, if we take $\delta=1/2$ and set $\alpha_{1}>1$ in Lemma III.1, we have

[TABLE]

when $N\geq C_{\alpha_{1}}\mu_{max}^{2}KJ\log^{2}(N)$ with probability at least $1-N^{-\alpha_{1}+1}$. Then, applying the same $\alpha_{1}$ in Lemma III.2 and setting $\gamma=\sqrt{M\log(MN/2)+\alpha_{1}\log(N)}$, we have that $||\mathcal{L}||\leq\gamma$ with probability at least $1-N^{-\alpha_{1}}\geq 1-N^{-\alpha_{1}+1}$. In Theorem III.1, we have proved that $\mathbf{Y}_{T}=\operatorname*{sign}(\mathbf{X}_{0,T})$ and $||\mathbf{Y}_{T^{C}}||_{2,\infty}\leq\frac{1}{2}$ when $N\geq 40\alpha_{2}\mu_{max}^{2}KJ$ with probability at least $1-(M-J)e^{-\alpha_{2}}$ and $\alpha_{2}\geq\log(M-J)$.

Note that if $\alpha_{2}\geq(\alpha_{1}-1)\log(N)+\log(M-J)$, we have $(M-J)e^{-\alpha_{2}}\leq N^{-\alpha_{1}+1}$. Combining the above requirements on $N$, all conditions in Proposition 3 are satisfied with probability at least $1-3N^{-\alpha_{1}+1}$ when $N\geq\max\{C_{\alpha_{1}},40\}((\alpha_{1}-1)\log(N)+\log(M-J))\mu_{max}^{2}KJ\log^{2}(N)$. Furthermore,

[TABLE]

if we set $C_{\alpha}=\max\{C_{\alpha_{1}},40\}\cdot\alpha_{1}$ and $\alpha=\alpha_{1}>1$, which yields the Theorem II.1 when $\mathbf{A}$ is a random Gaussian matrix.

III-E Constructing The Dual Certificate for The Fourier Case

In this section, we construct a certificate $\mathbf{Y}$ that satisfies the inexact duality condition in Proposition 3 when $\mathbf{A}$ is a random Fourier matrix. Specifically, we construct the dual certificate using the golfing scheme [29] which has been widely applied in compressive sensing [7, 25]. In the golfing scheme, a series of matrices in the range of $\mathcal{L}^{*}$ are constructed iteratively. In each iteration step, only some of the measurements are utilized to ensure independence between iterations. And the constructed matrices will converge to $\operatorname*{sign}(\mathbf{X}_{0,T})$ on support $T$ while entries on $T^{C}$ are small. The goal is to find the conditions under which the final constructed matrix can serve as the certificate matrix.

According to Section (4.2.1) in [9], there exists a partition of the $N$ measurements into $P$ disjoint subsets such that each subset, $\Gamma_{p}$, contains $Q$ elements and

[TABLE]

where $\mathbf{B}_{p}=\sum_{l\in\Gamma p}\bm{b}^{\prime}_{l}\bm{b}_{l}^{\prime H}$ and $Q>C\mu_{max}^{2}K\log(N)$. So

[TABLE]

Define $\mathcal{L}_{p}(\mathbf{X})=\{\bm{b}^{\prime H}_{l}\mathbf{X}\bm{a}^{\prime}_{l}\}_{l\in\Gamma_{p}}$ and [math] on entries $l\notin\Gamma_{p}$. $\mathcal{L}_{p}^{*}(\bm{x})=\sum_{l\in\Gamma_{p}}x_{l}\bm{b}^{\prime}_{l}\bm{a}^{\prime H}_{l}$. The golfing scheme iterates through

[TABLE]

Theorem III.2**.**

If $\mathbf{X}_{0}$ is the ground truth solution to problem (I.5), there exists a matrix $\mathbf{Y}\in\mathcal{L}^{*}$ such that

[TABLE]

with probability at least $1-2N^{-\alpha+1}$ for $\alpha>1$ when

[TABLE]

and

[TABLE]

where $C_{\alpha}$ a constant determined by $\alpha$.

Proof.

If we define $\mathbf{W}_{p}=\mathbf{Y}_{p,T}-\text{sign}(\mathbf{X}_{0,T})$, (III.7) gives

[TABLE]

where $\mathcal{L}_{p,T}(\mathbf{X})=\{\bm{b}^{\prime H}_{l}\mathbf{X}\bm{a}^{\prime}_{l,T}\}_{l\in\Gamma_{p}}$ with [math] on entries $l\notin\Gamma_{p}$ and $\mathcal{L}_{p,T}^{*}(\bm{x})=\sum_{l\in\Gamma_{p}}x_{l}\bm{b}^{\prime}_{l}\bm{a}^{\prime H}_{l,T}$ which are used to generate the sequence $\mathbf{Y}_{p,T}$. And we can obtain

[TABLE]

with probability at least $1-N^{-\alpha+1}$ when $Q\geq C_{\alpha,1}\mu_{max}^{2}KJ\log(N)$ with $\alpha>1$ applying Lemma 4.6 in [9]. Therefore,

[TABLE]

To ensure that $||\mathbf{W}_{P}||_{F}=||\mathbf{Y}_{P,T}-\operatorname*{sign}(\mathbf{X}_{0,T})||_{F}\leq\frac{1}{4\sqrt{2}\gamma}$ where $\mathbf{Y}_{P}=\mathbf{Y}$ is the final constructed dual certificate after $P$ iterations, we need

[TABLE]

We now turn to find the conditions such that $||\mathbf{Y}_{T^{C}}||_{2,\infty}\leq\frac{1}{2}$. Note that substituting $\mathbf{W}_{p}$ into equation (III.7) yields

[TABLE]

It is sufficient to show $||\Pi_{T^{C}}(\mathcal{L}_{p}^{*}\mathcal{L}_{p}\mathbf{W}_{p-1})||_{2,\infty}\leq 2^{-p-1}\frac{Q}{N}$, where $\Pi_{T^{C}}$ is the projection operator which projects a matrix on the support $T^{C}$, to make $||\mathbf{Y}_{T^{C}}||_{2,\infty}\leq\frac{1}{2}$ because

[TABLE]

Defining $\mathbf{\Phi}_{p}$ to be $\mathbf{\Phi}$ with non-zero rows indexed by $\Gamma_{p}$ and zero otherwise, we have $\mathcal{L}_{p}(\mathbf{X})=\mathbf{\Phi}_{p}\cdot\text{vec}(\mathbf{X})$ and for a vector $\bm{w}=\text{vec}(\mathbf{W})\in\mathbf{C}^{KM\times 1}$ where $\mathbf{W}\in\mathbf{C}^{K\times M}$ has support $T$,

[TABLE]

where $i$ is the column index and $\bm{e}_{j}$ is the $j$-th column of the identity matrix $\mathbf{I}_{KM}$. In addition,

[TABLE]

Furthermore, we have $\mathbf{E}(\bm{z}_{l,i})=0$ because

[TABLE]

following $\mathbf{E}(\bar{\bm{a}}^{\prime}_{l}\bar{\bm{a}}^{\prime H}_{l})=\mathbf{I}_{M}$ since $\bm{a}^{\prime}_{l}\in\mathbf{C}^{M\times 1}$ is the transpose of a random row of the $M\times M$ DFT matrix and $\bm{b}^{\prime}_{l}\bm{b}_{l}^{\prime H}\mathbf{W}$ has support $T$ and $\mathbf{0}$ on $T^{C}$. Therefore, for $i\in T^{C}$, $\mathbf{E}(\bm{z}_{l,i})=0$. Moreover,

[TABLE]

Because each entry of $\bm{z}_{l,i}$ can be bounded by

[TABLE]

where the third equality holds because $\bm{w}=\text{vec}(\mathbf{W})$ and $\mathbf{W}$ has support T. The variance of $\bm{z}_{l,i}$ is also bounded:

[TABLE]

because for each element of $||\bm{z}_{l,i}||_{2}^{2}$, we have

[TABLE]

and therefore

[TABLE]

As a result,

[TABLE]

The second inequality in (III.13) applies the inequality (III.6) and $||\mathbf{I}_{M}\otimes\mathbf{B}_{p}||=||\mathbf{I}_{M}||\cdot||\mathbf{B}_{p}||$. We then apply the matrix Bernstein inequality from Theorem 1.6 in [30]. If we set $\bm{w}=\text{vec}(\mathbf{W}_{p-1})$ and we know from (III.10) that $||\bm{w}||_{2}=||\mathbf{W}_{p-1}||_{F}\leq 2^{-p+1}\sqrt{J}$, we obtain

[TABLE]

where $t=2^{-p-1}\frac{Q}{N}$, for a particular $i\in T^{C}$ and $p$. We then take the union over all $i\in T^{C}$ and get

[TABLE]

To ensure $||\Pi_{T^{C}}(\mathcal{L}_{p}^{*}\mathcal{L}_{p}(\mathbf{W}_{p-1}))||_{2,\infty}\leq 2^{-p-1}\frac{Q}{N}$ for all $p$, we obtain

[TABLE]

when $Q\geq\frac{128\mu_{max}^{2}KJ\alpha}{3}(\log(M-J)+\log(K+1)+\log(N))$ using the same $\alpha$ as in deriving equation (III.9). Setting $C_{\alpha}=\max\{C,C_{\alpha,1},\frac{128}{3}\alpha\}$, where $C$ is a constant comes from equation (III.6), gives us Theorem III.2. ∎

III-F Proof of Theorem II.1 for Random Fourier Dictionary

We now complete the proof of Theorem II.1 in the case when $\mathbf{A}$ is a random Fourier matrix. First, combining the conditions and probabilities from Lemma III.1 and III.2, we know that the operator $\mathcal{L}$ satisfies the inequalities $||\mathcal{L}_{T}^{*}\mathcal{L}_{T}-\mathbf{I}_{T}||\leq\frac{1}{2}$ and $||\mathcal{L}||\leq\gamma=\sqrt{2M\log(2KM)+2M+1}$ with probability at least $1-(N+1)N^{-\alpha}\geq 1-2N^{-\alpha+1}$ when $N\geq C_{\alpha,1}\mu_{max}^{2}KJ\log(N)$ for some constant, $C_{\alpha,1}$, that grows linearly with $\alpha>1$.

Applying the same $\alpha$ in Theorem III.2, the desired dual matrix exists with probability at least $1-2N^{-\alpha+1}$ when $N\geq C_{\alpha,2}\mu_{max}^{2}KJ\log(4\sqrt{2J}\gamma)(\log(M-J)+\log(K+1)+\log(N))$. Merging the requirement on $N$ by setting $C_{\alpha}=\max\{C_{\alpha,1},C_{\alpha,2}\}$ and combining the probabilities, we complete the proof by applying Proposition 3.

IV Proof of Theorem II.2

To derive our recovery guarantee in the presence of measurement noise, the main ingredient of the proof is Theorem IV.1 which is a variation of the Theorem 4.33 in [26] from the infinity norm optimization to $\ell_{2,1}$ norm optimization problem.

Theorem IV.1**.**

Define $\mathbf{\Phi}\in\mathbf{C}^{N\times KM}$ and $\mathbf{\Phi}\cdot\text{vec}(\mathbf{X})=\mathcal{L}(\mathbf{X})$. Suppose the ground truth $\mathbf{X}_{0}$ to (I.11) has $J$ non-zero columns with support $T$ and the measurement vector $\bm{y}=\mathcal{L}(\mathbf{X}_{0})+\bm{n}$ with $||\bm{n}||_{2}\leq\eta$. For $\delta,\beta,\theta,\gamma,\tau>0$ and $\delta<1$, assume that

[TABLE]

and that there exists a matrix $\mathbf{Y}=\mathcal{L}^{*}(\bm{p})\in\mathbf{C}^{K\times M}$ such that

[TABLE]

If $\rho:=\theta+\frac{\beta}{4\sqrt{2}\gamma(1-\delta)}<1$, then the minimizer, $\hat{\mathbf{X}}$, to (I.11) satisfies

[TABLE]

where $C_{1}$ and $C_{2}$ are two constants depending on $\delta,\beta,\theta,\gamma,\tau$.

Proof.

Due to our assumption on the noise, $\mathbf{X}_{0}$ is a feasible solution. Assume the final minimizer to (I.11) is $\hat{\mathbf{X}}=\mathbf{X}_{0}+\mathbf{H}$, which implies

[TABLE]

where the second inequality comes from equation (III.1). Thus

[TABLE]

The last inequality comes from the Hölder inequality and our assumption $||\bm{n}||\leq\eta$, which tells us

[TABLE]

and

[TABLE]

Moreover, $||\mathbf{H}_{T}||_{F}$ can also be bounded as follows.

[TABLE]

because $||\mathbf{\Phi}_{T}^{H}\mathbf{\Phi}_{T}-\mathbf{I}_{T}||\leq\delta$ ensures that $||(\tilde{\mathbf{\Phi}}_{T}^{H}\tilde{\mathbf{\Phi}}_{T})^{-1}||\leq\frac{1}{1-\delta}$ and $||\mathbf{\Phi}_{T}^{H}||\leq\sqrt{1+\delta}$ according to Lemma A.12 and Proposition A.15 in [26] respectively. Furthermore,

[TABLE]

in which $\bm{h}_{i}$ is the $i$-th column of $\mathbf{H}$. By setting $\rho=\theta+\frac{\beta}{4\sqrt{2}\gamma(1-\delta)}$, $\mu=\frac{\sqrt{1+\delta}}{1-\delta}$ and substituting the inequality (IV.2) into (IV.1), we obtain

[TABLE]

Substituting inequality (IV.3) into (IV.2) yields

[TABLE]

Combining the above two inequalities, we obtain

[TABLE]

∎

Next, we specify the values of the variables $\theta$, $\tau$, $\delta$ and $\beta$ when $\mathbf{A}$ is a random Gaussian and Fourier matrix. The Orlicz-1 norm [7] and associated matrix Bernstein inequality are needed for determining the value of $\beta$ when $\mathbf{A}$ is Gaussian. Specifically, the Orlicz-1 norm is defined as [7]

[TABLE]

Its associated matrix Bernstein inequality is provided in Proposition 3 in [7] which can be rewritten as

Proposition 4**.**

Let $\mathbf{Z}_{1},...,\mathbf{Z}_{N}$ be independent $M\times M$ random matrices with $\mathbf{E}(\mathbf{Z}_{j})=0$. Suppose

[TABLE]

and define

[TABLE]

Then there exists a constant $C$ such that for $t>0$

[TABLE]

The following theorem utilizes the Proposition 4 and depicts the conditions under which $\beta=1$.

Theorem IV.2**.**

For $\mathbf{\Phi}$ defined in (I.8) and $\mathcal{L}(\mathbf{X})=\mathbf{\Phi}\text{vec}(\mathbf{X})$,

[TABLE]

with probability at least $1-N^{-\alpha+1}$

• •

if $\mathbf{A}$ is a random Gaussian matrix and

[TABLE]

• •

if $\mathbf{A}$ is a random Fourier matrix and

[TABLE]

where $C_{\alpha}$ is a constant that grows linearly with $\alpha>1$ and $C$ is a constant.

Proof.

We first prove the Gaussian case; the Fourier case is very similar. Note that for an arbitrary $i\in T^{C}$

[TABLE]

where $\mathbf{\Phi}_{i}\in\mathbf{C}^{N\times KM}$ is $\mathbf{\Phi}$ but only contains values in the $\left(K(i-1)+1\right)$-th to $\left(K(i-1)+K\right)$-th columns and is zero otherwise. $\mathbf{\Phi}_{i}$ can also be viewed as an extension of $[\mathbf{\Phi}_{T^{C},K(i-1)+1}\cdots\quad\mathbf{\Phi}_{T^{C},K(i-1)+K}]$ by padding zero columns. Moreover, $\bar{\bm{a}}_{j,i}^{\prime}$ is the conjugate of the $j$-th column of $\mathbf{A}^{T}$ who has only one non-zero value in the $i$-th entry. In addition, $\mathbf{E}(\mathbf{Z}_{j})=\mathbf{E}(\bar{\bm{a}}^{\prime}_{j,T}\bar{\bm{a}}_{j,i}^{\prime H}\otimes\bm{b}^{\prime}_{j}\bm{b}_{j}^{\prime H})=\mathbf{E}(\bar{\bm{a}}^{\prime}_{j,T}\bar{\bm{a}}_{j,i}^{\prime H})\otimes\bm{b}^{\prime}_{j}\bm{b}_{j}^{\prime H}=\mathbf{0}$ for $i\in T^{C}$. By applying the property of the Kronecker product, we estimate the spectral norm of $\mathbf{Z}_{j}$ which can be used to determine its Orlicz-1 norm:

[TABLE]

in which $\bar{\bm{a}}^{\prime}_{j,\{T,i\}}$ contains non-zero values in the entries indexed by $\{T,i\}$. Therefore, $||\bar{\bm{a}}^{\prime}_{j,\{T,i\}}||_{2}^{2}$ follows the Chi-squared distribution with $J+1$ degrees of freedom which implies that $||\mathbf{Z}_{j}||_{\psi_{1}}\leq\frac{C\mu_{max}^{2}K(J+1)}{2N}\leq\frac{C\mu_{max}^{2}K\cdot 2J}{2N}=\frac{C\mu_{max}^{2}KJ}{N}=R$ for some constant $C$ according to the proof of Lemma 4.7 in [9] and the definition of Orlicz-1 norm in (IV.5). Moreover,

[TABLE]

following from the fact that $\mathbf{E}\left(\bar{\bm{a}}^{\prime}_{j,i}\bar{\bm{a}}_{j,T}^{\prime H}\bar{\bm{a}}^{\prime}_{j,T}\bar{\bm{a}}_{j,i}^{\prime H}\right)=J\mathbf{I}_{M,i}$ for all $j$ and $\sum_{j=1}^{N}\bm{b}_{j}^{\prime}\bm{b}_{j}^{\prime H}=\mathbf{I}_{K}$ from the assumption. On the other hand,

[TABLE]

Therefore, $\max\{||\sum_{j=1}^{N}\mathbf{E}(\mathbf{Z}_{j}\mathbf{Z}_{j}^{H})||,||\sum_{j=1}^{N}\mathbf{E}(\mathbf{Z}_{j}^{H}\mathbf{Z}_{j})||\}=\frac{\mu_{max}^{2}KJ}{N}=\sigma^{2}$. Substituting the variables $R$ and $\sigma^{2}$ into Proposition 4 and taking the union bound over all $i\in T^{C}$ results in

[TABLE]

Define a variable $\alpha>1$ and set

[TABLE]

where $C_{\alpha}=C_{0}\alpha$. Simplifying the probability term gives

[TABLE]

Following the same procedures, when $\mathbf{A}$ is a random Fourier matrix and for any $i\in T^{C}$, we have $\mathbf{E}(\mathbf{Z}_{j})=\mathbf{E}(\bar{\bm{a}}^{\prime}_{j,T}\bar{\bm{a}}_{j,i}^{\prime H})\otimes\bm{b}^{\prime}_{j}\bm{b}_{j}^{\prime H}=0$, $||\mathbf{Z}_{j}||=\frac{\mu_{max}^{2}K}{N}||\bar{\bm{a}}^{\prime}_{j,T}||_{2}||\bar{\bm{a}}_{j,i}^{\prime}||_{2}=\frac{\mu_{max}^{2}K\sqrt{J}}{N}=R$ and $\sigma^{2}=\frac{\mu_{max}^{2}KJ}{N}$. The matrix Bernstein inequality implies

[TABLE]

Similarly, if we define a variable $\alpha>1$ and let

[TABLE]

by setting $C_{\alpha}=\frac{8}{3}\alpha$, simplifying the probability gives us

[TABLE]

∎

IV-A Proof of Theorem II.2 for Random Gaussian Dictionary

According to Section III-D, $||\mathbf{\Phi}_{T}^{H}\mathbf{\Phi}_{T}-\mathbf{I}_{T}||\leq\frac{1}{2}=\delta$, $||\mathbf{Y}_{T^{C}}||_{2,\infty}\leq\frac{1}{2}=\theta$ and $\gamma=\sqrt{M\log(MN/2)+\alpha\log(N)}$ with probability at least $1-3N^{-\alpha+1}$ when $\frac{N}{\log^{2}N}\geq C_{\alpha,1}\mu_{max}^{2}KJ(\log(N)+\log(M-J))$. Moreover, in Theorem III.1, where we construct the dual certificate matrix when $\mathbf{A}$ is a random Gaussian matrix, we define $\bm{p}=\tilde{\mathbf{\Phi}}_{T}(\tilde{\mathbf{\Phi}}_{T}^{H}\tilde{\mathbf{\Phi}}_{T})^{-1}\text{vec}(\operatorname*{sign}(\tilde{\mathbf{X}}_{0,T}))\in\mathbf{C}^{N\times 1}$ and $||\mathbf{\Phi}_{T}^{H}\mathbf{\Phi}_{T}-\mathbf{I}_{T}||\leq\frac{1}{2}$ leads to $||(\tilde{\mathbf{\Phi}}_{T}^{H}\tilde{\mathbf{\Phi}}_{T})^{-1}||\leq 2$. So

[TABLE]

which implies $\tau=\sqrt{2}$. If we use the same $\alpha$ in Theorem IV.2, we have $\beta=1$ with probability at least $1-N^{-\alpha+1}$ when

[TABLE]

Combining the requirement on $N$ and setting $C_{\alpha}=\max\{C_{\alpha,1},C_{\alpha,2}\}$ yield

[TABLE]

Therefore, the conditions in Theorem IV.1 are satisfied with probability at least $1-4N^{\alpha+1}$ when $N$ is as defined in equation (IV.6). In addition, after substituting the parameters $\rho=\theta+\frac{\beta}{4\sqrt{2}\gamma(1-\delta)}=\frac{1}{2}+\frac{1}{2\sqrt{2}\gamma}<1$ and $\mu=\frac{\sqrt{1+\delta}}{1-\delta}=\sqrt{6}$ into (IV.4), $2\mu+\frac{\mu}{2\sqrt{2}\gamma(1-\rho)}+\frac{\beta\mu}{2\sqrt{2}\gamma(1-\delta)(1-\rho)}=2\sqrt{6}+\frac{3\sqrt{6}}{\sqrt{2}\gamma-1}\leq 5\sqrt{6}=C_{1}$ and $\frac{2\tau}{1-\rho}+\frac{2\beta\tau}{(1-\delta)(1-\rho)}=\frac{24\gamma}{\sqrt{2}\gamma-1}\leq 24=C_{2}$.

IV-B Proof of Theorem II.2 for Random Fourier Dictionary

In the proof of Theorem III.2, we have derived $\mathbf{Y}=-\frac{N}{Q}\sum_{p=1}^{P}\mathcal{L}_{p}^{*}\mathcal{L}_{p}(\mathbf{W}_{p-1})$. Since the sets $\Gamma_{p}$ are disjoint, the indices of non-zero entries of $\mathcal{L}_{p}(\mathbf{W}_{p-1})$ for different $p$ are disjoint and $\mathbf{Y}=\mathcal{L}^{*}(-\frac{N}{Q}\sum_{p=1}^{P}\mathcal{L}_{p}(\mathbf{W}_{p-1}))=\mathcal{L}^{*}(\bm{p})$. Moreover, $\mathbf{W}_{p-1}$ has support $T$ from its definition in (III.8) which gives us

[TABLE]

because $||\mathbf{\Phi}_{p,T}^{H}\mathbf{\Phi}_{p,T}||\leq\frac{3Q}{2N}$ and $||\mathbf{W}_{p-1}||_{F}^{2}\leq 4^{-p+1}J$ following from Lemma 4.6 in [9] and equation (III.10) respectively. $\mathbf{\Phi}_{p,T}$ is $\mathbf{\Phi}$ constructed with $\mathbf{A}_{T}$ and only rows indexed by $\Gamma_{p}$ are non-zero. Therefore, $||\bm{p}||_{2}\leq\sqrt{2PJ}$ and $\tau=\sqrt{2P}$ with $P\geq\log(4\sqrt{2J}\gamma)/\log 2$ defined in equation (III.11). In addition, from Section III-F and Theorem II.1, we have $\delta=\frac{1}{2}$, $\theta=\frac{1}{2}$ and $\gamma=\sqrt{2M\log(2KM)+2M+1}$ with probability at least $1-4N^{-\alpha+1}$ when

[TABLE]

Applying the same $\alpha$ to Theorem IV.2, $\beta=1$ with probability at least $1-N^{-\alpha+1}$ when $N\geq C_{\alpha,2}\mu_{max}^{2}KJ\left(\log(KM)+\log(M-J)+\log(N)\right)$. One can easily examine that $\rho=\theta+\frac{\beta}{4\sqrt{2}\gamma(1-\delta)}=\frac{1}{2}+\frac{1}{2\sqrt{2}\gamma}<1$.

If we set $C_{\alpha}=\max\{C_{\alpha,1},C_{\alpha,2}\}$ and merge the requirements on $N$, we obtain

[TABLE]

Thus, the conditions in Theorem IV.1 are satisfied with probability at least $1-5N^{-\alpha+1}$ when $N$ satisfies (IV.7). Moreover, since $\mu=\frac{\sqrt{1+\delta}}{1-\delta}$, $2\mu+\frac{\mu}{2\sqrt{2}\gamma(1-\rho)}+\frac{\beta\mu}{2\sqrt{2}\gamma(1-\delta)(1-\rho)}=2\sqrt{6}+\frac{3\sqrt{6}}{\sqrt{2}\gamma-1}\leq 5\sqrt{6}=C_{1}$ and $\frac{2\tau}{1-\rho}+\frac{2\beta\tau}{(1-\delta)(1-\rho)}=\frac{24\gamma\sqrt{P}}{\sqrt{2}\gamma-1}\leq 24\sqrt{P}=C_{2}\sqrt{P}$ with $P\geq\log(4\sqrt{2J}\gamma)/\log 2$.

V Numerical Simulations

Here we present numerical simulations that illustrate and support our theoretical results. We set $\mathbf{B}\in\mathbf{C}^{N\times K}$ to be the first $K$ columns of the normalized DFT matrix $\frac{1}{\sqrt{N}}\mathbf{F}\in\mathbf{C}^{N\times N}$. The parameters $c_{j}$ and $\bm{h}_{j}$ are generated by sampling independently from the standard normal distribution, and the $J$ non-zero columns of the ground truth solution $\mathbf{X}_{0}=[c_{j}\bm{h}_{j}\cdots c_{M}\bm{h}_{M}]$ are selected uniformly. 40 simulations are run for each setting, based on which we compute the percentage of successful recovery. Both the dictionary, $\mathbf{A}$, and the ground truth solution, $\mathbf{X}_{0}$, including the support and its content, are sampled independently for each simulation. We solve problems (I.5) and (I.11) via CVX [31], and in the noiseless case if the relative error between the solution $\hat{\mathbf{X}}$ and the ground truth $\mathbf{X}_{0}$ is smaller than $10^{-5}$, $\frac{||\hat{\mathbf{X}}-\mathbf{X}_{0}||_{F}}{||\mathbf{X}_{0}||_{F}}\leq 10^{-5}$, we count it as a successful recovery.

V-A * The Sufficient Number of Measurement*

In the first noiseless simulation, we examine the recovery rate with respect to the parameters $K$ and $J$. We fix $M=200$ and $N=100$ and let $K$ and $J$ range from $1$ to $20$. The results are summarized in the phase transition plots of Fig. 1 for the random Gaussian dictionary and Fig. 2 for the random Fourier dictionary. The results for the two dictionaries are similar. The reciprocal nature of the phase transition boundary supports the linear scaling with $KJ$ in equations (II.1) and (II.2). Roughly when $KJ\leq 60$, the recovery success rate is satisfactory.

To further illustrate the linear scaling of the required number of measurements $N$ with respect to $K$ and $J$, we fix $M=200$ and $K=5$, and let $N$ and $J$ range from $30$ to $100$ and $1$ to $20$, respectively. The results are recorded in Figs. 3 and 4 for the random Gaussian and Fourier dictionaries, respectively. The same simulation but switching the roles of $K$ and $J$ is also implemented, and the results are shown in Figs. 5 and 6. These results support the linear scaling of Theorem II.1.

V-B The Recovery Error Bound with Noisy Measurement

To test the noisy case, we set $M=200$, $K=J=5$, and $N=100$, and we let $\bm{y}=\mathcal{L}(\mathbf{X}_{0})+\bm{n}$ with $||\bm{n}||_{2}\leq\eta$. Theorem II.2 gives a recovery guarantee of the form $||\hat{\mathbf{X}}-\mathbf{X}_{0}||_{F}\leq C\cdot\eta$ for a constant $C$ . Therefore, after dividing both sides by $||\mathbf{X}_{0}||_{F}$, setting $||\bm{n}||_{2}=\eta$ and changing the units to decibels (dB), we obtain

[TABLE]

We call $20\log_{10}\left(\frac{||\hat{\mathbf{X}}-\mathbf{X}_{0}||_{F}}{||\mathbf{X}_{0}||_{F}}\right)$ the relative error in dB and $20\log_{10}\left(\frac{||\bm{n}||_{2}}{||\mathbf{X}_{0}||_{F}}\right)$ the noise-to-signal ratio in dB. To examine the linear relation between the relative error and the noise-to-signal ratio in equation (V.1), we sample the real and complex components of the noise vector $\bm{n}$ independently from a standard Gaussian distribution and scale $||\bm{n}||_{2}$ to attain different noise-to-signal ratios. Similar to the previous plots, 40 independent simulations are run for each noise-to-signal ratio and the range of the standard deviation and mean (computed before transforming to dB) of the relative error in dB are recorded in Figs. 7 and 8. The dashed lines show the theoretical error bound from Theorem II.2 which are drawn by substituting the constants derived in Section IV-A and IV-B and the system parameters into equations (II.3) and (II.5). The slope of each dashed line are 1. We observe that when noise-to-signal ratio is smaller than 0 dB, the relative error scales linearly with respect to the noise-to-signal ratio with slope 1 for both random Gaussian and Fourier dictionaries. This confirms that $||\hat{\mathbf{X}}-\mathbf{X}_{0}||_{F}$ grows linearly with respect to $\eta$ in Theorem II.2. Moreover, if the noise dominates the observed signal, solving the problem (I.11) results in $\hat{\mathbf{X}}=0$ and the relative error becomes 0 dB.

V-C Direction of Arrival Estimation

In this section, we apply the proposed signal model to the direction of arrival estimation problem introduced in Section I-D1. Note that there exits thousands of different subspaces that the complex calibration could live in. To give a concrete example and compare to the related work, we adopt the setting from [9] where the calibration subspace $\mathbf{B}\in\mathbf{C}^{N\times K}$ is modeled by the first $K$ columns of the normalized DFT matrix $\frac{1}{\sqrt{N}}\mathbf{F}\in\mathbf{C}^{N\times N}$. The entries of $\bm{h}_{j}$ are sampled independently from the standard normal distribution and $\bm{h}_{j}$ is normalized to have unit norm. Moreover, we set $M=181$ and discretize the direction of arrival into $\theta_{j}=\{0,1,\cdots,180\}$ degrees. When the distance between array elements is half of the operating wavelength, we can obtain $\mathbf{A}$ by substituting $d=\frac{\lambda}{2}$ and $\theta_{j}$ into $\bm{a}(\theta_{j})$ defined in Section I-D1. Furthermore, we set $N=50$ and $K=J=5$ where the directions of arrival of the 5 sources are $\{67,75,92,127,133\}$ degrees and the signal magnitudes are sampled independently from the uniform distribution on $[0,1]$. The real and imaginary parts of the noise vector are independent random Gaussian vectors with [math] mean and identity covariance matrix. $SNR=30$ dB. By solving the $\ell_{2,1}$ norm minimization problem in (I.11), the index of the nonzero column in the solution $\hat{\mathbf{X}}$ indicates the direction of arrival and the norm of the nonzero column indicates the signal strength. The result is recorded in Fig .9 (a). As a comparison, we also apply the Sparselift method proposed in [9] to this problem, which assumes $\mathbf{D}_{j}$ for all $j$ are the same and solves an $\ell_{1}$ norm minimization problem. The result is recorded in Fig .9 (b).

V-D Single Molecule Imaging

Furthermore, we apply the proposed signal model to the single molecule imaging described in Section I-D2. All data comes from the Single-Molecule Localization Microscopy grand challenge organized by ISBI [32] which contains 12,000 low-resolution frames. Each low-resolution frame is 64 pixel $\times$ 64 pixel with pixel size 100 nm $\times$ 100 nm, so that $N=64\times 64=4096$. A typical, observed frame is shown in Fig. 10 (a). Superimposing all the observed frames leads to the low-resolution structure in Fig. 10 (b). The target of this experiment is to recover the high resolution image of size 320 pixel $\times$ 320 pixel, which implies that $M=320\times 320=102400$, whose pixel is of size 20 nm $\times$ 20 nm. In addition, according to the statistic of the dataset, the number of activated fluorophores in each frame is less or equal to $J=17$ and we use the Gaussian point spread functions to approximate the point spread functions of the microscope. By implementing the SVD on the Gaussian point spread functions with different variances, we obtain a $K=3$ dimension subspace that point spread functions live in. Then by solving an $\ell_{2,1}$ norm regularized least square minimization problem on each low-resolution frame, we get totally 12,000 high resolution images and superimposing all the high resolution images results in the super-resolution output recorded in Fig. 10 (c).

VI Conclusion

In this paper, we introduce the generalized sparse recovery and blind demodulation model and achieve sparse recovery and blind demodulation simultaneously. Under the assumption that the modulating waveforms live in a known common subspace, we employ the lifting technique and recast this problem as the recovery of a column-wise sparse matrix from structured linear measurements. In this framework, we accomplish sparse recovery and blind demodulation simultaneously by minimizing the induced atomic norm, which in this problem corresponds to $\ell_{2,1}$ norm minimization. In the noiseless case, we derive near optimal sampling complexity that is proportional to the number of degrees of freedom, and in the noisy case we bound the recovery error of the structured matrix. Numerical simulations support our theoretical results. In addition to extending the class of dictionaries we have considered, an interesting future direction would be to relax the constraint that each $\mathbf{D}_{j}$ is diagonal while preserving the low-dimensional subspace assumption.

Acknowledgment

This work was supported by NSF grant CCF$-1704204$.