Signal Reconstruction from Modulo Observations

TL;DR

This paper introduces a novel algorithm for reconstructing signals from under-determined modulo observations, inspired by phase retrieval techniques, demonstrating perfect recovery with enough measurements and superior performance over existing methods.

Contribution

The paper presents a new approach for signal reconstruction from modulo measurements limited to two periods, extending phase retrieval algorithms to this challenging setting.

Findings

Algorithm achieves perfect signal recovery with sufficient measurements.

Proposed method outperforms existing algorithms in experiments.

Validated on synthetic and real data showing superior performance.

Abstract

We consider the problem of reconstructing a signal from under-determined modulo observations (or measurements). This observation model is inspired by a (relatively) less well-known imaging mechanism called modulo imaging, which can be used to extend the dynamic range of imaging systems; variations of this model have also been studied under the category of phase unwrapping. Signal reconstruction in the under-determined regime with modulo observations is a challenging ill-posed problem, and existing reconstruction methods cannot be used directly. In this paper, we propose a novel approach to solving the inverse problem limited to two modulo periods, inspired by recent advances in algorithms for phase retrieval under sparsity constraints. We show that given a sufficient number of measurements, our algorithm perfectly recovers the underlying signal and provides improved performance over other existing algorithms. We also provide experiments validating our approach on both synthetic and real data to depict its superior performance.

Tables

Table 1. TABLE I: Comparison of MoRAM with existing modulo recovery methods.

	Unlimited Sampling [Bhandari]	OLS Method [Cucuringu2017]	multishot UHDR [ICCP15_Zhao]	MoRAM (our approach)
Assumption on structure of signal	Bandlimited	Bandlimited	No assumptions	Sparsity
Sampling scheme	uniform grid	uniform grid	(carefully chosen) linear measurements	random linear measurements
Sample complexity	oversampled, $𝒪\left(n\right)$	–	oversampled, $𝒪\left(n\right)$	undersampled, $𝒪\left(s\log \left(n\right)\right)$
Provides sample complexity bounds?	Yes	–	No	Yes
Leverages Sparsity?	No	No	No	Yes
(Theoretical) bound on dynamic range	Unbounded	Unbounded	Unbounded	$2R$

Equations

$$y_i=\mathrm{mod}(⟨{\mathbf{a}}_{\mathbf{i}}\cdot {\mathbf{x}}^{\ast }⟩,R)i=\{1,2,...,m\},$$

$$y_i=|⟨{\mathbf{a}}_{\mathbf{i}},{\mathbf{x}}^{\ast }⟩|,i=1,2,\dots ,m,$$

$$f(t):=t+\left(\frac{1-\mathrm{sgn}(t)}{2}\right)R.$$

$$y_i=⟨{\mathbf{a}}_{\mathbf{i}}\cdot {\mathbf{x}}^{\ast }⟩+\left(\frac{1-\mathrm{sgn}(⟨{\mathbf{a}}_{\mathbf{i}}\cdot {\mathbf{x}}^{\ast }⟩)}{2}\right)R,i=\{1,..,m\}.$$

$$p_i^{\ast }=\frac{1-\mathrm{sgn}(⟨{\mathbf{a}}_{\mathbf{i}}\cdot {\mathbf{x}}^{\ast }⟩)}{2},i=\{1,..,m\}.$$

$$y_i=⟨{\mathbf{a}}_{\mathbf{i}}\cdot {\mathbf{x}}^{\ast }⟩+p_i^{\ast }R=y_{c_i}+p_i^{\ast }R.$$

