Adaptive recovery of dictionary-sparse signals using binary measurements

One-bit compressive sensing (CS) is an advanced version of sparse recovery in which the sparse signal of interest can be recovered from extremely quantized measurements. Namely, only the sign of each measure is available to us. The ground-truth signal is not sparse in many applications yet can be represented in a redundant dictionary. A strong line of research has addressed conventional CS in this signal model, including its extension to one-bit measurements. However, one-bit CS suffers from an extremely large number of required measurements to achieve a predefined reconstruction error level. A common alternative to resolve this issue is to exploit adaptive schemes. We utilize an adaptive sampling strategy to recover dictionary-sparse signals from binary measurements in this work. A multi-dimensional threshold is proposed for this task to incorporate the previous signal estimates into the current sampling procedure. This strategy substantially reduces the required number of measurements for exact recovery. We show that the proposed algorithm considerably outperforms the state-of-the-art approaches through rigorous and numerical analysis.

counts the number of nonzero elements. 1It has been shown that O(s log( N s )) measure- ments are sufficient to guarantee exact recovery of the signal, by solving the convex program: with high probability (see [1,2]).
Practical limitations force us to quantize the measurements in (1) as y = Q(A x) where Q : R m → A m is a nonlinear operator that maps the measurements into a finite symbol alphabet A .It is an interesting question to ask: What is the result of extreme quantiza- tion?[3] addressed this question: Signal reconstruction is still feasible using only one-bit quantized measurements.In one-bit compressed sensing, samples are taken as the sign of a linear transform of the signal y = sign(Ax) .This sampling scheme discards mag- nitude information.Therefore, we can only recover the direction of the signal.Fortunately, we can keep the amplitude information by using nonzero threshold.Thus, the new sampling scheme is y = [sign](Ax − τ ) where τ is the threshold vector.In our work, each element of τ is generated via τ i ∼ N (0, 1) .While a great part of CS literature dis- cusses sparse signals, most natural signals are dictionary-sparse, i.e., sparse in a transform domain.For instance, sinusoidal signals and natural image are sparse in Fourier and wavelet domains, respectively [4][5][6][7].This means that the signal of interest f ∈ R n can be expressed as f = Dx where D ∈ R n×N is a redundant dictionary with the con- strain DD H = I 2 and x ∈ R N is a sparse vector.With this assumption, y = Af = ADx gives the measurement vector.A common approach for recovering such signals is to use the optimization problem which is called ℓ 1 analysis problem [5,6].
In this work, we investigate a more practical situation where the signal of interest f is effective s-analysis-sparse which means that f satisfies �D H f � 1 ≤ √ s�D H f � 2 .In fact, perfect dictionary sparsity is rarely satisfied in practice, since real-world signals of interest are only compressible in a domain.Our approach is adaptive which means that we incorporate previous signal estimates into the current sampling procedure.More explicitly, we solve the optimization problem where ϕ ∈ R m is a vector of thresholds chosen adaptively based on previous estimations via Algorithm 1.We propose a strategy to find a best effective s-analysis-sparse approximation to a signal in R n . (2) D) H stands for the Hermitian of the matrix D 1 The used notation is introduced in Sect.1.2.

Contributions
In this section, we state our novelties compared to the previous works.To highlight the contributions, we list them as below.
1 Proposing a novel algorithm for dictionary-sparse signals: We introduce an adaptive thresholding algorithm for reconstructing dictionary-sparse signals in case of binary measurements.The proposed algorithm provides accurate signal estimation even in case of redundant and coherent dictionaries.The required number of one-bit measurements considerably outperforms the non-adaptive approach used in [8]. 2 Exponential decay of reconstruction error: The reconstruction error of our algorithm reduces exponentially when the number of adaptive stages (i.e., T) increases.
To be more precise, we obtain a near-optimal relation between the reconstruction error and the required number of adaptive batches.Written in mathematical form, if one takes the output of our reconstruction algorithm by f T , then we show that , where f is the ground-truth signal and T is the number of stages in our adaptive algorithm (see Theorem 1 for more details).3 High-dimensional threshold selection: We propose an adaptive high-dimensional threshold to extract the most information from each sample, which substantially improves performance and reduces the reconstruction error (see Algorithm 1 for more explanations).

Prior works and key differences
In this section, we review prior works about applying quantized measurements to CS framework [3,[8][9][10][11][12][13]. In what follows, we explain some of them.The authors of [3] propose a heuristic algorithm to reconstruct the ground-truth sparse signal from extreme quantized measurements, i.e., one-bit measurements.In [9], it has been shown that conventional CS algorithms also work well when the measurements are quantized.In [10], an algorithm with simple implementation is proposed.This algorithm has less error in Hamming distance than the existing ones.Investigated from a geometric view, the authors of [11] exploit functional analysis tools to provide an almost optimal solution to the problem of one-bit CS.They show that the number of required one-bit measurements is O(s log 2 ( n s )).The work of [14] presents two algorithms for full (i.e., direction and norm) reconstruction with provable guarantees.The former approach takes advantage of the random thresholds, while the latter predicts the direction and magnitude separately.The authors of [12] introduce an adaptive thresholding scheme which utilizes a generalized approximate message passing algorithm (GAMP) [12] for recovery and thresholds update throughout sampling.In a different approach, the work [13] proposes an adaptive quantization and recovery scenario making an exponential error decay in one-bit CS frameworks.The authors of [15] use an adaptive process to take measurements around the estimated signal in each iteration.While [15] only changes the bias (mean) of the estimated signal, our algorithm also generates random dithers with adaptive variance.The variance value of thresholds initializes from an overestimation of desired signal and is divided by a factor of 2 in each iteration.This adaptive thresholding scheme forces the algorithm estimate to concentrate around the optimal solution and makes the feasible set smaller in each iteration (random part with reduced variance).
Many of the techniques mentioned for adaptive sparse signal recovery do not generalize (at least not in an obvious strategy) to dictionary-sparse signal.For example, determining a surrogate of f that is supposed to be of lower complexity with respect to D H is non-trivial and challenging.We should emphasize that while the proofs and main parts in [13] rely on hard thresholding operator, it could not be used for either effective or exact dictionary-sparse signals.This is due to that given a vector x in the analysis domain 3 R N , one cannot guarantee the existence of a signal f in R n such that D H f = x .Recently, the work [8] shows both direction and magnitude of a dictionary-sparse signal can be recovered by a convex program with strong guarantees.The work [8] has inspired our work for recovering dictionary-sparse signal in an adaptive manner.In contrast to the existing method [8] for binary dictionary-sparse signal recovery which takes all of the measurements in one step with fixed settings, we solve the problem in an adaptive multistage way.In each stage, regarding the estimate from previous stage, our algorithm is propelled to the desired signal.In the non-adaptive work [8], the error rate is poorly large, while in our work, the error rate exponentially decays with increased number of adaptive steps.
Notation.Here, we introduce the notation used in the paper.Vectors and matrices are denoted by boldface lowercase and capital letters, respectively.A T and A H stand for transposition and Hermitian of A , respectively.C and c denote positive absolute con- stants which can be different from line to line.We use We write S n−1 := {v ∈ R n : �v� 2 = 1} for the unit Euclidean sphere in R n .For x ∈ R n , we define x S as the sub-vector in R |S| consisting of the entries indexed by the set S.

Method
In this section, we present the general procedure of our adaptive scheme in details.Our system model is built upon the optimization problem (4).A major part of this problem is to choose an efficient threshold ϕ ∈ R m .To this end, we propose a closed- loop feedback system (see Fig. 1) which exploits prior information from previous stages.The superscript (i) corresponding to an element represents the i-th stage of that element in the adaptive algorithm.Our adaptive approach consists of three algorithms high-dimensional threshold generator ( HDTG ), adaptive sampling ( AS ) and 1 Block diagram of adaptive sampling procedure 3 The domain in which the desired signal is being analyzed due to having a particular structure, e.g., sparsity.
adaptive recovery ( AR ).HDTG specifies the method of choosing the thresholds based on previous estimates.AS takes one-bit measurements of the signal by using HDTG and returns the one-bit samples and corresponding thresholds.Lastly, AR recovers the signal given the one-bit measurements.We provide a mathematical framework to guarantee our algorithm results in the following.
Theorem 1 (Main theorem) Let r, c ǫ , γ , ǫ > 0 .Consider f ∈ R n as the desired effective s-analysis-sparse signal with �f � 2 ≤ r and A ∈ R m×n is the measurement matrix with standard normal entries where m is the total number of measurements divided into T stages.Assume that be the sampling and recovery algorithms introduced later in Algorithms 2 and 3, respectively, where ϕ is determined by Algorithm 1.Then, with probability at least 1 − γ exp( −c ǫ m T ) over the choice of ϕ and A , the output of Algorithm 3 satisfies A remarkable note is that if we only consider one stage, i.e., T = 1 , the exponential behavior of our error bound disappears and reaches the state-of-the-art error bound [8,Theorem 8].In fact, the results of [8] are a special case of our work when the thresholds are non-adaptive.The term adaptivity is related to the threshold updating, and the measurement matrix ( A ) is fixed.
In what follows, we provide rigorous explanations about our proposed algorithms 1, 2 and 3 in three items.
1. High-Dimensional Threshold: Our algorithm for high-dimensional threshold selection is given in Algorithm 1.The algorithm output consists of two parts: deterministic (i.e., Af i ) and random (i.e., τ ∈ R q ) .The former transfers the origin to the pre- vious signal estimate ( f i ), while the latter creates measurement dithers from the origin.(The variance parameter σ 2 controls the threshold distance from f i .) (5) 2. Adaptive Sampling: Our adaptive sampling algorithm is given in Algorithm 2. To implement this algorithm, we need the dictionary D ∈ R n×N , the measurement matrix A ∈ R m×n , linear measurement Af and an over estimation of signal power r ( f 2 ≤ r ).At the first stage, we initialize signal estimation to zero vector.We choose T stages for our algorithm.At the i-th stage, the measurement matrix A (i) is taken from the i-th row subset (of size q := m T ) of A .The adaptive sampling pro- cess consists of four essential parts.First, in step 2 of the pseudocode, we generate the high-dimensional thresholds using Algorithm 1 by the parameters σ 2 = 2 1−i r and f i .Second, we compare the linear measurement block A (i) f with the thresh- old vector ϕ (i) and obtain the sample vector y (i) (step 3 of Algorithm 2).Third, we implement a second-order cone (steps 4 and 5) to find an approximate for f .However, this strategy does not often lead to an effective dictionary-sparse sig- nal.So, in step 6, we devise a strategy to find a low-complexity approximation of f with respect to operator D H .In other words, we project the resulting signal onto the nearest effective s-analysis-sparse signal.We refer to the third part of the adaptive sampling algorithm (steps 4-6) as single-step recovery (SSR).The estimated signal at each stage ( f i ) builds the deterministic part of our high-dimensional threshold in step 2. Finally, Algorithm 2 returns binary vectors {y (i) } T i=1 and the threshold vectors {ϕ (i) } T i=1 to the output.

Adaptive Recovery:
In the recovery procedure (Algorithm 3), we need the dictionary D , the measurement matrix A , binary measurements vector y , high-dimensional threshold vector ϕ and an upper norm estimation of signal r.In the adaptive recovery algorithm, we first divide the inputs ( y and A ) into T blocks (i.e., y (i) and A (i) ).Then, we simply implement SSR on each block.The output of SSR in the last stage is the final recovered signal (i.e., f T ).

Intuitive description of our algorithms
We provide here an example to clarify how our algorithm works.In this example, we consider a 2D point denoted by point DS (Fig. 2a) as the desired signal and depict the algorithm procedure for two consecutive iterations.At the first iteration of the algorithm, we do not have any estimation of the desired signal.Therefore, we set the first estimate to zero (equivalently [0, 0, 1] T in the 3D domain or point P in Fig. 2a).We gen- erate high-dimensional thresholds and take four measurements (step 2 of AS ).These measurements are equivalently four hyperplane crossing the center of the hemisphere (Fig. 2b, step 3 of AS ).In the next step, we implement the augmented random hyper- plane tessellation on the hemisphere and compute the first estimate of the desired signal project it to the signal domain (steps 3-6 of AS ).In this example, we consider the point B as the first estimate ( f 1 in first iteration of AS ).In the next iteration, we transfer the hemisphere to the adjacent point B and reduce the radius by a factor of two (Fig. 2c).The point B is the high-dimensional threshold of the second iteration of the algorithm.We perform our algorithm with four measurements.The main advantage of the reduction in the radius of the hemisphere (equivalently the variance argument of the threshold generator) is to get more precise measurements in the neighborhood of the desired signal.In the second stage of the algorithm, we take measurements (or equivalently hyperplanes) around the threshold, leading to locate measurements more concentrated around the desired signal.The intersection of hyperplanes and signal space in this Fig. 2 3D illustration of random hyperplane tessellation for the first example of intuitive description example creates a new feasible set border which is depicted by the red line in Fig. 2d.In this intuitive description, we tried to clarify the random hyperplane tessellation usage in our proof sketch. 4.
In the second example, we follow our algorithm steps in a 2D example and compare it with the non-adaptive algorithm.In Fig. 3a, we perform the non-adaptive algorithm for eight measurements.It is obvious that the solid lines ( L 5 − L 8 ) or equivalently measure- ments do not have any impact on the final result.However, if we get these measurements in two iterations (i.e., stages) of the adaptive algorithm, the second bunch of samples ( L 5 − L 8 in Fig. 3(b)) is more concentrated toward the desired signal ( f ).

Experimental results and discussion
In this section, we investigate the performance of our algorithm and compare it with the two previous one-bit dictionary-sparse recovery given by [8].The first algorithm solves linear programming optimization (LP) [8, Subsection 4.1], and the second algorithm solves a second-order cone programming (CP) optimization [8,Subsection 4.2].First, we construct a matrix where its columns are drawn randomly and independently from S n−1 .Then, the dictionary D ∈ R n×N ( N = 1000, n = 50 ) is defined as an orthonormal basis of this matrix.Then, the signal f is generated as f = Bc where B is a basis of null(D S ) and c is drawn from standard normal distribu- tion.Here, S is used to denote the complement of the support of D H f .The elements of A are chosen from i.i.d.standard normal distribution.We set r = 2�f � 2 , σ = r and the number of stages to T = 10 .Define the normalized reconstruction error as �f − f � 2 �f � 2 ( f denotes the estimated signal in each algorithm, and it is represented by f T ) espe- cially in our algorithm.The results in Fig. 4a are obtained by implementing Algorithms 2 and 3 100 times and taking the average of the normalized reconstruction error.As it is clear from Fig. 4a, LP algorithm outperforms CP by 2 dB on average.Our algorithm with few measurements behaves slightly weaker than others.The proposed approach requires a minimum number of measurements to act perfectly.The main reason for this behavior is our thresholding scheme.The variance of dithers is reduced in each iteration (i.e., stage).We also take the measurements around the previous estimate of the desired signal.On the contrary, the CP and LP algorithms only reduce the signal space uniformly.When the number of measurements is below the required bound, our first estimate gets away from the desired signal, so the adaptive algorithm tries to find it around a false estimate.However, it seems that there is a phase transition behavior in our algorithm when the number of measurements increases.In fact, after a certain number of measurements, our proposed algorithm substantially outperforms both LP and CP.Table 1 shows the average execution time (CPU time) of the three algorithms per the specified number of measurements.
In the second experiment, we examine the performance of our algorithm for multiple degrees of sparsity.Figure 4b shows the result for s = 20, 30, 40.(Other parameters are assumed similar to the first experiment.)As it is clear from the figure, by increasing the sparsity level, the accuracy is decreased.
In the third experiment, we consider the Shepp-Logan phantom image as the groundtruth signal.Since the required number of measurements in one-bit CS is significantly larger than that of the conventional CS (see, e.g., [13,Section 5]), we split the picture into multiple blocks of size 32 × 32 and process each block separately, merely due to compu- tational complexity reasons.This block processing allows us to implement the algorithm in multiple threads and reduce the execution time.We use a redundant wavelet dictionary and 10 5 Gaussian measurements to recover each block.We evaluate the reconstruc- tion quality of final result in terms of the peak signal-to-noise ratio (PSNR) given by  where X and X are the true and estimated images.As shown in Fig. 5, LP (Fig. 5b) and CP (Fig. 5c) algorithms clearly fail with poor performances, but as it is evident in Fig. 5d, the output of our algorithm is almost similar to the desired picture.

Conclusion
This paper proposed a one-bit compressed sensing algorithm that recovers dictionarysparse signals from binary measurements by adaptively exploiting the inherent information in the sampling procedure.This scheme helps to reduce the number of needed measurements substantially.In addition, our algorithm exhibits an exponential decaying behavior in the reconstruction error.The proof approach is based on geometric theories in the high-dimensional estimation.In this work, we used geometric intuition to explain our result, which also can be used in other areas of signal processing.We believe our analysis can be extended to the multi-bit setting.We used a fixed reduction pattern in the thresholds dithers throughout this work.We believe that this reduction can be chosen smartly by extracting the geometric features in each algorithm step.

Proof of Theorem 1 1 Proof
By induction law, we show that (8) PSNR X, X = 20 log 10 holds with high probability for any i ∈ {1, . . ., T } .Consider the first step, i.e., i = 1 in Algorithm 2. At this step, our initial estimate f 0 is equal to 0 .Thus, the output of Algo- rithm 1 only contains the random part of high-dimensional thresholds (step 2 of Algorithm 2).Then, we obtain f 1 by using steps 4-9 of Algorithm 2 (except that we assume σ = r in Algorithm 2 and q = m T ).As a result, to verify �f − f 1 � 2 ≤ ǫr , we use [8, The- orem 8] to prove our result: Theorem 2 [8, Theorem 8] Let ǫ, r, σ > 0 , let and let A ∈ R m×n be populated by independent standard normal random variables.Fur- thermore, let τ 1 , . . ., τ m be independent normal random variables with mean zero and variance σ 2 that are also independent from A .Then, with failure probability at most Since f is effective s-analysis-sparse with �f � 2 ≤ 1 , we have that Further, the output of Algorithm 2 at the (i − 1)-th stage sat- isfies, i.e., �D H f i−1 � 1 ≤ √ sr , and it then holds that the signal f  (12) ϕ (i) ← HDTG(A (i) , q, 2 1−i r, f i−1 ), (13) y (i) = sign A (i) f − ϕ (i) .( 14) in Theorem 2. As a result, by exploiting this theorem and some simplifications, we shall have that with probability at least 1 − γ exp(−c ′ ǫ q) where in the latter expression we used ǫ 1 = 1 4 and i is the output vector of the optimization problem in step 4 of Algorithm 2. c ′ ǫ is a certain constant dependent on ǫ .Now, suppose that (16) occurs.Consider After applying step 6 of Algorithm 2, we obtain f i that has the property Then, we have Finally, by using ( 16), (17), and the fact that �f i − f tmp � 2 ≤ �f − f tmp � 2 (step 6), the lat- ter equation becomes Since we consider T stage in Algorithms 2 and 3, we reach the error bound: (15)

Fig. 3
Fig. 3 This figure compares the sampling procedure of the non-adaptive algorithm (a) with our adaptive algorithm (b)

Fig. 4
Fig.4 Reconstruction of a dictionary-sparse signal x ∈ R 50 in the dictionary D ∈ R 50×1000 .Both figures show the reconstruction error averaged over 100 Monte Carlo simulations.Image a compares the performance of the algorithms in s = 20 .Image b examines the performance of our algorithm over different sparsity levels

Fig. 5
Fig. 5 Simulation results for reconstruction of the Shepp-Logan phantom image (phantom(256) in MATLAB) in which the picture split to 32 × 32 blocks and each algorithm get m = 10 5 measurements per block observed via y = sign(Af − τ ) is approximated by with error Now, suppose that the result (9) holds in the (i − 1)-th step, i.e., Consider i-th stage of Algorithm 2 where the high-dimensional thresholds and the measurements are obtained as By substituting (12) in (13), we reach:

Table 1
Average execution time (millisecond)