Blind Super-resolution via Projected Gradient Descent

The point source signal can be represented as a superposition of $r$ spikes

where $d_{k}$ and $\tau_{k}$ are the amplitude and location of $k$-th point source respectively. Let $\{g_{k}(t)\}_{k=1}^{r}$ be the unknown point spread functions. The observation is a convolution between $x(t)$ and $\{g_{k}(t)\}_{k=1}^{r}$,

After taking Fourier transform and sampling, we obtain the measurements as

Let $\bm{g}_{k}=\begin{bmatrix}\hat{g}_{k}[0]&\cdots&\hat{g}_{k}[n-1]\end{bmatrix}^{\mathsf{T}}$ be a vector corresponding to the $k$-th unknown point spread function. The goal is to estimate $\{d_{k},\tau_{k}\}_{k=1}^{r}$ as well as $\{\bm{g}_{k}\}_{k=1}^{r}$ from (2.1). Since the number of unknowns is larger than $n$, this problem is an ill-posed problem without any additional assumptions. Following the same route as that in [4, 5, 6, 8], we assume $\{\bm{g}_{k}\}_{k=1}^{r}$ belong to a known subspace spanned by the columns of $\bm{B}\in\mathbb{C}^{n\times s}$, i.e.,

Then under the subspace assumption and applying the lifting technique [9], the measurements (2.1) can be rewritten as a linear observations of $\bm{X}^{\natural}:=\sum_{k=1}^{r}d_{k}\bm{h}_{k}\bm{a}_{\tau_{k}}^{\mathsf{T}}\in\mathbb{C}^{s\times n}$:

where $\bm{b}_{j}$ is the $j$-th row of $\bm{B}$, $\bm{e}_{j}$ is the $(j+1)$-th standard basis of $\mathbb{R}^{n}$, and $\bm{a}_{\tau}\in\mathbb{C}^{n}$ is the vector defined as

The measurement model (2.2) can be rewritten succinctly as

where ${\cal A}:\mathbb{C}^{s\times n}\rightarrow\mathbb{C}^{n}$ is the linear operator. Therefore, blind super-resolution can be cast as the problem of recovering the data matrix $\bm{X}^{\natural}$ from its linear measurements (2.3).

Let ${\cal H}$ be the vectorized Hankel lifting operator which maps a matrix $\bm{X}\in\mathbb{C}^{s\times n}$ into an $sn_{1}\times n_{2}$ matrix,

where $\bm{x}_{i}$ is the $(i+1)$-th column of $\bm{X}$ and $n_{1}+n_{2}=n+1$. It is shown that ${\cal H}(\bm{X}^{\natural})$ is a rank-$r$ matrix [8] and thus the matrix ${\cal H}(\bm{X}^{\natural})$ admits low rank structure when $r\ll\min(sn_{1},n_{2})$. It is natural to recover $\bm{X}$ by solving the constrained least squares problem

To introduce our algorithm, we assume that $\bm{X}^{\natural}$ is $\mu_{1}$-incoherent which is defined as below.

Letting $\bm{L}^{\natural}=\bm{U}\bm{\Sigma}^{1/2}$ and $\bm{R}^{\natural}=\bm{V}\bm{\Sigma}^{1/2}$, we have

where $\sigma_{1}=\left\|{\cal H}(\bm{X}^{\natural})\right\|$. Since that target data matrix $\bm{X}^{\natural}$ is $\mu_{1}$-incoherence and the low rank structure of the vectorized Hankel matrix can be promoted by ${\cal H}(\bm{X}^{\natural})=\bm{L}^{\natural}{\bm{R}^{\natural}}^{\mathsf{H}}$, it is natural to recover the low rank factors of the ground truth matrix ${\cal H}(\bm{X}^{\natural})$ by solving an optimization problem in the form of

where $\bm{M}=\begin{bmatrix}\bm{L}^{\mathsf{H}}&\bm{R}^{\mathsf{H}}\end{bmatrix}^{\mathsf{H}}\in\mathbb{C}^{(sn_{1}+n_{2})\times r}$, ${\cal H}^{\dagger}$ is the Moore-Penrose pseudoinverse of ${\cal H}$ obeying that ${\cal H}^{\dagger}{\cal H}={\cal I}$, the second term in objective guarantees that $\bm{L}\bm{R}^{\mathsf{H}}$ admits vectorized Hankel structure. Last term penalizes the mismatch between $\bm{L}$ and $\bm{R}$, which is widely used in rectangular low rank matrix recovery [12, 13, 14]. The convex feasible set ${\cal M}$ is defined as follows

where $\bm{L}_{\ell}$ is the $\ell$-th block of $\bm{L}$, $\mu$ and $\sigma$ be two absolute constants such that $\mu\geq\mu_{1}$ and $\sigma\geq\sigma_{1}$.

Inspired by [15], we design a projected gradient descent method for the problem (2.1), which is summarized in Algorithm 1.

The initialization involves two steps: (1) computes the best rank $r$ approximation of ${\cal H}({\cal A}^{\ast}(\bm{y}))$ via one step hard thresholding ${\cal P}_{r}(\cdot)$, where ${\cal A}^{\ast}$ is the adjoint of ${\cal A}$ ;(2) projects the low rank factors of best rank-$r$ approximated matrix onto the set ${\cal M}$. Given a matrix $\bm{M}=\begin{bmatrix}\bm{L}^{\mathsf{H}}&\bm{R}^{\mathsf{H}}\end{bmatrix}^{\mathsf{H}}$, the projection onto ${\cal M}$, denoted by $\begin{bmatrix}\widehat{\bm{L}}^{\mathsf{H}}&\widehat{\bm{R}}^{\mathsf{H}}\end{bmatrix}^{\mathsf{H}}$, has a closed form solution:

for $0\leq\ell\leq n_{1}-1$ and

for $0\leq j\leq n_{2}-1$. Let $\bm{M}_{t}$ be the current estimator. The algorithm updates $\bm{M}_{t}$ along gradient descent direction $-\nabla f(\bm{M}_{t})$ with step size $\eta$, followed by projection onto the set ${\cal M}$. The gradient of $f(\bm{M})$ is computed with respect to Wirtinger calculus given by $\nabla f=\begin{bmatrix}\nabla^{\mathsf{H}}_{\bm{L}}f&\nabla^{\mathsf{H}}_{\bm{R}}f\end{bmatrix}^{\mathsf{H}}$ where

To obtain the computational cost of $\nabla f$, we first introduce some notations. Let ${\cal H}_{v}$ be the Hankel operator which maps a vector $\bm{x}\in\mathbb{C}^{1\times n}$ into an $n_{1}\times n_{2}$ matrix,

where $x_{i}$ is the $i$-th entry of $\bm{x}$. The adjoint of ${\cal H}_{v}$, denoted by ${\cal H}_{v}^{\ast}$, is a linear mapping from $n_{1}\times n_{2}$ to $1\times n$. It is known that the computational complexity of both ${\cal H}_{v}^{\ast}(\bm{L}_{\ell}\bm{R}^{\mathsf{H}})$ and $\left({\cal H}_{v}(\bm{x})\right)\bm{R}$ is ${\cal O}(rn\log n)$ flops [16]. Moreover, the authors in [8] show that ${\cal H}(\bm{X})=\bm{P}\widetilde{{\cal H}}(\bm{X})$, where $\widetilde{{\cal H}}(\bm{X})$ is a matrix constructed by stacking all $\{{\cal H}_{v}(\bm{e}_{\ell}^{\mathsf{T}}\bm{X})\}_{\ell=1}^{s}$ on top of one another, and $\bm{P}$ is a permutation matrix. Therefore we can compute ${\cal H}^{\dagger}(\bm{L}\bm{R}^{\mathsf{H}})$ and ${\cal H}(\bm{X})\bm{R}$ by using ${\cal O}(srn\log n)$ flops. Thus the implementation of our algorithm is very efficient and the main computational complexity in each step is ${\cal O}(sr^{2}n+srn\log(n))$.