Signal Recovery on Incoherent Manifolds

Introduction

Estimation of an unknown signal from linear observations is a core problem in signal processing, statistics, and information theory. Particular energy has been invested in problem instances where the available information is limited and noisy and where the signals of interest possess a low-dimensional geometric structure. Indeed, focused efforts on certain instances of the linear inverse problem framework have spawned entire research subfields, encompassing both theoretical and algorithmic advances. Examples include signal separation and morphological component analysis ; sparse approximation and compressive sensing ; affine rank minimization ; and robust principal component analysis .

where ${\bf a}^{*}\in\mathcal{A},{\mathbf{b}}^{*}\in\mathcal{B}$. This expression for ${\bf x}$ contains $2N$ unknowns but only $N$ observations and hence is fundamentally ill-posed. Unless we make additional assumptions on the geometric structure of the component manifolds $\mathcal{A}$ and $\mathcal{B}$, a unique decomposition of ${\bf x}$ into its constituent signals $({\bf a}^{*},{\mathbf{b}}^{*})$ may not exist.

(Identifiability II) To complicate matters, in more general situations the linear operator $\mathbf{\Phi}$ in (1) might have fewer rows that columns, so that $M<N$. Thus, $\mathbf{\Phi}$ possesses a nontrivial nullspace. Indeed, we are particularly interested in cases where $M\ll N$, in which case the nullspace of $\Phi$ is extremely large relative to the ambient space. This further obscures the issue of identifiability of the ordered pair $({\bf a}^{*},{\mathbf{b}}^{*})$, given the available observations ${\bf z}$.

(Nonconvexity) Even if the above two identifiability issues were resolved, the manifolds $\mathcal{A},\mathcal{B}$ might be extremely nonconvex, or even non-differentiable. Thus, classical numerical methods, such as Newton’s method or steepest descent, cannot be successfully applied; neither can the litany of convex optimization methods that have been specially designed for linear inverse problems with certain types of signal priors .

In this paper, we propose a simple method to recover the component signals $({\bf a}^{*},{\mathbf{b}}^{*})$ from ${\bf z}$ in (1). We dub our method Successive Projections onto INcoherent manifolds (SPIN) (see Algorithm 1) Despite the highly nonconvex nature of the problem and the possibility of underdetermined measurements, SPIN provably recovers the signal components $({\bf a}^{*},{\mathbf{b}}^{*})$. For this to hold true, we will require that (i) the signal manifolds $\mathcal{A},\mathcal{B}$ are incoherent in the sense that the secants of $\mathcal{A}$ are almost orthogonal to the secants of $\mathcal{B}$; and (ii) the measurement operator $\mathbf{\Phi}$ satisfies a certain restricted isometry property (RIP) on the secants of the direct sum manifold $\mathcal{C}=\mathcal{A}\oplus\mathcal{B}$. We will formally define these conditions in Section 2. We prove the following theoretical statement below in Section 3.

Our proposed algorithm (SPIN) is iterative in nature. Each iteration consists of three steps: computation of the gradient of the error function $\psi({\bf a},{\mathbf{b}})=\frac{1}{2}\left\|{\bf z}-\mathbf{\Phi}({\bf a}+{\mathbf{b}})\right\|^{2}$, forming signal proxies for ${\bf a}$ and ${\mathbf{b}}$, and orthogonally projecting the proxies onto the manifolds $\mathcal{A}$ and $\mathcal{B}$. The projection operators onto the component manifolds play a crucial role in algorithm stability and performance; some manifolds admit stable, efficient projection operators while others do not. We discuss this in detail in Section 3. Additionally, we demonstrate that SPIN is stable to measurement noise (the quantity ${\bf e}$ in (1)) as well as numerical inaccuracies (such as finite precision arithmetic).

2 Prior Work

The core essence of our proposed approach has been extensively studied in a number of different contexts. Methods such as Projected Landweber iterations , iterative hard thresholding (IHT) , and singular value projection (SVP) are all instances of the same basic framework. SPIN subsumes and generalizes these methods. In particular, SPIN is an iterative projected gradient method with the same basic approach as two recent signal recovery algorithms — Gradient Descent with Sparsification (GraDeS) , and Manifold Iterative Pursuit (MIP) . We generalize these approaches to situations where the signal of interest is a linear mixture of signals arising from a pair of nonlinear manifolds. Due to the particular structure of our setting, SPIN consists of two projection steps (instead of one), and the analysis is more involved (see Section 4). We also explore the interplay between the geometric structure of the component manifolds, the linear measurement operator, and the stability of the recovery algorithm.

SPIN exhibits a strong geometric convergence rate comparable to many state-of-the-art first-order methods , despite the nonlinear and nonconvex nature of the reconstruction problem. We duly note that, for the case of certain special manifolds, sophisticated higher-order recovery methods with stronger stability guarantees have been proposed (e.g., approximate message passing (AMP) for sparse signal recovery and augmented Lagrangian multiplier (ALM) methods for low-rank matrix recovery ); see also . However, an appealing feature of SPIN is its conceptual simplicity plus its ability to generalize to mixtures of arbitrary nonlinear manifolds, provided these manifolds satisfy certain geometric properties, as detailed in Section 2.

3 Setup

Geometric Assumptions

In linear inverse problems such as sparse signal approximation and compressive sensing, the assumption of incoherence between linear subspaces, bases, or dictionary elements is common. We introduce a nonlinear generalization of this concept.

The secant manifold $\mathcal{S}(\mathcal{A})$ is the family of unit vectors ${\bf u}$ generated by all pairs ${\bf a},{\bf a}^{\prime}$ in $\mathcal{A}$.

where $\mathcal{S}(\mathcal{A}),\mathcal{S}(\mathcal{B})$ are the secant manifolds of $\mathcal{A},\mathcal{B}$ respectively. Then, $\mathcal{A}$ and $\mathcal{B}$ are called $\epsilon$-incoherent manifolds.

Informally, any point on the secant manifold $\mathcal{S}(\mathcal{A})$ represents a direction that aligns with a difference vector of $\mathcal{A}$, while the incoherence parameter $\epsilon$ controls the extent of “perpendicularity” between the manifolds $\mathcal{A}$ and $\mathcal{B}$. We define $\epsilon$ in terms of a supremum over sets $\mathcal{S}(\mathcal{A}),\mathcal{S}(\mathcal{B})$. Therefore, a small value of $\epsilon$ implies that each (normalized) secant of $\mathcal{A}$ is approximately orthogonal to all secants of $\mathcal{B}$. By definition, the quantity $\epsilon$ is always non-negative; further, $\epsilon\leq 1$, due to the Cauchy-Schwartz inequality.

We prove that any signal ${\bf x}$ belonging to the direct sum $\mathcal{A}\oplus\mathcal{B}$ can be uniquely decomposed into its constituent signals when the upper bound on $\epsilon$ holds with strict inequality.

Suppose that $\mathcal{A},\mathcal{B}$ are $\epsilon$-incoherent with $0<\epsilon<1$. Consider ${\bf x}={\bf a}+{\mathbf{b}}={\bf a}^{\prime}+{\mathbf{b}}^{\prime}$, where ${\bf a},{\bf a}^{\prime}\in\mathcal{A}$ and ${\mathbf{b}},{\mathbf{b}}^{\prime}\in\mathcal{B}$. Then, ${\bf a}={\bf a}^{\prime},{\mathbf{b}}={\mathbf{b}}^{\prime}$.

Proof. It is clear that $\left\|{\bf a}+{\mathbf{b}}-({\bf a}^{\prime}+{\mathbf{b}}^{\prime})\right\|^{2}=0$, i.e.,

However, due to the manifold incoherence assumption, the (unnormalized) secants ${\bf a}-{\bf a}^{\prime}$, ${\mathbf{b}}-{\mathbf{b}}^{\prime}$ obey the relation:

where the last inequality follows from the relation between arithmetic and geometric means (henceforth referred to as the AM-GM inequality). Therefore, we have that

for $\epsilon<1$, which is impossible unless ${\bf a}={\bf a}^{\prime},{\mathbf{b}}={\mathbf{b}}^{\prime}$. $\Box$

We can also prove the following relation between secants and direct sums of signals lying on incoherent manifolds.

Suppose that $\mathcal{A},\mathcal{B}$ are $\epsilon$-incoherent with $0<\epsilon<1$. Consider ${\bf x}_{1}={\bf a}_{1}+{\mathbf{b}}_{1},{\bf x}_{2}={\bf a}_{2}+{\mathbf{b}}_{2},$ where ${\bf a}_{1},{\bf a}_{2}\in\mathcal{A}$ and ${\mathbf{b}}_{1},{\mathbf{b}}_{2}\in\mathcal{B}$. Then

Rearranging terms, we obtain the desired result. $\Box$

2 Restricted isometry

The notion of restricted isometry (and its generalizations) is an important component in the analysis of many algorithms in sparse approximation, compressive sensing, and low-rank matrix recovery . While the RIP has traditionally been studied in the context of sparse signal models, (5) generalizes this notion to arbitrary nonlinear manifolds. The restricted isometry condition is of particular interest when the range space of the matrix $\mathbf{\Phi}$ is low-dimensional. A key result states that, under certain upper bounds on the curvature of the manifold $\mathcal{C}$, there exist probabilistic constructions of matrices $\mathbf{\Phi}$ that satisfy the RIP on $\mathcal{C}$ such that the number of rows of $\mathbf{\Phi}$ is proportional to the intrinsic dimension of $\mathcal{C}$, rather than the ambient dimension $N$ of the signal space. We will discuss this further in Section 5.

3 Projections onto manifolds

The projection operator $\mathcal{P}_{\mathcal{A}}(\cdot)$ plays a crucial role in the development of our proposed signal recovery algorithm in Section 3. Note that in a number of applications, $\mathcal{P}_{\mathcal{A}}(\cdot)$ may be quite difficult to compute exactly. The reasons for this might be intrinsic to the application (such as the nonconvex, non-differentiable structure of $\mathcal{A}$), or might be due to extrinsic constraints (such as finite-precision arithmetic). Therefore, following the lead of , we also define a $\gamma$-approximate projection operator onto $\mathcal{A}$:

so that $\mathcal{P}_{\mathcal{A}}^{\gamma}({\bf x})$ yields a vector ${\bf x}^{\prime}\in\mathcal{A}$ that approximately minimizes the squared distance from ${\bf x}$ to $\mathcal{A}$. Again, $\mathcal{P}_{\mathcal{A}}^{\gamma}({\bf x})$ need not be uniquely defined for a particular input signal ${\bf x}$.

The SPIN Algorithm

We now describe an algorithm to solve the linear inverse problem (1). Our proposed algorithm, Successive Projections onto INcoherent manifolds (SPIN), can be viewed as a generalization of several first-order methods for signal recovery for a variety of different models . SPIN is described in pseudocode form in Algorithm 1.

The key innovation in SPIN is that we formulate two proxy vectors for the signal components $\widetilde{{\bf a}}_{k}$ and $\widetilde{{\mathbf{b}}}_{k}$ and project these onto the corresponding manifolds $\mathcal{A}$ and $\mathcal{B}$.

We demonstrate that SPIN possesses strong uniform recovery guarantees comparable to existing state-of-the-art algorithms for sparse approximation and compressive sensing, while encompassing a very broad range of nonlinear signal models. The following theoretical result describes the performance of SPIN for signal recovery.

then SPIN (Algorithm 1) with step size $\eta=1/(1+\delta)$ with exact projections $\mathcal{P}_{\mathcal{A}},\mathcal{P}_{\mathcal{B}}$ outputs ${\bf a}_{T}\in\mathcal{A}$ and ${\mathbf{b}}_{T}\in\mathcal{B}$, such that $\left\|{\bf z}-\mathbf{\Phi}({\bf a}_{T}+{\mathbf{b}}_{T})\right\|^{2}\leq\beta\left\|{\bf e}\right\|^{2}+\nu$ in no more than $T=\lceil\frac{1}{\log(1/\alpha)}\log{\frac{\left\|{\bf z}\right\|^{2}}{2\nu}}\rceil$ iterations for any $\nu>0$.

Here, $\alpha<1$ and $\beta$ are moderately-sized positive constants that depend only on $\delta$ and $\epsilon$; we derive explicit expressions for $\alpha$ and $\beta$ in Section 4. For example, when $\epsilon=0.05,~{}\delta=0.5$, we obtain $\alpha\approx 0.812,~{}\beta\approx 5.404$.