$$p_i^{init}=\{\begin{array}{ll}0,& \text{if }0\le y_i<R/2\\ 1,& \text{if }R/2\le y_i\le R\end{array}$$

$${\mathbf{y}}_{\mathbf{c}}^{\mathbf{init}}=\mathbf{y}+{\mathbf{p}}^{\mathbf{init}}\mathbf{R}.$$

$${\mathbf{x}}^0=H_s\left(\frac{1}{N}\sum _{i=1}^N{y}_{c,i}^{init}{a}_i\right),$$

$${\mathbf{p}}^{\mathbf{t}}=\frac{1-\mathrm{sgn}(⟨\mathbf{A}\cdot {\mathbf{x}}^{\mathbf{t}}⟩)}{2}.$$

$${\mathbf{y}}_{\mathbf{c}}^{\mathbf{t}}=⟨\mathbf{A}{\mathbf{x}}^{\mathbf{t}+1}⟩=\mathbf{y}-{\mathbf{p}}^{\mathbf{t}}R,$$

$${\mathbf{x}}^{\mathbf{t}+1}=\underset{\mathbf{x}\in {\mathcal{M}}_s}{\mathrm{argmin}}\Vert \mathbf{Ax}-{\mathbf{y}}_{\mathbf{c}}^{\mathbf{t}}{\Vert }_2^2,$$

$${\mathbf{x}}^{\mathbf{t}+1}=\underset{[\mathbf{x}\mathbf{d}{]}^{\top }\in {\mathcal{M}}_{s+s_{dt}}}{\mathrm{argmin}}\Vert \left[\begin{array}{cc}\mathbf{A}& \mathbf{I}\end{array}\right]\left[\begin{array}{c}\mathbf{x}\\ \mathbf{d}\end{array}\right]-{\mathbf{y}}_{\mathbf{c}}^{\mathbf{t}}{\Vert }_2^2,$$

$$⟹{\mathbf{x}}^{\mathbf{t}+1}=JP(\frac{1}{\sqrt{m}}\left[\begin{array}{cc}\mathbf{A}& \mathbf{I}\end{array}\right],\frac{1}{\sqrt{m}}{\mathbf{y}}_{\mathbf{c}}^{\mathbf{t}},[{\mathbf{x}}^{\mathbf{t}}{\mathbf{p}}^{\mathbf{t}}{]}^{\top }).$$

$${\mathbf{p}}^0=\frac{1-\mathrm{sgn}(⟨\mathbf{A}\cdot {\mathbf{x}}^0⟩)}{2}.$$

$${\mathbf{y}}_{\mathbf{c}}^0=⟨\mathbf{A}{\mathbf{x}}^0⟩=\mathbf{y}-{\mathbf{p}}^{0}R,$$

$$d_S(\mathbf{x},\mathbf{y})-ϵ\le d_H(F(\mathbf{x}),F(\mathbf{y}))\le d_S(\mathbf{x},\mathbf{y})+ϵ;$$

$$F(\mathbf{x}):=\mathrm{sgn}(\mathbf{Ax});$$

$$d_H(F({\mathbf{x}}^{\ast }),F({\mathbf{x}}^0))\le d_S({\mathbf{x}}^{\ast },{\mathbf{x}}^0)+ϵ.$$

$$2d_s(\mathbf{p},\mathbf{q})\le \Vert \mathbf{p}-\mathbf{q}{\Vert }_2\le \pi d_s(\mathbf{p},\mathbf{q})$$

$$d_H(\mathrm{sgn}(\mathbf{A}{\mathbf{x}}^{\ast }),\mathrm{sgn}(\mathbf{A}{\mathbf{x}}^0))$$

Full text

I Introduction

I-A Motivation

The problem of reconstructing a signal (or image) from (possibly) nonlinear observations is widely encountered in standard signal acquisition and imaging systems. Our focus in this paper is the problem of signal reconstruction from modulo measurements, where the modulo operation with respect to a positive real valued parameter $R$ returns the (fractional) remainder after division by $R$. See Fig. 1 for an illustration.

Formally, we consider a high dimensional signal (or image) $\mathbf{x}^{*}\in{\mathbb{R}}^{n}$. We are given modulo measurements of $\mathbf{x^{*}}$, that is, for each measurement vector $\mathbf{a_{i}}\in{\mathbb{R}}^{n}$, we observe:

[TABLE]

The task is to recover $\mathbf{x^{*}}$ using the modulo measurements $\mathbf{y}$ and knowledge the measurement matrix $\mathbf{A}=\left[\mathbf{a_{1}~{}a_{2}~{}...~{}a_{m}}\right]^{\top}$.

This specific form of signal recovery is gaining rapid interest in recent times. Recently, the use of a novel imaging sensor that wraps the data in a periodical manner has been shown to overcome certain hardware limitations of typical imaging systems [Bhandari, ICCP15_Zhao, Shah, Cucuringu2017]. Many image acquisition systems suffer from the problem of limited dynamic range; however, real-world signals can contain a large range of intensity levels, and if tuned incorrectly, most intensity levels can lie in the saturation region of the sensors, causing loss of information through signal clipping. The problem gets amplified in the case of multiplexed linear imaging systems (such as compressive cameras or coded aperture systems), where required dynamic range is very high because of the fact that each linear measurement is a weighted aggregation of the original image intensity values.

The standard solution to this issue is to improve sensor dynamic range via enhanced hardware; this, of course, can be expensive. An intriguing alternative is to deploy special digital modulo sensors [rheejoo, kavusi2004quantitative, sasagawa2016implantable, yamaguchi2016implantable]. As the name suggests, such a sensor wraps each signal measurement around a scalar parameter $R$ that reflects the dynamic range. However, this also makes the forward model (1) highly nonlinear and the reconstruction problem highly ill-posed. The approach of [Bhandari, ICCP15_Zhao] resolves this problem by assuming overcomplete observations, meaning that the number of measurements $m$ is higher than the ambient dimension $n$ of the signal itself. For the cases where $m$ and $n$ are large, this requirement puts a heavy burden on computation and storage.

In contrast, our focus is on solving the the inverse problem (1) with very few number of samples, i.e., we are interested in the case $m\ll n$. While this makes the problem even more ill-posed, we show that such a barrier can be avoided if we assume that the underlying signal obeys a certain low-dimensional structure. In this paper, we focus on the sparsity assumption on the underlying signal, but our techniques could be extended to other signal structures. Further, for simplicity, we assume that our forward model is limited to only two modulo periods, as shown in the Fig. 2(a). Such a simplified variation of the modulo function already inherits much of the challenging aspects of the original recovery problem. Intuitively, this simplification requires that the value of dynamic range parameter $R$ should be large enough so that all the measurements $\langle\mathbf{a_{i}}\cdot\mathbf{x^{*}}\rangle$ can be covered within the domain of operation of the modulo function, i.e., $\langle\mathbf{a_{i}}\cdot\mathbf{x^{*}}\rangle\in[-R,R]~{}\forall i\in\{1,2,..,m\}$.

I-B Our contributions

In this paper, we propose a recovery algorithm for exact reconstruction of signals from modulo measurements of the form (1). We refer our algorithm as MoRAM, short for Modulo Recovery using Alternating Minimization. The key idea in our approach is to identify and draw parallels between modulo recovery and the problem of phase retrieval. Indeed, this connection enables us to bring in algorithmic ideas from classical phase retrieval, which also helps in our analysis.

Phase retrieval has its roots in several classical imaging problems, but has attracted renewed interest of late. There, we are given observations of the form:

[TABLE]

and are tasked with reconstructing $\mathbf{x^{*}}$. While these two different class of problems appear different at face value, the common theme is the need of undoing the effect of a piecewise linear transfer function applied to the observations. See Fig. 2 for a comparison. Both the functions are identical to the identity function in the positive half, but differ significantly in the negative half. Solving the phase retrieval problem is essentially equivalent to retrieving the phase ($\text{sign}\left(y_{i}\right)$) corresponding to each measurement $y_{i}$. However, the phase can take only two values: $1$ if $t\geq 0$, or $-1$ if $t<0$. Along the same lines, for modulo recovery case, the challenge is to identify the bin-index for each measurement. Estimating the bin-index correctly lets us “unravel” the modulo transfer function, thereby enabling signal recovery.

At the same time, several essential differences between the two problems restrict us from using phase retrieval algorithms as-is for the modulo reconstruction problem. The absolute value function can be represented as a multiplicative transfer function (with the multiplying factors being the signs of the linear measurements), while the modulo function adds a constant value ($R$) to negative inputs. Therefore, the estimation procedures propagate very differently in the two cases. In the case of phase retrieval, a wrongly estimated phase induces an error that increases linearly with the magnitude of each measurement. On the other hand, for modulo recovery problem, the error induced by an incorrect bin-index is $R$ (or larger), irrespective of the measurement. Therefore, existing algorithms for phase retrieval perform rather poorly for our problem (both in theory and practice).

We resolve this issue by making non-trivial modifications to existing phase retrieval algorithms that better exploit the structure of modulo reconstruction. We also provide analytical proofs for recovering the underlying signal using our algorithm, and show that such a recovery can be performed using an (essentially) optimal number of observations, provided certain standard assumptions are met. To the best of our knowledge we are the first to pursue this type of approach for modulo recovery problems with generic linear measurements, distinguishing us from previous work [ICCP15_Zhao, Bhandari].

I-C Techniques

The basic approach in our proposed (MoRAM) algorithm is similar to several recent non-convex phase retrieval approaches. We pursue two stages.

In the first stage, we identify a good initial estimated signal $\mathbf{x^{0}}$ that that lies (relatively) close to the true signal $\mathbf{x^{*}}$. A commonly used initialization technique for phase retrieval is spectral initialization as described in [netrapalli2013phase]. However, that does not seem to succeed in our case, due to markedly different behavior of the modulo transfer function. Instead, we introduce a novel approach of measurement correction by comparing our observed measurements with typical density plots of Gaussian observations. Given access to such corrected measurements, $\mathbf{x^{0}}$ can be calculated simply by using a first-order estimator. This method is intuitive, yet provides a provable guarantee for getting an initial vector that is close to the true signal.

In the second stage, we refine this coarse initial estimate to recover the true underlying signal. Again, we follow an alternating-minimization (AltMin) approach inspired from phase retrieval algorithms (such as [netrapalli2013phase]) that estimates the signal and the measurement bin-indices alternatively. However, as mentioned above, any estimation errors incurred in the first step induces fairly large additive errors (proportional to the dynamic range parameter $R$.) We resolve this issue by using a robust form of alternating-minimization (specifically, the Justice Pursuit algorithm [Laska2009]). We prove that AltMin, based on Justice Pursuit, succeeds provided the number of wrongly estimated bin-indices in the beginning is a small fraction of the total number of measurements. This gives us a natural radius for initialization, and also leads to provable sample-complexity upper bounds.

I-D Paper organization

The reminder of this paper is organized as follows. In Section II, we briefly discuss the prior work. Section III contains notation and mathematical model used for our analysis. In Section IV, we introduce the MoRAM algorithm and provide a theoretical analysis of its performance. We demonstrate the performance of our algorithm by providing series of numerical experiments in Section LABEL:sec:exp. Section LABEL:sec:disc provides concluding remarks.

II Prior work

We now provide a brief overview of related prior work. At a high level, our algorithmic development follows two (hitherto disconnected) streams of work in the signal processing literature.

Phase retrieval

As stated earlier, in this paper we borrow algorithmic ideas from previously proposed solutions for phase retrieval to solve the modulo recovery problem. Being a classical problem with a variety of applications, phase retrieval has been studied significantly in past few years. Approaches to solve this problem can be broadly classified into two categories: convex and non-convex.

Convex approaches usually consist of solving a constrained optimization problem after lifting the true signal $\mathbf{x^{*}}$ in higher dimensional space. The seminal PhaseLift formulation [candes2013phaselift] and its variations [gross2017improved], [candes2015phasediff] come under this category. Recently, [Bahmani2016PhaseRM, Goldstein2018PhaseMaxCP] proposed a convex optimization based algorithm for solving the phase retrieval problem that doesn’t lift the underlying signal and hence does not square the number of variables. The precise analysis of it has been derived in [Dhifallah2017PhaseRV, Salehi2018APA]. Typical non-convex approaches involve finding a good initialization, followed by iterative minimization of a loss function. Approaches based on Wirtinger Flow [candes2015phase, zhang2016reshaped, chen2015solving, cai2016optimal] and Amplitude flow [wang2016sparse, wang2016solving] come under this category.

In recent works, extending phase retrieval algorithms to situations where the underlying signal exhibits a sparse representation in some known basis has attracted interest. Convex approaches for sparse phase retrieval include [ohlsson2012cprl, li2013sparse, bahmani2015efficient, jaganathan2012recovery]. Similarly, non-convex approaches for sparse phase retrieval include [netrapalli2013phase, cai2016optimal, wang2016sparse]. Our approach in this paper towards solving the modulo recovery problem can be viewed as a complement to the non-convex sparse phase retrieval framework advocated in [Jagatap2017].

Modulo recovery

The modulo recovery problem is also known in the classical signal processing literature as phase unwrapping. The algorithm proposed in [bioucas2007phase] is specialized to images, and employs graph cuts for phase unwrapping from a single modulo measurement per pixel. However, the inherent assumption there is that the input image has very few sharp discontinuities, and this makes it unsuitable for practical situations with textured images. Our work is motivated by the recent work of [ICCP15_Zhao] on high dynamic range (HDR) imaging using a modulo camera sensor. For image reconstruction using multiple measurements, they propose the multi-shot UHDR recovery algorithm, with follow-ups developed further in [Lang2017]. However, the multi-shot approach depends on carefully designed camera exposures, while our approach succeeds for non-designed (generic) linear observations; moreover, they do not include sparsity in their model reconstructions. In our previous work [Shah], we proposed a different extension based on [ICCP15_Zhao, soltani2017stable] for signal recovery from quantized modulo measurements, which can also be adapted for sparse measurements, but there too the measurements need to be carefully designed.

In the literature, several authors have attempted to theoretically understand the modulo recovery problem. Given modulo-transformed time-domain samples of a band-limited function, [Bhandari] provides a stable algorithm for perfect recovery of the signal and also proves sufficiency conditions that guarantees the perfect recovery. [Cucuringu2017] formulates and solves an QCQP problem with non-convex constraints for denoising the modulo-1 samples of the unknown function along with providing a least-square based modulo recovery algorithm.. However, both these methods relay on the smoothness of the band-limited function as a prior structure on the signal, and as such it is unclear how to extend their use to more complex modeling priors (such as sparsity in a given basis).

In recent works, [Bhandari2018UnlimitedSO] proposed unlimited sampling algorithm for sparse signals. Similar to [Bhandari], it also exploits the bandlimitedness by considering the low-pass filtered version of the sparse signal, and thus differs from our random measurements setup. In [Musa2018GeneralizedAM], modulo recovery from Gaussian random measurements is considered, however, it assumes the true signal to be distributed as mixed Bernoulli-Gaussian distribution, which is impractical in real world imaging scenarios.

For a qualitative comparison of our MoRAM method with existing approaches, refer Table I. The table suggests that the previous approaches varied from the Nyquist-Shannon sampling setup only along the amplitude dimension, as they rely on band-limitedness of the signal and uniform sampling grid. We vary the sampling setup along both the amplitude and time dimensions by incorporating sparsity in our model, which enables us to work with non-uniform sampling grid (random measurements) and achieve a provable sub-Nyquist sample complexity.

III Preliminaries

III-A Notation

Let us introduce some notation. We denote matrices using bold capital-case letters ($\mathbf{A,B}$), column vectors using bold-small case letters ($\mathbf{x,y,z}$ etc.) and scalars using non-bold letters ($R,m$ etc.). We use letters $C$ and $c$ to represent constants that are large enough and small enough respectively. We use $\mathbf{x}^{\top},\mathbf{A}^{\top}$ to denote the transpose of the vector $\mathbf{x}$ and matrix $\mathbf{A}$ respectively. The cardinality of set $S$ is denoted by $\mathrm{card}(S)$. We define the signum function as $\mathrm{sgn}(x):=\frac{x}{|x|}$ for every $x\in{\mathbb{R}},x\neq 0$, with the convention that $\mathrm{sgn}(0)=1$. The $i^{th}$ element of the vector $\mathbf{x}\in{\mathbb{R}}^{n}$ is denoted by $x_{i}$. Similarly, $i^{th}$ row of the matrix $\mathbf{A}\in{\mathbb{R}}^{m\times n}$ is denoted by $\mathbf{a_{i}}$, while the element of $\mathbf{A}$ in the $i^{th}$ row and $j^{th}$ column is denoted as $a_{ij}$. The projection of $\mathbf{x}\in{\mathbb{R}}^{n}$ onto a set of coordinates $S$ is represented as $\mathbf{x}_{S}\in{\mathbb{R}}^{n}$, i.e., $\mathbf{x}_{S_{j}}=\mathbf{x}_{j}$ for $j\in S$, and [math] elsewhere.

III-B Mathematical model

As depicted in Fig. 2(a), we consider the modulo operation within 2 periods (one in the positive half and one in the negative half). We assume that the value of dynamic range parameter $R$ is large enough so that all the measurements $\langle\mathbf{a_{i}}\cdot\mathbf{x^{*}}\rangle$ are covered within the domain of operation of modulo function. Rewriting in terms of the signum function, the (variation of) modulo function under consideration can be defined as:

[TABLE]

One can easily notice that the modulo operation in this case is nothing but an addition of scalar $R$ if the input is negative, while the non-negative inputs remain unaffected by it. If we divide the number line in these two bins, then the coefficient of $R$ in above equation can be seen as a bin-index, a binary variable which takes value [math] when $\mathrm{sgn}(t)=1$, or $1$ when $\mathrm{sgn}(t)=-1$. Inserting the definition of $f$ in the measurement model of Eq. 1 gives,

[TABLE]

We can rewrite Eq. 2 using a bin-index vector $\mathbf{p}\in\{0,1\}^{m}$. Each element of the true bin-index vector $\mathbf{p}^{*}$ is given as,

[TABLE]

If we ignore the presence of modulo operation in above formulation, then it reduces to a standard compressive sensing reconstruction problem. In that case, the compressed measurements $y_{c_{i}}$ would just be equal to $\langle\mathbf{a_{i}}\cdot\mathbf{x^{*}}\rangle$.

While we have access only to the compressed modulo measurements $\mathbf{y}$, it is useful to write $\mathbf{y}$ in terms of true compressed measurements $\mathbf{y}_{c}$. Thus,

[TABLE]

It is evident that if we can recover $\mathbf{p^{*}}$ successfully, we can calculate the true compressed measurements $\langle\mathbf{a_{i}}\cdot\mathbf{x^{*}}\rangle$ and use them to reconstruct $\mathbf{x^{*}}$ with any sparse recovery algorithm such as CoSaMP [needell2010cosamp] or basis-pursuit [chen2001atomic, spgl1:2007, BergFriedlander:2008].

IV Sparse signal recovery

Of course, the major challenge is that we do not know the bin-index vector. In this section, we describe our algorithm to recover both $\mathbf{x^{*}}$ and $\mathbf{p^{*}}$, given $\mathbf{y,A,s,R}$. Our algorithm MoRAM (Modulo Reconstruction with Alternating Minimization) comprises of two stages: (i) an initialization stage, and (ii) descent stage via alternating minimization.

IV-A Initialization by re-calculating the measurements

Similar to other non-convex approaches, MoRAM also requires an initial estimate $\mathbf{{x}^{0}}$ that is close to the true signal $\mathbf{{x}^{*}}$. We have several initialization techniques available; in phase retrieval, techniques such as spectral initialization are often used. However, the nature of the problem in our case is fundamentally different due to the non-linear additive behavior of the modulo transfer function. To overcome this issue, we propose a method to re-calculate the true Gaussian measurements ($\mathbf{y_{c}}=\mathbf{Ax^{*}}$) from the available modulo measurements.

The high level idea is to undo the nonlinear effect of modulo operation in a significant fraction of the total available measurements. To understand the method for such re-calculation, we will first try to understand the effect of modulo operation on the linear measurements.

IV-A1 Effect of the modulo transfer function

To provide some intuition, let us first examine the distribution of the $\mathbf{Ax^{*}}$(Fig. 3) and $\mathbf{\mod(\mathbf{Ax^{*}})}$ (shown in Fig. 4) to understand what information can be obtained from the modulo measurements. We are particularly interested in the case where the elements of $\mathbf{Ax^{*}}$ are small compared to the modulo range parameter $R$.

Note that the compressed measurements $\mathbf{y_{c}}$ follow the standard normal distribution, as $\mathbf{A}$ is Gaussian random matrix. These plots essentially depict the distribution of our observations before and after the modulo operation.

With reference to Fig. 3, we divide the compressed observations $\mathbf{y_{c}}$ in two sets: $\mathbf{y_{c,+}}$ contains all the non-negative observations (orange) with bin-index$=0$, while $\mathbf{y_{c,-}}$ contains all the negative ones (green) with bin-index$=1$.

As shown in Fig. 4, after modulo operation, the set $\mathbf{y}_{c,-}$ (green) shifts to the right by $R$ and gets concentrated in the right half ($[R/2,R]$); while the set $\mathbf{y}_{c,+}$ (orange) remains unaffected and concentrated in the left half ($[0,R/2]$). Thus, for some of the modulo measurements, their correct bin-index can be identified just by observing their magnitudes relative to the midpoint $R/2$. This leads us to obtain following maximum likelihood estimator for bin-indices ($\mathbf{p}$):

[TABLE]

The $\mathbf{p^{init}}$ obtained with above method contains the correct values of bin-indices for many of the measurements, except for the ones concentrated within the ambiguous region in the center.

Once we identify the initial values of bin-index for the modulo measurements, we can calculate corrected measurements as,

[TABLE]

We use these corrected measurements $\mathbf{y^{init}_{c}}$ to calculate the initial estimate $\mathbf{{x}^{0}}$ with first order unbiased estimator.

[TABLE]

where $H_{s}$ denotes the hard thresholding operator that keeps the $s$ largest absolute entries of a vector and sets the other entries to zero.

IV-B Alternating Minimization

Using Eq. 5, we calculate the initial estimate of the signal $\mathbf{{x}^{0}}$ which is relatively close to the true vector $\mathbf{x^{*}}$. Starting with $\mathbf{{x}^{0}}$, we calculate the estimates of $\mathbf{p}$ and $\mathbf{x}$ in an alternating fashion to converge to the original signal $\mathbf{x^{*}}$. At each iteration of alternating-minimization, we use the current estimate of the signal ${\mathbf{x^{t}}}$ to get the value of the bin-index vector $\mathbf{{p}^{t}}$ as following:

[TABLE]

Given that $\mathbf{x^{0}}$ is close to $\mathbf{x^{*}}$, we expect that $\mathbf{p^{0}}$ would also be close to $\mathbf{p^{*}}$. Ideally, we would calculate the correct compressed measurements $\mathbf{y^{t}_{c}}$ using $\mathbf{p^{t}}$, and use $\mathbf{y^{t}_{c}}$ with any popular compressive recovery algorithms such as CoSaMP or basis pursuit to calculate the next estimate $\mathbf{{x}^{t+1}}$. Thus,

[TABLE]

[TABLE]

where $\mathcal{M}_{s}$ denotes the set of $s$-sparse vectors in $\mathbb{R}^{n}$. Note that sparsity is only one of several signal models that can be used here, and in principle a rather similar formulation would extend to cases where $\mathcal{M}$ denotes any other structured sparsity model [modelcs].

However, it should be noted that the “bin” error $\mathbf{d^{t}}=\mathbf{p^{t}-p^{*}}$, even if small, would significantly impact the correction step that constructs $\mathbf{y^{t}_{c}}$, as each incorrect bin-index would add a noise of the magnitude $R$ in $\mathbf{y^{t}_{c}}$. Our experiments suggest that the typical sparse recovery algorithms are not robust enough to cope up with such large errors in $\mathbf{y^{t}_{c}}$. To tackle this issue, we employ an outlier-robust sparse recovery method [Laska2009]. We consider the fact that the nature of the error $\mathbf{d^{t}}$ is sparse with sparsity $s_{dt}=\|\mathbf{d^{t}}\|_{0}$; and each erroneous element of $\mathbf{p}$ adds a noise of the magnitude $R$ in $\mathbf{y^{t}_{c}}$.

Rewriting in terms of Justice Pursuit, the recovery problem now becomes problem becomes,

[TABLE]

However, the sparsity of $\mathbf{d^{t}}$ is unknown, suggesting that greedy sparse recovery methods cannot be directly used without an additional hyper-parameter. Thus, we employ basis pursuit [candes2006compressive] which does not rely on sparsity. The robust formulation of basis pursuit is referred as Justice Pursuit (JP) [Laska2009], specified in Eq. 7.

[TABLE]

Proceeding this way, we repeat the steps of bin-index calculation (as in Eq. 6) and sparse recovery (Eq. 7) altenatingly for $\mathrm{T}$ iterations. Our algorithm is able to achieve convergence to the true underlying signal, as supported by the results in the experiments section.

V Mathematical Analysis

Before presenting experimental validation of our proposed MoRAM algorithm, we now perform a theoretical analysis of the descent stage of our algorithm. We assume the availability of an initial estimate $\mathbf{x^{0}}$ that is close to $\mathbf{x^{*}}$, i.e. $\|\mathbf{x^{0}}-\mathbf{x^{*}}\|_{2}\leq\delta\|\mathbf{x^{*}}\|_{2}$. In our case, our initialization step (in Alg. 2) provide such $\mathbf{x^{0}}$.

We perform alternating-minimization as described in 2, starting with $\mathbf{x^{0}}$ calculated using Alg. 1. For simplicity, we limit our analysis of the convergence to only one AltMin iteration. In fact, according to our theoretical analysis, if initialized closely enough, one iteration of AltMin suffices for exact signal recovery with sufficiently many measurements. However, in practice we have observed that our algorithm requires more than one AltMin iterations.

The first step is to obtain the initial guess of the bin-index vector (say $\mathbf{p^{0}}$) using $\mathbf{x^{0}}$.

[TABLE]

If we try to undo the effect of modulo operation by adding back $R$ for the affected measurements based on the bin-index vector $\mathbf{p^{0}}$, it would introduce an additive error equal to $R$ corresponding to each of the incorrect bin-indices in $\mathbf{p^{0}}$.

[TABLE]

We show the guaranteed recovery of the true signal as the corruption in the first set of corrected measurements $\mathbf{y_{c}^{0}}$ can be modeled as sparse vector with sparsity less than or equal to $\lambda m$, with $c$ being a fraction that can be explicitly bounded.

To prove this, we first introduce the concept of binary $\epsilon$-stable embedding as proposed by [Jacques2013]. Let $\mathcal{B}^{m}$ be a Boolean cube defined as $\mathcal{B}^{m}:=\{-1,1\}^{m}$ and let $S^{n-1}:=\{\mathbf{x}\in{\mathbb{R}}^{n}:\|\mathbf{x}\|_{2}=1\}$ be the unit hyper-sphere of dimension $n$.

Definition 1** (Binary $\epsilon$-Stable Embedding).**

A mapping $F:{\mathbb{R}}^{n}\rightarrow\mathcal{B}^{m}$ is a binary $\epsilon$-stable embedding (B$\epsilon$SE) of order $s$ for sparse vectors if:

[TABLE]

for all $\mathbf{x,y}\in S^{n-1}$ with $|supp(x)\cup supp(y)|\leq s$.

In our case, let us define the mapping $F:{\mathbb{R}}^{n}\rightarrow\mathcal{B}^{m}$ as:

[TABLE]

with $\mathbf{A}\sim\mathcal{N}^{m\times n}(0,1)$. We obtain: