Sparse signal recovery from modulo observations

We consider the problem of reconstructing a signal from under-determined modulo observations (or measurements). This observation model is inspired by a relatively new 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 signal recovery problem under sparsity constraints for the special case to modulo folding limited to two periods. We show that given a sufficient number of measurements, our algorithm perfectly recovers the underlying signal. We also provide experiments validating our approach on toy signal and image data and demonstrate its promising performance.

a The modulo transfer function. b The modulo function limited to 2 modulo periods-a specific case that we consider in this paper range; however, real-world signals can contain a large range of intensity levels, and if tuned incorrectly, signal measurements 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 the 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 the sensor dynamic range via enhanced hardware; this, of course, can be untenably expensive. An intriguing alternative is to deploy special digital modulo sensors [2][3][4][5]. 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 [6,7] 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 inverse problem (1) with very few samples, i.e., we are interested in the case m 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 special case when the underlying signal is sparse, but our techniques could be extended to other signal structures. Further, for simplicity, we assume that our observation model is limited to only two modulo fold periods, one for positive-and one for negative-valued coefficients. This does not reflect practice, but we will see that such a 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 finite, but large enough that most of the measurements fall within the interval [ −R, R]. We emphasize that such requirement does not put a hard constraint in our algorithm, and successful recovery is possible even if some measurements lie outside the interval covered by two modulo periods.

Our contributions
In this paper, we propose a recovery algorithm for exact reconstruction of sparse signals from modulo measurements of the form (1) with the modulo fold operation being limited to two periods. We refer to our algorithm as MoRAM, short for Modulo Recovery using (2021)  Alternating Minimization. We observe that the modulo operation with respect to parameter R can be seen as the addition or subtraction of an integer multiple of R from the input. We refer to this integer multiple as the bin-index p. It is not hard to observe that a successful recovery of bin-index can lead to successful recovery of the input using standard sparse recovery methods. Concretely, the forward model in (1) can be written as follows: As we restrict our modulo operation to only two periods, the bin-index can only take two values: 0 for non-negative inputs, and 1 for negative inputs. We discuss extending this to multiple modulo fold periods in Section 4.
As mentioned above, 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. However, due to the highly nonlinear nature of the forward model, each incorrect bin-index introduces a gross additive error of magnitude R irrespective of the magnitude of the measurement. Commonly used signal recovery algorithms fail to cope up with such gross corruptions of large magnitudes, and perform poorly on our problem (both in theory and practice). Thus, it becomes even more important to be able to find a good estimate of the bin-index and to employ a robust signal recovery algorithm that can recover the signal even in the presence of grossly corrupted measurements.
To this end, our proposed algorithm follows two steps. We first leverage the structure of the modulo folding operator to obtain a promising initial estimate of the bin-indices. We then plug this estimate into a robust sparse recovery formulation that gives us the final signal estimate. We provide analytical correctness proofs for both states and show that signal recovery can be performed using an (essentially) optimal number of observations provided certain assumptions are met. To the best of our knowledge, we are the first to pursue this type of approach for modulo recovery problems with Gaussian linear measurements, distinguishing us from previous work [6,7].

Techniques
Our proposed algorithm (MoRAM) consists of two stages.
In the first stage, we identify a good initial estimated bin-index p 0 that is relatively close to the true bin-index p * . To obtain this estimate, we employ a simple bin-index estimation step by comparing our observed measurements with typical density plots of Gaussian observations. This method is able to recover a large number of bin-indices correctly, and also provides a provable upper bound on the number of erroneous bin-indices.
In the second stage, we use the initial estimate of the bin-indices of the measurements to recover the true underlying signal. We follow an alternating minimization (AltMin) approach inspired from phase retrieval algorithms (such as [9]) that estimates the signal and the measurement of 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, using the Justice Pursuit algorithm [10]). We prove that given enough number of measurements, our Justice Pursuit-based AltMin approach succeeds provided the number of wrongly estimated bin-indices in the beginning is a sufficiently small fraction of the total number of measurements. This gives us a natural radius of initialization for the initial bin-index estimate and also leads to provable sample-complexity upper bounds.

Prior work
Since signal recovery from nonlinear measurements is a very large and classical area of study, our review of prior work will unfortunately be incomplete.

Modulo recovery
The modulo recovery problem shares some similarities with the problem of phase unwrapping from the classical signal processing literature that would allow one to use apply phase unwarapping methods to modulo recovery; however, these methods do not provide proven guarantees. For example, the algorithm proposed in [11] 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 [7] 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 [12]. However, the multi-shot approach depends on carefully designed camera exposures, while our approach succeeds for non-designed random Gaussian linear observations; moreover, they do not include sparsity in their model reconstructions. In our previous work [8], we proposed a different extension based on [7,13] 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, [6,14,15] provide a stable algorithm for signal recovery and also proves sufficiency conditions that guarantees the recovery. Ordentlich et al. [14] use non-modulo data for the initialization to exploit the statistical structure in order to undo the effects of the modulo operation. Cucuringu and Tyagi [16] formulate and solve a 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 rely 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). On the contrast, our approach does not rely on such smoothness assumptions.
In recent works, [17] proposed unlimited sampling algorithm for sparse signals and images. Similar to [6], it also exploits the band-limitedness by considering the low-pass filtered version of the sparse signal and thus differs from our random measurements setup. In [18], modulo recovery from Gaussian random measurements is considered. However, it assumes the true signal to be distributed as a mixed Bernoulli-Gaussian distribution which is not a standard assumption.
For a qualitative comparison of our MoRAM method with existing approaches, we refer the reader to 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 a non-uniform sampling grid (random measurements) and achieve a provable sub-Nyquist sample complexity.

Preliminaries
Let us introduce some notation. We denote matrices using bold capital-case letters (A, B), column vectors using bold-small case letters (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 x and A to denote the transpose of the vector x and matrix A respectively. The cardinality of set S is denoted by card(S). We define the signum function as sgn(x) := x |x| for every x ∈ R, x = 0, with the convention that sgn(0) = 1. The ith element of the vector x ∈ R n is denoted by x i . Similarly, the ith row of the matrix A ∈ R m×n is denoted by a i , while the element of A in the ith row and jth column is denoted as a ij .

Mathematical model
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 most of the measurements a i , x * 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: 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 0 when sgn(t) = 1, or 1 when sgn(t) = −1.
Inserting the definition of f in the measurement model of Eq. 1 gives, We can rewrite Eq. 2 using a bin-index vector p ∈ {0, 1} m . Each element of the true bin-index vector p * is given as, If we ignore the presence of the modulo operation in the 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 a i , x * . While we have access only to the compressed modulo measurements y, it is useful to write y in terms of true compressed measurements y c . Thus, It is evident that if we can recover p * successfully, we can calculate the true compressed measurements a i , x * and use them to reconstruct x * with any sparse recovery algorithm such as CoSaMP [19] or basis pursuit [20][21][22].

Signal recovery
The major barrier to signal recovery is that the bin-index vector is unknown. In this section, we describe our algorithm to recover both x * and p * , given the modulo measurements y, measurement matrix A, sparsity of underlying signal s, and the modulo parameter R. In this work, we rely on the assumption that our signal is sparse in a known domain with sparsity being s. Our algorithm MoRAM (Modulo Reconstruction with Alternating Minimization) comprises two stages: (i) an bin-index initialization stage and (ii) a descent stage via alternating minimization.

Bin-index initialization
As stated earlier, if we recover true bin-index p * successfully, x * can be recovered easily using any sparse recovery algorithm as we can obtain the true compressed measurements a i , x * from p * . Thus, in the absence of p * , we propose to estimate a fraction of the values from the p * correctly. To understand the rationale for such a procedure, we will first try to understand the effect of the modulo operation on the linear measurements.

Effect of the modulo transfer function
To provide some intuition, let us first examine the relation between the distributions of Ax * and mod (Ax * ). It is easy to see that the compressed measurements y c follow a normal distribution.
We can now divide the compressed observations y c into two sets: y c,+ , which contains all the non-negative observations with bin-index= 0, and y c,− , which contains all the negative observations with bin-index= 1. As shown in Fig. 2, after the modulo operation, the set y c,− (green) shifts to the right by R and gets concentrated in the right half ([ R/2, R]), while the set 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 by observing their magnitudes relative to the midpoint R/2. This leads us to the following estimator for bin-indices (p): The vector p 0 obtained with the above method contains the correct values of binindices for many of the measurements, except for the ones concentrated within the ambiguous region in the center. We should highlight that the procedure in Eq. 3 will succeed only for the specific case of modulo fold operations limited to two periods, one for the positive and one for the negative cycle.
Once we identify the initial values of bin-index for the modulo measurements, we can calculate corrected measurements as:

Alternating minimization
Using Eq. 3, we calculate the initial estimate of the bin-index p 0 in which significant fraction of the total values are estimated correctly. Starting with p 0 , we calculate the estimates of x and p in an alternating fashion to converge to the original signal x * and true bin-index p * . With p t being close to p * , we would calculate the correct compressed measurements y t c using p t and use y t c with any popular compressive recovery algorithms (such as CoSaMP, or basis pursuit) to calculate the signal estimate x t . Therefore: where M s denotes the set of s-sparse vectors in R n . Note that sparsity is only one of several signal models that can be used here, and a rather similar formulation would extend to cases where M denotes any other structured sparsity model [23,24]. However, the bin-index estimation error, d t = p t − p * , even if small, would significantly impact the correction step that constructs y t c since each incorrect bin-index would add a noise of the magnitude R in y t c . Our experiments suggest that the typical sparse recovery algorithms are not robust enough to cope up with such large errors in y t c . To tackle this issue, we employ an outlier-robust sparse recovery method known as Justice Pursuit [10].
At a high level, Justice Pursuit tackles the problem of sparse signal recovery from the measurements that are corrupted by a sparse but large (unbounded) corruptions. Justice Pursuit leverages the fact that the corruptions are also sparse, and reformulates the problem to recover both the sparse signal and sparse corruptions together in the form of a concatenated sparse vector. In our case, the error d t is sparse with sparsity s dt = d t 0 , and each erroneous element of p adds a corruption of magnitude R in y t c . Following [10], we augment the measurement matrix A with an identity matrix I m×m and introduce an intermediate vector u ∈ R n+m to represent our measurements at iteration t as: and solve for the (s + s dt )−sparse estimate u: x t d t = u = arg min Here, the signal estimate x t is obtained by selecting the first n elements of u, while an estimate of the corruptions can be obtained by selecting the last m elements of u. The problem in Eq. 8 can be solved by any stable sparse recovery algorithm such as CoSaMP or IHT; however, note that the sparsity of d t is unknown, suggesting that greedy sparse recovery methods cannot be directly used without an additional hyper-parameter. Therefore, we employ basis pursuit [25] which does not heavily depend on a priori knowledge of the sparsity level. We refer to the routine that solves the program in Eq. 8 using basis pursuit as JP. Given A, y t c , JP returns x t . Thus, Once the signal estimate x t is obtained at each iteration of alternating minimization, we use it to calculate the value of the bin-index vector p t+1 as follows: Proceeding this way, we repeat the steps of sparse recovery (Eq. 8) and bin-index calculation (Eq. 10) in alternating fashion for T iterations. Under certain conditions (described in Section 2.4 below), our algorithm is able to achieve convergence to the true underlying signal.

Mathematical analysis
In this section, we provide correctness proofs for both steps of Algorithm 1. For the first stage, we derive an upper bound on the number of incorrect estimations in p 0 obtained in the bin-index initialization step. This upper bound essentially provides an upper bound on the permissible sparsity of d 0 . For the second stage, we calculate a sufficient number of measurements required such that the augmented matrix used in the Justice Pursuit formulation in (8) satisfies the Restricted Isometry Property (RIP), which would in turn enable a recovery guarantee.

Bin-index initialization
In this step, we initialize the bin-index vector p 0 according to Eq. 3. We can also quantify the number of correctly estimated bin-indices by calculating the area under the curve of the density plots of the measurements before and after the modulo operation. An illustration is provided in Fig. 3.
In this analysis, our goal is to characterize the distribution of total number of measurements for which we can estimate the correct bin-index through Eq. 3. Such a random variable is denoted by M c . From M c , we can calculate the sparsity of d 0 as d 0 0 = m−M c . The following lemma presents a bound on the sparsity of d 0 .

Lemma 1 Let the entries of the measurement matrix be generated as
Here, φ(·) is a Gaussian density with mean μ = 0 and variance σ 2 = x * 2 2 . Proof Observe that each element of A is i.i.d. standard normal, i.e., μ A ij = 0 and σ 2 A ij = 1. Recall that Therefore, we have Thus, each element of y c follows a zero-mean Gaussian distribution with variance σ 2 . Let E i be the event that the random variable y c,i lies in the interval [ −R/2, R/2]; this event indicates that the corresponding measurement is appropriately corrected using Eq. 4 Clearly, E i is a Bernoulli random variable with probability q Elementary probability calculations give us: where Q 0,σ 2 (·) is the usual Q-function. This is not calculable in closed form; however, it can be lower bounded using the following identity (where φ 0,σ 2 (·) is a Gaussian density function with mean zero and variance σ 2 : The random variable M c = m i=1 E i denotes the number of corrected measurements. By an application of the Chernoff bound, for any δ ∈ (0, 1),where μ is the mean of M c . Plugging in μ = mq gives the desired result.
We now perform a theoretical analysis of the descent stage of our algorithm. We assume the availability of an initial estimate of bin-index vector p 0 that is close to p * . In our case, our initialization step (in Alg. 1) provides such p 0 .
We perform alternating minimization (AltMin) as described in 1, starting with p 0 calculated using Eq. 3. For simplicity, we limit our analysis of the convergence to only one AltMin iteration. In fact, according to our theoretical analysis, if initialized well enough, one iteration of AltMin suffices for exact signal recovery with sufficiently many measurements; however, in practice, we have observed that our algorithm performs better with multiple AltMin iterations. Proof In the estimation step, Algorithm 1 recasts the problem of recovering the true signal x * as a special case of sparse signal recovery from sparsely corrupted compressive measurements. The presence of modulo operation modifies the compressive measurements by adding a constant noise of the value R in fraction of total measurements. However, once we identify correct bin-index for some of the measurements using Eq. 3, the remaining noise can be modeled as sparse corruptions d according to the formulation: Here, the 0 -norm of d 0 gives us the number of noisy measurements in y 0 c . If the initial bin-index vector p 0 is close to the true bin-index vector p * , then d 0 0 is small enough with respect to total number of measurements m; thus, d 0 can be treated as sparse corruption. If we model this corruption as a sparse noise, then we can employ JP for a guaranteed recovery of the true signal given sufficiently large number of measurements are available. Denote d 0 0 = m − M c as number of measurements for which the bin-index estimates were incorrect. Then, using Lemma 1, with probability at least 1 − e −O(mδ 2 ) : Algorithm 1 is essentially the Justice Pursuit (JP) formulation as described in [10]. Exact signal recovery from sparsely corrupted measurements is a well-studied problem with uniform recovery guarantees available in the existing literature. We use the guarantee proved in [10] for Gaussian observations, which states that provided enough measurements, the augmented matrix [ A I] satisfies the Restricted Isometry Property. As stated in [26], one can recover a sparse signal exactly by tractable 1minimization if the measurement matrix is known to satisfy the RIP. Thus, provided m ≥ C x * 0 + d 0 0 log (n + m)/ x * 0 + d 0 0 , we invoke Theorem 1.1 from [10] and replace d 0 0 with m (1 − U + δU) as stated above to complete the proof.
From the theorem, we see that the number of measurements required for guaranteed recovery depends on the ratio of σ (standard deviation of the measurements) and R. In practical applications, choosing a sufficiently large R such that the interval [ −R, R] covers multiple standard deviations on both sides of origin enables successful recovery.

Experiments
In this section, we present the results of simulations of signal reconstruction using our algorithm. All numerical experiments were conducted using MATLAB R2020b on a Windows system with an Intel CPU and 16GB RAM. Our experiments explores the performance of the MoRAM algorithm on both synthetic data and real images.
We perform experiments on a synthetic sparse signal x * ∈ R n with n = 1000. The sparsity level of the signal is chosen in steps of 6 starting from 6 with a maximum value of 24. The non-zero elements of the test signal x * are generated using a zero-mean Gaussian distribution N (0, 1) and normalized such that x * = 1. The elements of the Gaussian measurement matrix A ∈ R m×n , a ij are also generated using the standard normal distribution N (0, 1/m). The number of measurements m is varied from m = 100 to m = 1000 in steps of 100. Using A, x * , and R, we first obtain the compressed modulo measurements y by passing the signal through forward model described by Eq. 2. For reconstruction, algorithm 1 is employed. We plot the variation of the relative reconstruction error x * −x T x * with the number of measurements m for our AltMin-based sparse recovery algorithm MoRAM.
For each combination of R, m, and s, we run 10 independent Monte Carlo trials and calculate mean of the relative reconstruction error over these trials. Figure 4a and b illustrate the performance of our algorithm for two values of R respectively. Additionally, in Fig. 5, we also show reconstruction performance results for higher sparsity values s = 20, 30, 40, 50 as they vary with the number of measurements. It is evident that for each combination of R and s, our algorithm recovers the true signal (with zero relative error) provided enough measurements. In all such cases, the minimum number of measurements required for exact recovery is well below the ambient dimension (n) of the underlying signal. From the analysis provided in the previous section, it is clear that the success of our algorithm depends on the ratio of the standard deviation σ of the Gaussian measurements and modulo period R. When the matrix A is correctly normalized, σ directly depends on the Euclidean norm of the signal, which in turn depends on the sparsity. Thus, as discussed earlier, if R is chosen to be large enough, success of our algorithm is guaranteed.
To put this into perspective, we provide additional results where we fix the number of measurements m and sparsity s, and vary the modulo parameter R. These values suggest that even when the R becomes small and the ratio of σ and R increases, our algorithm is able to recover the underlying signal perfectly as shown in Fig. 6. As R increases, recovery becomes easier to achieve for our algorithm. It also shows that the performance of our algorithm decays gracefully when varying R.
In practical scenarios, we may encounter cases where a few measurements exceed the modulo period R. If we define γ as max (| A, x * |), then the value of γ would vary across each Monte Carlo run of our experiment (due to the randomness in the measurement matrix A), and we may encounter cases with γ > R as well. To analyze such cases, we provide an additional experiment where we fix the number of measurements m to 400, sparsity s to 12, and modulo parameter R to 3.2 and run our recovery algorithm multiple times with different realizations of A. We note down the values of γ during these experiments and note summary statistics as mean γ = 3.45, variance γ = 0.2065, max γ = 4.41, min γ = 2.82 across 50 random trials. In all these 50 cases, our algorithm recovered the true signal perfectly, i.e., resulting in zero relative error in each case. As the γ varies widely and also takes values higher than R, perfect recovery achieved in all cases shows that our algorithm is robust to measurements with magnitudes exceeding R, as far as such measurements are low in number (as guaranteed by Gaussian tail bounds on A, x * ). In this sense, our algorithm follows a "graceful decay" with respect to the assumption that most of the measurements are contained in the interval [ −R, R] in this sense.
We also evaluated the performance of our algorithm on a real image. We obtain sparse representation of the real image by transforming the original image in the wavelet basis (db1). The image used in our experiment is a 128 × 128 (n = 16384) image (Fig. 7a), We reconstruct the image with MoRAM using m = 4000 and m = 6000 compressed modulo measurements, for 3 different values of R, 4, 4.25, and 4.5. As expected, the reconstruction performance increases with increasing value of R. As shown in Fig. 7 (bottom), for m = 6000, the algorithm produces near-perfect recovery for all 3 values of R with high PSNR. Here, let us note that the blocky artifacts appearing in the recovered image are actually due to sparsification of the original image using the Haar wavelet transform. Since we use s = 800 to obtain the ground truth image, the ground truth itself contains some compression artifacts as depicted in Fig. 7 (a, bottom). The PSNR values are calculated with respect to the sparse ground truth image and not with respect to the original image as our algorithm aims to recover only the former. We expect the effect of compression artifacts to decrease with a better choice of sparsity basis for the underlying image.

Conclusions
In this paper, we presented a novel algorithmic approach for sparse signal recovery from compressed modulo measurements, inspired by techniques from phase retrieval. We also support our proposed algorithm via mathematical analysis and several experimental results. Our work points the way to a few directions for further research. While in this paper we considered only two modulo periods, extending the proposed approach for more periods (up to a theoretically infinite number) is a significant and interesting research direction. Also, instead of relying on a sparsity prior for compressed recovery, employing richer priors such as GANs [27][28][29] is an additional direction. Moreover, our analysis is limited to the case of Gaussian measurements schemes, which may or may not be physically realizable. Extending our results to more practical measurement schemes such as Fourier-based sampling or ptychography [30] can be an interesting problem for future study.