For the special case when there is no measurement noise (i.e., ${\bf e}=0$), Theorem 2 states that, after a finite number of iterations, SPIN outputs signal component estimates $(\widehat{{\bf a}},\widehat{{\mathbf{b}}})$ such that $\left\|{\bf z}-\mathbf{\Phi}(\widehat{{\bf a}}+\widehat{{\mathbf{b}}})\right\|<\nu$ for any desired precision parameter $\nu$. From the restricted isometry assumption on $\mathbf{\Phi}$ and Lemma 1, we immediately obtain Theorem 1. Since we can set $\nu$ to an arbitrarily small value, we have that the SPIN estimate $(\widehat{{\bf a}},\widehat{{\mathbf{b}}})$ converges to the true signal pair $({\bf a}^{*},{\mathbf{b}}^{*})$. Exact convergence of the algorithm might potentially take a very large number of iterations, but convergence to any desired positive precision constant $\beta$ takes only a finite number of iterations. For the rest of the paper, we will informally denote signal “recovery” to imply convergence to a sufficiently fine precision.

SPIN assumes the availability of the exact projection operators $\mathcal{P}_{\mathcal{A}},\mathcal{P}_{\mathcal{B}}$. In certain cases, it might be feasible to numerically compute only $\gamma$-approximate projections, as in (7). In this case, the bound on the norm of the error ${\bf z}-\mathbf{\Phi}({\bf a}_{T}+{\mathbf{b}}_{T})$ is only guaranteed to be upper bounded by a positive multiple of the approximation parameter $\gamma$. The following theoretical guarantee (with a near-identical proof mechanism as Theorem 2) captures this behavior.

Under the same suppositions as Theorem 2, SPIN (Algorithm 1) with $\gamma$-approximate projections and step size $\eta=1/(1+\delta)$ outputs ${\bf a}_{T}\in\mathcal{A}$ and ${\mathbf{b}}_{T}\in B$ such that $\left\|{\bf z}-\mathbf{\Phi}({\bf a}_{T}+{\mathbf{b}}_{T})\right\|^{2}\leq\beta\left\|{\bf e}\right\|^{2}+\frac{1+\delta}{1-\alpha}\gamma+\nu$, in no more than $T=\lceil\frac{1}{\log(1/\alpha)}\log{\frac{\left\|{\bf z}\right\|^{2}}{2\nu}}\rceil$ iterations.

We note some implications of Theorem 2. First, suppose that $\Phi$ is the identity operator, i.e., we have full measurements of the signal ${\bf x}^{*}={\bf a}^{*}+{\mathbf{b}}^{*}$. Then, $\delta=0$ and the lower bound on the restricted isometry constant holds with equality. However, we still require that $\epsilon<1/11$ for guaranteed recovery using SPIN. We will discuss this condition further in Section 5.

Second, suppose that the one of the component manifolds is the trivial (zero) manifold; then, we have that $\epsilon=0$. In this case, SPIN reduces to the Manifold Iterative Pursuit (MIP) algorithm for recovering signals from a single manifold . Moreover, the condition on $\delta$ reduces to $0\leq\delta<1/3$, which exactly matches the condition required for guaranteed recovery using MIP.

Lastly, the condition (8) in Theorem 2 automatically implies that $\epsilon<1/11$. This represents a mild tightening of the condition on $\epsilon$ required for a unique decomposition (Lemma 1), even with full measurements (i.e., when $\mathbf{\Phi}$ is the identity operator or, more generally, when $\delta=0$).

Analysis

It is clear that $\psi({\bf a}^{*},{\mathbf{b}}^{*})=\frac{1}{2}\left\|{\bf e}\right\|^{2}$. The following lemma bounds the error of the estimated signals output by SPIN at the $(k+1)$-st iteration in terms of the error incurred at the $k$-th iteration, and the norm of the measurement error.

Define $({\bf a}_{k},{\mathbf{b}}_{k})$ as the intermediate estimates obtained by SPIN at the $k$-th iteration. Let $\delta,\epsilon$ be as defined in Theorem 2. Then,

Proof. Fix a current estimate of the signal components $({\bf a}_{k},{\mathbf{b}}_{k})$ at iteration $k$. Then, for any other pair of signals $({\bf a},{\mathbf{b}})\in\mathcal{A}\times\mathcal{B}$, we have

where ${\bf x}_{k}\triangleq{\bf a}_{k}+{\mathbf{b}}_{k},~{}{\bf x}\triangleq{\bf a}+{\mathbf{b}}$. Since $\mathbf{\Phi}$ is a linear operator, we can take the adjoint within the inner product to obtain

The last inequality occurs due to the RIP of $\mathbf{\Phi}$ applied to the secant vector ${\bf x}-{\bf x}_{k}\in\mathcal{S}(\mathcal{C})$. To the right hand side of (11), we further add and subtract $\frac{1}{2(1+\delta)}\left\|\mathbf{\Phi}^{T}({\bf z}-\mathbf{\Phi}{\bf x}_{k})\right\|^{2}$ to complete the square:

Define ${\mathbf{g}}_{k}\triangleq\frac{1}{1+\delta}\mathbf{\Phi}^{T}({\bf z}-\mathbf{\Phi}({\bf a}_{k}+{\mathbf{b}}_{k}))$. Then,

Next, define the function $\zeta$ on $\mathcal{A}\times\mathcal{B}$ as $\zeta({\bf a},{\mathbf{b}})\triangleq\left\|{\bf a}+{\mathbf{b}}-({\bf a}_{k}+{\mathbf{b}}_{k}+{\mathbf{g}}_{k})\right\|^{2}$. Then, we have

But, as specified in Algorithm 1, ${\bf a}_{k+1}=\mathcal{P}_{\mathcal{A}}({\bf a}_{k}+{\mathbf{g}}_{k})$, and hence $\left\|{\bf a}_{k+1}-({\bf a}_{k}+{\mathbf{g}}_{k})\right\|\leq\left\|{\bf a}-({\bf a}_{k}+{\mathbf{g}}_{k})\right\|$ for any ${\bf a}\in\mathcal{A}$. An analogous relation can be formed between ${\mathbf{b}}_{k+1}$ and ${\mathbf{b}}^{*}$. Hence, we have

Substituting for $({\bf a}_{k+1},{\mathbf{b}}_{k+1})$, we obtain

The last term on the right hand side equals zero, and so we obtain

Combining this inequality with (12), we obtain the series of inequalities

Again, ${\bf x}^{*}-{\bf x}_{k}$ is a secant on the direct sum manifold $\mathcal{C}$. By the RIP property of $\mathbf{\Phi}$, we have

By definition, we have that $\psi({\bf a}_{k},{\mathbf{b}}_{k})=\frac{1}{2}\left\|{\bf z}-\Phi{\bf x}_{k}\right\|^{2}$. Further, we can substitute ${\bf a}={\bf a}^{*},{\mathbf{b}}={\mathbf{b}}^{*}$ in (10) to obtain

via the same technique used to obtain (16). Similarly,

To ensure that the value of $\psi({\bf a}_{k},{\mathbf{b}}_{k})$ does not diverge, the leading coefficient $\alpha$ must be smaller than 1, i.e.,

Rearranging, we obtain the upper bound on $\delta$ as in (8):

By choosing $\beta=\frac{C}{1-\alpha}$, and $k\geq T$ such that $T=\lceil\frac{1}{\log(1/\alpha)}\log{\frac{\left\|{\bf z}\right\|^{2}}{2\nu}}\rceil$, the result follows. $\Box$

The proof mechanism of Theorem 3 follows a near-identical procedure as in Lemma 3 and we omit the details for brevity. Also, we observe that Theorem 2 represents merely a sufficient condition for signal recovery; the constants in (8) could likely be improved, but we will not pursue that direction in this paper.

Applications

The two-manifold signal model described in this paper is applicable to a wide variety of problems that have attracted considerable interest in the literature over the last several years. We discuss a few representative instances and show how SPIN can be utilized for efficient signal recovery in each of these instances. We also present several numerical experiments that indicate the kind of gains that SPIN can offer in practice.

The problem is to recover $({\bf a}^{*},{\mathbf{b}}^{*})$ given ${\bf x}^{*}$. This problem has been studied in many different forms in the literature, and several algorithms have been proposed in order to solve it efficiently . See for an in-depth study of the various state-of-the-art methods. All these methods assume a certain notion of incoherence between the two bases, most commonly referred to as the mutual coherence $\mu$, which is defined as

We establish the following simple relation between $\mu$ and the manifold incoherence between $\mathcal{A}$ and $\mathcal{B}$.

Let $\mathcal{A}$ be the set of all $K_{1}$-sparse signals in $\mathbf{\Lambda}$, and $\mathcal{B}$ be the set of all $K_{2}$-sparse signals in $\mathbf{\Lambda}^{\prime}$. Let $\epsilon$ denote the manifold incoherence between $\mathcal{A}$ and $\mathcal{B}$. Then,

where the last relation follows from the triangle inequality. We can further bound the right hand side of (19). We have

since ${\bf u}$ is a unit vector. Similarly, $\sum_{i=1}^{2K_{2}}|b_{i}|\leq\sqrt{2K_{2}}$. Inserting these upper bounds in (19), we have

The lemma follows by considering the supremum over all vectors ${\bf u}\in\mathcal{S}(\mathcal{A}),{\bf u}^{\prime}\in\mathcal{S}(\mathcal{B})$. $\Box$

We show how SPIN can be used to solve the linear inverse problem of recovering $({\bf a}^{*},{\mathbf{b}}^{*})$ from ${\bf x}^{*}$. The restricted isometry assumption is not relevant in this case, since we assume that we have full measurements of the signal; therefore $\delta=0$. An upper bound for the manifold incoherence parameter $\epsilon$ is specified in Lemma 4. The (exact) projection operators $\mathcal{P}_{\mathcal{A}},\mathcal{P}_{\mathcal{B}}$ can be easily implemented; we simply perform a coefficient expansion in the corresponding orthornormal basis and retain the coefficients of largest magnitude. Mixing these ingredients together, we can guarantee that, given any signal ${\bf x}^{*}$, SPIN will return the true components ${\bf a}^{*},{\mathbf{b}}^{*}$. This guarantee is summarized in the following result.

Let ${\bf x}^{*}={\bf a}^{*}+{\mathbf{b}}^{*}$, where ${\bf a}^{*}$ is $K_{1}$-sparse in $\mathbf{\Lambda}$ and ${\mathbf{b}}^{*}$ is $K_{2}$-sparse in $\mathbf{\Lambda}^{\prime}$. Let $\mu$ denote the mutual coherence between $\mathbf{\Lambda}$ and $\mathbf{\Lambda}^{\prime}$. Then, SPIN exactly recovers (${\bf a}^{*},{\mathbf{b}}^{*}$) from ${\bf x}^{*}$ provided

Proof. If (20) holds, then from Lemma 4 we know that the manifold incoherence $\epsilon$ between $\mathcal{A}$ and $\mathcal{B}$ is smaller than 1/11. But this is exactly the condition required for guaranteed convergence of SPIN to the true signal components $({\bf a}^{*},{\mathbf{b}}^{*})$. $\Box$

Therefore, SPIN yields a recovery guarantee that is off the best possible method by a factor of 10. Once again, it is possible that the constant $1/11$ in Corollary 1 can be tightened by a more careful analysis of SPIN specialized to the case when the component signal manifolds correspond to a pair of incoherent bases, but we will not pursue this direction here. It is also possible to generalize SPIN to the case where the sparsifying dictionary comprises a union of more than two orthonormal bases ; see Section 6 for a short discussion.

2 Articulation manifolds

for some constant $C$ that depends only on the smoothness and volume of the manifold $\mathcal{M}$. Therefore, the dimension of the range space of $\mathbf{\Phi}$ is proportional to the number of degrees of freedom $K$, but is only logarithmic in the ambient dimension $N$. Moreover, given such a measurement matrix $\mathbf{\Phi}$ with isometry constant $\delta<1/3$ and a projection operator $\mathcal{P}_{\mathcal{M}}(\cdot)$ onto $\mathcal{M}$, any signal ${\bf x}\in\mathcal{M}$ can be reconstructed from its compressive measurements ${\bf y}=\mathbf{\Phi}{\bf x}$ using Manifold Iterative Pursuit (MIP) .

We generalize this setting to the case where the unknown signal of interest arises as a mixture of signals from two manifolds $\mathcal{A}$ and $\mathcal{B}$. For instance, suppose we are interested in the space of images, where $\mathcal{A}$ and $\mathcal{B}$ comprise of translations of fixed template images $f({\bf t})$ and $g({\bf t})$, where ${\bf t}$ denotes the 2D domain over which the image is defined. Then, the signal of interest is an image of the form

where $\theta_{1}$ and $\theta_{2}$ denote the unknown translation parameters. The problem is to recover $({\bf a}^{*},{\mathbf{b}}^{*})$, or equivalently $(\theta_{1},\theta_{2})$, given compressive measurements ${\bf z}=\mathbf{\Phi}({\bf a}^{*}+{\mathbf{b}}^{*})$.

We demonstrate that SPIN offers an easy, efficient technique to recover the component images. This example also demonstrates that SPIN is robust to practical considerations such as noise. Figure 1 displays the results of SPIN recovery of a $64\times 64$ image from very limited measurements. The unknown image consists of the linear sum of arbitrary translations of template images $f({\bf t})$ and $g({\bf t})$, that are smoothed binary images on a black background of a white disk and a white square, respectively. Further, the image has been contaminated with significant Gaussian noise (SNR = 14dB) prior to measurement (Fig. 1(a)). From Figs. 1(b) and 1(c), we observe that SPIN is able to perfectly recover the original component signals from merely $M=50$ random linear measurements.

For guaranteed SPIN convergence, we require that the manifolds $\mathcal{A},\mathcal{B}$ are incoherent. Informally, the condition of incoherence on the secants of $\mathcal{A}$ and $B$ is always valid when the template images $f({\bf t}),g({\bf t})$ are “sufficiently” distinct. This intuition is made precise using the bound in (3). More generally, we can state the following theoretical guarantee for SPIN performance in the case of general higher-dimensional manifolds.

Here, $C_{\mathcal{A}},C_{\mathcal{B}}$ are constants that depend only on certain intrinsic geometric parameters (such as the volume) of $\mathcal{A},\mathcal{B}$ respectively.

Proof. It is easy to see that if the matrix $\mathbf{\Phi}$ satisfies the RIP on the secants of the direct sum $\mathcal{C}=\mathcal{A}\oplus\mathcal{B}$, then SPIN recovery follows from Theorem 2. We show that a randomized construction of $\mathbf{\Phi}$ with number of rows specified by (21) satisfies the RIP on $\mathcal{C}$ with high probability. Essentially, our proof combines the techniques used in Section 3.2 of with Lemma 1 of .

where $C_{\mathcal{A}}$ is a constant that depends only on the intrinsic geometry of $\mathcal{A}$. However, in our setting we are interested in the direct sum of manifolds; correspondingly, we can construct finite sets $\mathcal{R}_{\mathcal{A}}$ and $\mathcal{R}_{\mathcal{B}}$ and apply Lemma 1 of , that specifies a lower bound on the number of measurements required to preserve the norms of linear sums of finite point sets:

where $K,K^{\prime}$ are the dimensions of $\mathcal{A},\mathcal{B}$ respectively. Corollary 2 follows. $\Box$

3 Signals in impulsive noise

In some situations, the signal of interest ${\bf x}$ might be corrupted with impulsive noise (or shot noise) prior to signal acquisition via linear measurements. For example, consider Fig. 2(a), where the Gaussian pulse is the signal of interest, and the spikes indicate the undesirable noise. In this case, the linear observations are more accurately modeled as:

and ${\bf n}$ is a $K^{\prime}$-sparse signal in the canonical basis. Therefore, SPIN can be used to recover ${\bf x}$ from ${\bf z}$, provided that the manifold $\mathcal{M}$ is incoherent with the set of sparse signals $\Sigma_{K^{\prime}}$ and $\mathbf{\Phi}$ satisfies the RIP on the direct sum $\mathcal{M}+\Sigma_{K^{\prime}}$.

We apply SPIN to recover ${\bf x}$ from ${\bf z}$. The projection operator $\mathcal{P}_{\mathcal{M}}(\cdot)$ consists of a matched filter with the template pulse ${\mathbf{g}}_{0}$, while the projection operator $\mathcal{P}_{\Sigma_{K^{\prime}}}(\cdot)$ simply returns the best $K^{\prime}$-term approximation in the canonical basis. Assuming that we have knowledge of the number of nonzeros in the noise vector ${\bf n}$, we can use SPIN to reconstruct both ${\bf x}$ and ${\bf n}$. We observe from Fig. 2(b) that SPIN recovers the true signal ${\bf x}$ with near-perfect accuracy. Further, this recovery is possible with only a small number $M=150$ linear measurements of ${\bf x}$, which constitutes but a fraction of the ambient dimension of the signal space.

Figure 3 plots the number of measurements $M$ vs. the signal reconstruction error (normalized relative to the signal energy and plotted in dB). We observe that, by increasing $M$, SPIN can tolerate an increased number $K^{\prime}$ of nuisance spikes. Further, by Corollary 2, we observe that this relationship between $M$ and $K^{\prime}$ is in fact linear. This result can be extended to any situation where the signals of interest obey a “hybrid” model that is a mixture of a nonlinear manifold and the set of sparse signals.

Discussion

We have proposed and rigorously analyzed an algorithm, which we dub Successive Projections onto INcoherent Manifolds (SPIN), for the recovery of a pair of signals given a small number of measurements of their linear sum. For SPIN to guarantee signal recovery, we require two main geometric criteria to hold: (i) the component signals should arise from two disjoint manifolds that are in a specific sense incoherent, and (ii) the linear measurement operator should satisfy a restricted isometry criterion on the secants of the direct sum of the two manifolds. The computational efficiency of SPIN is determined by the tractability of the projection operators onto either component manifold. We have presented indicative numerical experiments demonstrating the utility of SPIN, but defer a thorough experimental study of SPIN to future work.

Practical considerations. SPIN is an iterative gradient projection algorithm and requires as input parameters the number of iterations $T$ and the gradient step size $\eta$. The iteration count $T$ can be chosen using one of many commonly-used stopping criteria. For example, convergence can be declared if the norm of the error $\psi({\bf a}_{k},{\mathbf{b}}_{k})$ at the $(k+1)$-th time step does not differ significantly from the error at $k$-th time step. The choice of optimal step size $\eta$ is more delicate. Theorem 2 relates the step size to the restricted isometry constant $\delta$ of $\mathbf{\Phi}$, but this constant is not easy to calculate. In our preliminary findings, a step size in the range $0.5\leq\eta\leq 0.7$ consistently gave good results. See for a discussion on the choice of step size for hard thresholding methods.

In practical scenarios, the signal of interest rarely belongs exactly to a low-dimensional submanifold $\mathcal{M}$ of the ambient space, but is only well-approximated by $\mathcal{M}$. Interestingly, in such situations the effect of this mismatch can be studied using the concept of $\gamma$-approximate projections (7). Theorem 3 rigorously demonstrates that SPIN is robust to such approximations. Further, our main result (Theorem 2) indicates that SPIN is stable with respect to inaccurate measurements, owing to the fact that the reconstruction error is bounded by a constant times the norm of the measurement noise vector ${\bf e}$.

More than two manifolds. For clarity and brevity, we have focused our attention on signals belonging to the direct sum of two signal manifolds. However, SPIN (and its accompanying proof mechanism) can be conceptually extended to sums of any $Q$ manifolds. In such a scenario, the conditions of convergence of SPIN would require that the component manifolds are $Q$-wise incoherent, and the measurement operator $\mathbf{\Phi}$ satisfies a restricted isometry on the $Q$-wise direct sum of the component manifolds.

Connections to matrix recovery. An intriguing open question is whether SPIN (or a similar first-order projected gradient algorithm) is applicable to situations where either of the component manifolds is the set of low-rank matrices. The problem of reconstructing, from affine measurements, matrices that are a sum of low-rank and sparse matrices has attracted significant attention in the recent literature . The key stumbling block is that the manifold of low-rank matrices is not incoherent with the manifold of sparse matrices; indeed, the two manifolds share a nontrivial intersection (i.e., there exist low rank matrices that are also sparse, and vice versa). Phenomena such as these make the analysis of SPIN (or similar algorithms) quite challenging, and it may be that higher-order techniques will be needed for signal recovery.

References