Learning Dynamic Compressive Sensing Models
Learning Dynamic Compressive Sensing Models Dongeun Lee and Jaesik Choi
arXiv:1502.04538v2 [cs.IT] 3 Jul 2015
School of Electrical and Computer Engineering Ulsan National Institute of Science and Technology (UNIST), Ulsan 689-798, Korea [email protected], [email protected]
Abstract Random sampling in compressive sensing (CS) enables the compression of large amounts of input signals in an efficient manner, which is useful for many applications. CS reconstructs the compressed signals exactly with overwhelming probability when incoming data can be sparsely represented with a fixed number of components, which is one of the drawbacks of CS frameworks because the signal sparsity in many dynamic systems changes over time. We present a new CS framework that handles signals without the fixed sparsity assumption by incorporating the distribution of signal sparsity. We show that the signal recovery success in our beta distribution modeling is more accurate than the success probability analysis in the CS framework. Alternatively, the success or failure of signal recovery can be relaxed, and the numbers of components included in signal recoveries can be represented with a probability distribution. We show this distribution is skewed to the right and naturally represented by the gamma distribution. Experimental results confirm the accuracy of our modeling in the context of dynamic signal sparsity. Keywords: Compressive sensing, random sampling, dynamic signal sparsity, sparse signal recovery.
1
Introduction
Continuous flows of big data are generated by many sources nowadays. Among these, resource limited devices occupy a significant portion. For these devices, sensing and transmitting massive data are important challenges, as they are concerned with saving resources. Compressive sensing (CS) [1, 2, 3, 4, 5, 6, 7] is a well suited choice for resource limited devices because it enables the sensing and compression of massive data without the complexity burden imposed by conventional schemes. Recent advances in CS reduce the complexity burden even further by random sampling, so that CS schemes have been successfully applied to broader application areas [8, 9, 10]. More specifically, CS reconstructs the exact signals from the compressed codes with overwhelming probability when incoming data can be sparsely represented (i.e., small numbers of components). In other words, the probability of recovery failure (CS failing to reconstruct the exact input signal) can be bounded when the compressed signals have enough length. However, most CS frameworks are built based on the assumption that incoming data can be sparsely represented by a fixed number of components. This assumption does not hold in practice with many dynamic systems where the numbers of components change over time. The assumption also implies that the reconstruction would fail when input signals have more components (denser) than the predefined threshold. This prevents deriving a tight probabilistic model which exploits 1
the number of components in signal recovery. Recently introduced dynamic CS frameworks [11, 12, 13, 14, 15] enable tracking the changes of signal support and amplitude. However, an amortized error analysis on the dynamic CS framework has not been provided by existing CS frameworks yet. In this paper, we present a flexible CS framework that handles signals without the fixed sparsity assumption by incorporating the distribution of signal sparsity. Our framework models the distribution of the (changing) numbers of components of an input signal in a dynamic system. By incorporating the beta distribution, our model presents the signal recovery success more accurately than the success probability analysis in the CS framework. Furthermore, we relax the concept of signal recovery success and present the number of components included in the signal recovery as a varying quantity, for which we propose gamma distribution modeling. Since our new model can handle signals generated from general dynamic systems, it will pave the way to advance in CS frameworks. The rest of this paper is organized as follows. Section 2 reviews the CS frameworks and random sampling. Section 3 provides a new CS framework which incorporates the distribution of dynamic signal sparsity, where we prove our beta distribution modeling is more accurate than existing CS framework. In Section 4, we consider compressible signals where the success or failure of signal recovery does not exist. Section 5 exhibits experimental results, followed by concluding remarks in Section 6.
2
Compressive Sensing and Random Sampling
Compressive sensing, or compressed sampling (CS), is an efficient signal processing framework that incorporates signal acquisition and compression simultaneously [8, 16]. If a signal can be represented by only a few (significant) components with or without the help of a sparsifying basis, CS allows it to be efficiently acquired with a number of samples that is far fewer than the signal dimension and of the same order as the number of components.
2.1
Compressing While Sensing
In CS, a signal is projected onto random vectors whose cardinality is far below the dimension of the signal. Consider a signal x ∈ RN is compactly represented with a sparsifying basis Ψ having just a few components: x = Ψs, where s ∈ RN is the vector of transformed coefficients with a few significant coefficients. Here, Ψ could be a basis that makes x sparse in a transform domain such as the DCT, wavelet transform domains, or even the canonical basis, i.e., the identity matrix I, if x is sparse itself without the help of transform. DefinitionPA signal x is called K-sparse if it is a linear combination of only K N basis vectors such that K i=1 sni ψni , where {n1 , . . ., nK } ⊂ {1, . . ., N }; sni is a coefficient in s; and ψni is a column of Ψ. In practice, some signals may not be exactly K-sparse; rather they can be closely approximated with K basis vectors ignoring many small coefficients close to zero. This type of signal is called compressible [8, 9]. CS projects x onto a random sensing basis Φ ∈ RM ×N as follows (M 0, where sK is the vector s with all but the largest K components set to 0. When an original signal is exactly K-sparse, then s = sK with M = O(K log N ) measurements, which implies that the recovery is exact, i.e., s? = s. The random sensing basis Φ have RIP if (1 − δ)ksk22 ≤ kΦΨsk22 ≤ (1 + δ)ksk22 for small δ ≥ 0, and this condition applies to all K-sparse s. 2 The two bases Φ and Ψ are incoherent when the rows of Φ cannot sparsely represent the columns of Ψ and vice versa. 3 Note that greedy methods are not always fast. 1
3
3
Learning Recovery Success
The success of signal reconstruction in compressive sensing (CS) is not deterministic. For instance, when we say an exact recovery of a K-sparse signal is achievable with overwhelming probability, it implies there is also the chance of recovery not being exact. Most existing CS literature assumes a sufficient number of measurements M so that an exact recovery is almost always achievable [16, 9], which is based on the assumption that the sparsity K is already known or does not exceed a certain bound. However, the signal sparsity in dynamic systems can change over time and an excessive number of measurements can waste resources such as network bandwidth and storage space. Therefore, we propose a new dynamic CS framework for the random sampling of CS and provide a new theoretical analysis on signal recovery.
3.1
Compressive Sensing Framework
In the random sampling of CS, the number of required measurements M = O(K log N ) can be detailed as follows [9]: M ≥ C · K ln(N ) ln(−1 ) (4) for some constant C > 0, where ∈ (0, 1) denotes the probability of an inexact recovery of the K-sparse signal. In particular, the signal recovery succeeds with a probability of at least 1 − if (4) holds. We can then express (4) with regard to the probability of failure , which is given by M IP(s? 6= s | M, N, K) := ≤ exp − . C · ln(N )K
(5)
Thus, the probability of failure (inexact recovery) IP(s? 6= s | M, N, K) is conditional upon M , N , and K. Since we are interested in the dynamic signal sparsity K, we model K as a random variable with M and N as fixed quantities. If we denote an arbitrary probability density function (pdf) of K as fK (k), we can marginalize over fK (k) and find the upper bound of failure probability as follows: Z Z M ? ? IP(s 6= s | M, N ) = IP(s 6= s | M, N, K) · fK (k) dk ≤ exp − fK (k) dk. (6) C · ln(N )K k k The upper bound in (6) can have an analytic solution depending on the form of fK (k). In particular, this is the case when K follows certain distributions such as the inverse Gaussian distribution and 4 Therefore, we can state that a signal recovery succeeds with a probability the gamma distribution. R of at least 1 − k exp(−M/(C · ln(N )K))fK (k)dk, given the distribution of signal sparsity fK (k).
3.2
Beta Distribution Modeling
Unfortunately, the probability of signal recovery failure given in (5) does not hold in practice because there is a discrepancy between the failure probabilities in the CS framework and actual random sampling, as will be explained further in Section 5.1. Thus we have to model the success or failure probability of signal recovery from a new perspective. As an illustrative example, suppose that we examine the success probability by generating many different signed spike (±1) vectors for each signal sparsity and then performing experiments 4
Since K ≥ 0, probability distributions supported on semi-infinite intervals, i.e., (0, ∞), are rational choices.
4
K=20
K=21
K=22
K=23
K=24
K=25
K=26
K=27
K=28
K=29
K=30
35 30
Density
25 20 15 10 5 0 0
0.1
0.2
0.3
0.4 0.5 0.6 Success probability
0.7
0.8
0.9
1
Figure 1: Histograms of success probability for various K’s. Success probability distribution for each K was obtained with 300 different random signed spike vectors for N = 512 and M = 100. A single success probability for each signed spike vector was calculated with 300 experiments. for each signed spike vector.5 Fig. 1 shows histograms of success probability for various signal sparsities. We model the success probability shown in Fig. 1 by the beta distribution with its parameters α and β depending on signal sparsity, i.e., IP(s? = s | M, N, K) ∼ Beta(αK , βK ). We use the beta distribution since it is suitable for modeling the success/failure probability in that the beta distribution is supported on (0, 1) and the conjugate prior for the Benoulli and the binomial distributions. Moreover, in Fig. 1, we can see that the mean of success probability distribution constantly decreases as K increases; the variance of success probability distribution keeps increasing till K = 25 and then keeps decreasing. The maximum-likelihood estimates of the parameters α’s and β’s of the beta distributions in Fig. 1 are found by the least-square method. In particular, we consider the quantity M/(ln(N )K) in (5) because the signal sparsity K should be considered together with the number of measurements M and the signal length N . Thus we model both α and β as functions of this quantity, which is given by M M α(M, N, K) = aα ln + bα and β(M, N, K) = aβ + bβ , (7) ln(N )K ln(N )K where aα > 0, bα > 0; and aβ < 0, bβ > 0 are model parameters for α(·) and β(·) respectively. The accuracy of resulting model functions is presented in Appendix A. Incorporating extreme cases where signal recovery always succeeds or always fails (e.g., K < 20 or K > 30 in Fig. 1), we introduce Kmin and Kmax to denote the minimum and the maximum signal sparsities respectively that yield stochastic probability, not a deterministic result which can be represented as a success or a failure. Then we can model the new pdf of signal recovery success using the mixture of the Dirac delta function and the beta distribution. 5
Detailed settings are explained in Section 5.1.
5
Definition Let IP(s? = s | M, N ) := Π. The pdf of Π given K is given by6 k < Kmin δ(π − 1) fΠ|K (π | k) := Beta(αK , βK ) Kmin ≤ k ≤ Kmax . δ(π) Kmax < k
(8)
Combining this definition with an arbitrary pdf fK (k) of the dynamic signal sparsity K, we can find the success probability distribution marginalized over fK (k) as follows: Z fΠ (π) = fΠ|K (π | k)fK (k) dk k Z Kmin Z Kmax Z ∞ = δ(π − 1)fK (k) dk + Beta(αK , βK ) · fK (k) dk + δ(π)fK (k) dk 0
Kmin
Kmax
Z
Kmax
= δ(π − 1)FK (Kmin ) + δ(π)(1 − FK (Kmax )) +
Beta(αK , βK ) · fK (k) dk,
(9)
Kmin
where FK (·) is the cumulative distribution function (CDF) of K. The two Dirac delta function terms in (9) can be interpreted as probability masses. Since R Kmax Kmin Beta(αK , βK ) · fK (k) dk does not have an analytic solution, we compute the values numerically.
3.3
Modeling Accuracy
Here, we present the main theoretical contribution: the recovery success probability modeled by our beta distribution is tighter than the lower bound of that in the existing CS framework. We show the failure probability in the CS framework is incapable of reflecting the actual failure probability of signal recovery. It is not only that the inequality IP(s? 6= s | M, N, K) ≤ exp(−M/(C · ln(N )K)) cannot provide tight probability of failure, but the inequality itself is inaccurate. The following theorem is proven in Appendix B. Theorem 1 The recovery success probability modeled by the beta distribution is tighter than the lower bound of that given by the CS framework.
4
Learning Recovery Quality
When a signal of interest is not exactly K-sparse, but compressible as discussed in Section 2.1, the recovery of signal in Section 2.3 can be treated from a different perspective [20]. In particular, the inequality (3) is considered differently. If an original signal is compressible, then the quality of a recovered signal is proportional to that of the K most significant pieces of information. We get progressively better results as we compute more measurements M , since M = O(K log N ) [16]. Therefore, Ψs? ∈ RN also makes progress on its quality as M increases.7 From this viewpoint, the success or failure of signal recovery no longer exists. Rather, we can view the number of components included in the signal recovery as a varying quantity. Specifically, if a signal recovery is about to fail with a given K, then K can be lowered to make the recovery 6 7
Beta(αK , βK ) here is used to denote the pdf of the beta distribution. The error bound follows (3) as well if Ψ is an orthogonal matrix, which is usually the case.
6
eventually succeed. Here the number of included components K varies for different recoveries and signals, as analogous to the success probability in Section 3.2 that can be calculated with different recoveries and varies for different signals. In this regard, (3) can be utilized to infer varying K’s over different recoveries and signals. Here our assumption is that the upper bound in (3) is tight such that we solve the following optimization problem: max K subject to ks? − sk2 ≤ C1 · ks − sK k1 . (10) In (10), C1 has to be determined, where the maximum signal sparsity Kmax introduced in Section 3.2 plays a key role to set the upper limit on how large K can be, since K > Kmax is not reasonable. In particular, we can generate a compressible signal si ∈ S such that ksi k1 = C`1 for all i, where S is the set containing many different signals; C`1 > 0 is some constant. For each si , we have a set Si? that contains many different recoveries s?ij . Then C1 can be found as follows: C1 =
min ks?ij − si k2 max ksi − sK k1 i
,
(11)
max denotes the compressible signal si with all but the largest Kmax components set to 0. where sK i
Varying K’s obtained through (10) can be represented by a pdf, which has been empirically shown to follow the gamma distribution [20]. We are interested in the shape of this pdf. The following proposition is proven in Appendix C. Proposition 1 The pdf of K, the number of components included in the signal recovery of a compressible signal, is skewed to the right, i.e., right tailed.
4.1
Error Analysis with Dynamic Signal Sparsity
Since the success or failure of signal recovery does not exist in this framework, we instead investigate the amount of error occurring during the recovery procedure in expected value sense. In particular, the best K-term approximation ks − sK k1 in (3) is known to be bounded as follows [21]: ks − sK k1 ≤
2G , K
(12)
where the constant G can be learned by the power-law decay such that each magnitude of components in s sorted in decreasing order is upper bounded by G/i2 . (i = 1, . . . , N is the sorted index.) Then we can analyze the `2 error E of signal recovery assuming fK (k) = Gamma(κK , θK ), which is given by Z 2G 2C1 G E = C1 · fK (k) dk = B(κK − 1, 1), (13) k θK k where B(·, ·) is the beta function [20]. Here the pdf fK (k) is employed to represent varying K’s. In this framework, the dynamic signal sparsity should be included using an alternative approach. Since the mean of K would vary depending on the changing sparsity of a compressible input signal, the parameters κK and θK would vary as well. Specifically, we empirically found that θK does not change much with the varying mean. Thus we assume θK > 0 is a constant and κK > 1 is a random variable. Substituting κK − 1 for lK , we have the following dynamic error Edyn : Z 2C1 G Edyn = B(lK , 1)fLK (lK ) dlK , (14) lK θK 7
where fLK (lK ) is the pdf of LK that represents the dynamic signal sparsity. Note that the dynamic error Edyn in (14) has analytic solutions when LK follows the inverse Gaussian distribution and the gamma distribution, which are ((2C1 G)/θK ) · (1/µL + 1/λL ) and ((2C1 G)/θK )/((κL − 1)θL ) respectively.8
5
Experimental Results
5.1
Recovery Success
In Section 3, we discussed the discrepancy between the failure probabilities in the CS framework and actual random sampling. In order to show this discrepancy, we artificially generated signed spikes ±1 at random locations in proportion to desired sparsities and densified these spikes using Ψ9 to perform the random sampling. For each signal sparsity K, the actual failure probability can be calculated for different recovery experiments. To this end, we adopted an optimization method (basis pursuit) to solve the optimization problem in (2) [22]. Specifically, the primal-dual algorithm based on the interior point method was employed to solve (2) [17]. Fig. 2a shows that the actual failure probability of signal recovery with varying signal sparsity does not follow the failure probability given in the CS framework. The failure probability in (5) cannot model the actual failure probability of signal recovery, regardless of the value chosen for constant C. This result confirms Lemma 1 and Corollary 2. Moreover, in Section 3.2 we modeled the new pdf of signal recovery success fΠ (π) in (9). We compared this new pdf with the upper bound of failure in (6), given a dynamic signal sparsity K. Specifically, we employed the inverse Gaussian distribution such that fK (k) = IG(30, 200). Fig. 2b exhibits the efficacy of our beta distribution modeling, where the lower bounds of success probability given in the CS framework fail to capture actual success probability in random sampling. This result confirms Theorem 1. We also employed real-world environmental data sets obtained from wireless sensor network deployments for comparison [23]: (i) humidity (%), (ii) luminosity (V), and (iii) temperature (◦ C). Random numbers representing the dynamic signal sparsity K were drawn from the inverse Gaussian distribution and we used this K to randomly choose components sorted in decreasing order; other components were set to zero. Fig. 3 displays the success probability of signal recovery follows the shape of Fig. 2b. Specifically, two probability masses are evident for all of three data types.
5.2
Recovery Quality
When a signal is compressible and not exactly K-sparse, every signal is basically dense. In Section 4, we regarded the number of components included in the signal recovery as a varying quantity. We are interested in the general shape of this quantity in distribution. In order to verify Proposition 1, we performed experiments using real data sets as well as artificially generated random signed spikes. We can identify that these distributions are skewed to the right and naturally represented by the gamma distributions, as presented in Appendix D. Furthermore, we analyzed the error with dynamic signal sparsity in Section 4.1. In order to show the efficacy of Edyn in (14), we compared the solutions of (14) with real data sets. In particular, we used the inverse Gaussian distribution such that fLK (lK ) = IG(30/θK − 1, 200) for 8 9
µL and λL are parameters for the inverse Gaussian distribution; κL and θL for the gamma distribution. We used DCT as the sparsifying basis Ψ throughout experiments.
8
1 0.9
0.45
0.35
0.6
0.3
0.5
0.2
0.3
0.15
0.2
0.1
0.1
0.05 15
0.1879
0.25
0.4
0
0.4254
0.4
0.7 Density
Failure probability ( )
0.8
0.5 Actual probability C=0.5 C=1 C=2
20
25 Signal sparsity K
30
35
0
40
Distribution by (9)
0
0.1
0.2
0.3
0.4 0.5 0.6 Success probability
(a)
0.7
0.8
0.9
1
(b)
Figure 2: (a) Comparison between actual failure probability and failure probabilities given in (5) with varying C’s. Actual failure probability of each signal sparsity K was obtained with 300 experiments for N = 512 and M = 100. (b) Comparison between our new success probability distribution in (9) and the lower bounds of success probability obtained by (6) with varying C’s. The inverse Gaussian distribution was used for fK (k). Two probability masses are shown by vertical arrows, where solid boxes atop the arrows denote their probabilities. Three vertical dashed/dotted lines represent the lower bounds by (6): C = 0.5 at 0.6781; C = 1 at 0.4450; and C = 2 at 0.2596. Here, Kmin = 20 and Kmax = 30. 700
500
600
700 600
400
500
400 300
300
Density
Density
Density
500
200
200
300 200
100
100 0 0
400
0.2
0.4 0.6 Success probability
(a) humidity
0.8
1
0 0
100 0.2
0.4 0.6 Success probability
(b) luminosity
0.8
1
0 0
0.2
0.4 0.6 Success probability
0.8
1
(c) temperature
Figure 3: Histograms of success probability for (a) humidity (%) data, (b) luminosity (V) data, and (c) temperature (◦ C) data. Histogram was obtained with 1,500 random number generations to choose different signals and 100 different experiments for each signal with N = 512 and M = 100. the calculation of (14). Random numbers representing the dynamic signal sparsity were also drawn from the inverse Gaussian distribution IG(30/θK , 200).10 These numbers were used to find closely matching signals in data sets. Then these matching signals were used for experiments.11 Fig. 4 shows that Edyn is a useful estimator for the `2 error of signal recovery. 10 11
Note that this setting makes IE[K] = κK θK = 30, the same configuration as in Section 5.1. Other numbers without matching signals were discarded.
9
300
80
350 300
250 60
250
150
Density
Density
Density
200 40
100
150 100
20 50 0 0
200
50 10
20 L2 error
(a) humidity
30
40
0 0
2
4 L2 error
(b) luminosity
6
8
0 0
10
20
30 L2 error
40
50
60
(c) temperature
Figure 4: Histograms of `2 error for (a) humidity (%) data, (b) luminosity (V) data, and (c) temperature (◦ C) data. Histogram was obtained with 1,500 random number generations for different dynamic signal sparsities of signals with N = 512 and M = 100. The vertical dashed lines represent Edyn at (a) 41.5698, (b) 3.8354, and (c) 33.3732.
6
Conclusion
We have presented a new dynamic CS framework that handles signals without the fixed sparsity assumption. The success probability of signal recovery in random sampling was investigated when the signal sparsity can vary. The success probability analysis in the existing CS framework was shown to be incapable of reflecting actual success probability by both theoretical analysis and experiments. On the contrary, our beta distribution modeling could closely reflect the actual success probability. We also considered signals that cannot be exactly represented with a sparse representation, where we could alternatively view the number of components included in the signal recovery as a varying quantity. This quantity was shown by both theoretical analysis and experiments to follow a right-tailed distribution such as the gamma distribution. We modeled the dynamic signal sparsity for these signals as well and provided an error analysis whose usefulness was demonstrated by experiments.
References [1] W. Bajwa, J. Haupt, A. Sayeed, and R. Nowak, “Compressive wireless sensing,” in Proc. IPSN, 2006, pp. 134–142. [2] S. Ji and L. Carin, “Bayesian compressive sensing and projection optimization,” in Proc. ICML, 2007, pp. 377–384. [3] M. W. Seeger and H. Nickisch, “Compressed sensing and bayesian experimental design,” in Proc. ICML, 2008, pp. 912–919. [4] C. Luo, F. Wu, J. Sun, and C. W. Chen, “Compressive data gathering for large-scale wireless sensor networks,” in Proc. MobiCom, 2009, pp. 145–156. [5] D. Hsu, S. Kakade, J. Langford, and T. Zhang, “Multi-label prediction via compressed sensing.” in Proc. NIPS, 2009, pp. 772–780.
10
[6] M. Lopes, “Estimating unknown sparsity in compressed sensing,” in Proc. ICML, 2013, pp. 217–225. [7] D. Malioutov and K. Varshney, “Exact rule learning via boolean compressed sensing,” in Proc. ICML, 2013, pp. 765–773. [8] R. G. Baraniuk, “Compressive sensing [lecture notes],” IEEE Signal Process. Mag., vol. 24, no. 4, pp. 118–121, Jul. 2007. [9] S. Foucart and H. Rauhut, A Mathematical Introduction to Compressive Sensing. 2013.
Springer,
[10] D. Lee and J. Choi, “Low complexity sensing for big spatio-temporal data,” in Proc. BigData, 2014, pp. 323–328. [11] D. Sejdinovic, C. Andrieu, and R. Piechocki, “Bayesian sequential compressed sensing in sparse dynamical systems,” in Proc. Allerton, 2010, pp. 1730–1736. [12] B. Shahrasbi, A. Talari, and N. Rahnavard, “TC-CSBP: Compressive sensing for timecorrelated data based on belief propagation,” in Proc. CISS, 2011, pp. 1–6. [13] J. Ziniel and P. Schniter, “Dynamic compressive sensing of time-varying signals via approximate message passing,” IEEE Trans. Signal Process., vol. 61, no. 21, pp. 5270–5284, Nov. 2013. [14] S. Ganguli and H. Sompolinsky, “Short-term memory in neuronal networks through dynamical compressed sensing,” in Proc. NIPS, 2010, pp. 667–675. [15] D. M. Malioutov, S. R. Sanghavi, and A. S. Willsky, “Sequential compressed sensing,” IEEE J. Sel. Top. Signal Process., vol. 4, no. 2, pp. 435–444, Apr. 2010. [16] E. J. Cand`es and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag., vol. 25, no. 2, pp. 21–30, Mar. 2008. [17] S. Boyd and L. Vandenberghe, Convex Optimization.
Cambridge University Press, 2004.
[18] Y. C. Pati, R. Rezaiifar, and P. S. Krishnaprasad, “Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition,” in Proc. ACSSC, 1993, pp. 40–44. [19] T. Blumensath and M. E. Davies, “Iterative thresholding for sparse approximations,” J. Fourier Anal. Appl., vol. 14, no. 5-6, pp. 629–654, Dec. 2008. [20] D. Lee and J. Choi, “Learning compressive sensing models for big spatio-temporal data,” in Proc. SDM, 2015, pp. 667–675. [21] R. G. Baraniuk, V. Cevher, M. F. Duarte, and C. Hegde, “Model-based compressive sensing,” IEEE Trans. Inf. Theory, vol. 56, no. 4, pp. 1982–2001, Apr. 2010. [22] S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition by basis pursuit,” SIAM J. Sci. Comput., vol. 20, no. 1, pp. 33–61, Jan. 1998. [23] G. Quer, R. Masiero, G. Pillonetto, M. Rossi, and M. Zorzi, “Sensing, compression, and recovery for WSNs: Sparse signal modeling and monitoring framework,” IEEE Trans. Wireless Commun., vol. 11, no. 10, pp. 3447–3461, Oct. 2012. 11
Appendix A The model functions in (7) are plotted in Fig. 5 with α’s and β’s obtained by the maximumlikelihood estimation for comparison. 250
250
200
200
150
Parameter Model estimation by (7) Parameter Model estimation by (7)
100
50
0 0.5
150
100
50
0.55
0.6
0.65 0.7 M/(ln(N)K)
0.75
0.8
0 0.85
Figure 5: Two beta distribution parameters as functions of M , N , and K estimated by (7). Fig. 6 also presents (sample) means and variances of the success probability distributions shown in Fig. 1, along with means and variances calculated using the model functions in (7). We can see that the beta distribution modeling with parameters in (7) effectively follows the actual success probability distributions. −3
Mean
0.5
0.5
Variance
x 10 1
1
Actual mean Mean estimation Actual variance Variance estimation 0 0.5
0.55
0.6
0.65 0.7 M/(ln(N)K)
0.75
0.8
0 0.85
Figure 6: Mean and variance of success probability distributions in Fig. 1 estimated with the beta distribution modeling α(·) and β(·).
Appendix B: Proof of Theorem 1 The inaccuracy of the failure probability in the CS framework results from the slowly decaying lower bound of success probability, that is, 1 − exp(−M/(C · ln(N )K)). In fact, we can show this lower bound decays slower than a power-law decay by the following lemma.
12
Lemma 1 There exists K0 > 0 such that for all K > K0 , the lower bound of recovery success probability is greater than the value of a power-law-decay function. Proof We need to show the following inequality M 1 − exp − > K −α C · ln(N )K
(15)
holds if K > K0 for some K0 > 0, where α > 0. Adding, subtracting, and taking the power K on both sides yields M −α K (1 − K ) > exp − . (16) C · ln(N ) We now use the binomial approximation on the left-hand side: (1 − K −α )K ≥ 1 − K · K −α . Thus we instead prove the following inequality M 1 − K 1−α > exp − . (17) C · ln(N ) holds if K > K0 for some K0 > 0. If we assume α > 1, then adding, subtracting, and taking the power 1/(1 − α) on both sides of (17) yields 1/(1−α) M 1 − exp − < K. (18) C · ln(N ) Setting K0 = (1 − exp(−M/(C · ln(N ))))1/(1−α) , we can argue that for all K > K0 , the lower bound of recovery success probability is greater than the value of a power-law-decay function. Corollary 2 In the CS framework, there is always a chance of succeeding at signal recovery however large K is. Proof The power-law-decay function K −α in (15) slowly converges to zero as K → ∞: its value is noticeably greater than zero even with large K. As the lower bound of recovery success probability is greater than the value of the power-law-decay function for all K > K0 , we can say there is always a chance of recovery success however large K is. We can now show that our beta distribution modeling provides more accurate success probability. The claim of the CS framework in Corollary 2 is in fact implausible because it says we can even set K > M and there is still a chance of success. We cannot expect signal recovery with a number of measurements M less than K. On the contrary, our beta distribution modeling can yield IP(s? = s | M, N, K) = 0 with a bounded Kmax . In particular, we show that the mean of Beta(αK , βK ), αK /(αK + βK ), converges to zero with a bounded Kmax . If we state αK /(αK + βK ) = 0 using (7), we obtain M aα ln + bα = 0. (19) ln(N )K Representing (19) with respect to K, we have Kmax
M = exp ln(N ) 13
bα aα
< ∞.
(20)
Similarly, we show this mean converges to one (IP(s? = s | M, N, K) = 1) with Kmin which is not so close to zero, whereas the lower bound of the recovery success probability given by the CS framework converges to one only if K is very close to zero. If we state αK /(αK + βK ) = 1 using (7), we have aβ
M + bβ = 0. ln(N )K
(21)
Representing (21) with respect to K, we have Kmin = −(aβ M )/(bβ ln(N )). Plugging Kmin into 1 − exp(−M/(C · ln(N )K)) shows that bβ 1 − exp > 0, (22) C · aβ as C > 0, aβ < 0, and bβ > 0: that is, 1 − exp(−M/(C · ln(N )K)) → 1 if and only if K → 0. Since 0 < Kmin < Kmax < ∞, we can argue that the beta distribution modeling can provide tighter recovery success probability.
Appendix C: Proof of Proposition 1 Since ksi k1 = C`1 for all i, we can conceive the same sequence {sn } of elements (absolute values) in si for all i. Then we have NX −K K ksi − si k1 = sn . (23) n=1
P −K Without loss of generality, we consider the partial sum N n=1 sn in (23) to be an arithmetic series that can be represented by a quadratic function in terms of K. We also assume the inequality constraint in (10) is the equality constraint such that ks? − sk2 = C1 · ks − sK k1 . q If we take the (partial) inverse function of the quadratic function, we have K ∼ Kmax − ks? − sk2 − (min ks?ij − si k2 ). Assuming the distribution of ks? − sk2 is symmetric (zero skewness), this asymptotic relation says ks? − sk2 will be compressed as it becomes large, which in turn makes the pdf of K right tailed. P −K A similar claim can be made if we consider the partial sum N n=1 sn to be a geometric series, ? where K ∼ N − log(ks − sk2 ). In this case, the pdf of K is skewed to the right as well.
Appendix D First, random signed spikes were artificially generated in different magnitudes at random locations and densified to perform random sampling. In particular, we considered an arithmetic sequence of length 50 (2, 4, 6, . . ., 98, 100), whose elements were placed at random locations in each vector. These signals are dense enough to be used for experiments because signal recovery always fails when K > 30 in our case, as shown in Fig. 2a. Fig. 7 displays the histogram of K, the number of components included in each signal recovery, which was obtained using the method explained in Section 4. We can identify that Proposition 1 actually holds here, as this distribution is skewed to the right. Furthermore, we empirically found that the distribution follows the gamma distribution, which is also natural since the gamma distribution has positive skewness, i.e., right tailed. 14
Actual distribution Gamma fitting
0.25
Density
0.2
0.15
0.1
0.05
0
20
22
24
26 K
28
30
32
Figure 7: Distribution of K fitted with a gamma distribution Gamma(242.81, 0.09), using the maximum likelihood estimation. Histogram was obtained with 300 different signals and 300 different experiments for each signal, for N = 512 and M = 100. We also provide results with real-world sensor data sets to verify Proposition 1 and the gamma distribution fitting. Figs. 8, 9, and 10 display the histograms of K and their gamma distribution fittings, where we can see that K follows the gamma distribution as well. Actual distribution Gamma fitting
0.09 0.08 0.07
Density
0.06 0.05 0.04 0.03 0.02 0.01 0
5
10
15
K
20
25
30
Figure 8: Distribution of K fitted with a gamma distribution Gamma(8.43, 1.71) for humidity (%) data. Histogram was obtained with 1,000 experiments over the same signal for N = 512 and M = 100.
15
0.04 Actual distribution Gamma fitting
0.035 0.03
Density
0.025 0.02 0.015 0.01 0.005 0
0
10
20
30 K
40
50
60
Figure 9: Distribution of K fitted with a gamma distribution Gamma(2.55, 7.94) for luminosity (V) data. Histogram was obtained with 1,000 experiments over the same signal for N = 512 and M = 100.
Actual distribution Gamma fitting
0.06
Density
0.05 0.04 0.03 0.02 0.01 0
0
5
10
15
K
20
25
30
35
Figure 10: Distribution of K fitted with a gamma distribution Gamma(5.75, 2.73) for temperature (◦ C) data. Histogram was obtained with 1,000 experiments over the same signal for N = 512 and M = 100.
16
arXiv:1502.04538v2 [cs.IT] 3 Jul 2015
School of Electrical and Computer Engineering Ulsan National Institute of Science and Technology (UNIST), Ulsan 689-798, Korea [email protected], [email protected]
Abstract Random sampling in compressive sensing (CS) enables the compression of large amounts of input signals in an efficient manner, which is useful for many applications. CS reconstructs the compressed signals exactly with overwhelming probability when incoming data can be sparsely represented with a fixed number of components, which is one of the drawbacks of CS frameworks because the signal sparsity in many dynamic systems changes over time. We present a new CS framework that handles signals without the fixed sparsity assumption by incorporating the distribution of signal sparsity. We show that the signal recovery success in our beta distribution modeling is more accurate than the success probability analysis in the CS framework. Alternatively, the success or failure of signal recovery can be relaxed, and the numbers of components included in signal recoveries can be represented with a probability distribution. We show this distribution is skewed to the right and naturally represented by the gamma distribution. Experimental results confirm the accuracy of our modeling in the context of dynamic signal sparsity. Keywords: Compressive sensing, random sampling, dynamic signal sparsity, sparse signal recovery.
1
Introduction
Continuous flows of big data are generated by many sources nowadays. Among these, resource limited devices occupy a significant portion. For these devices, sensing and transmitting massive data are important challenges, as they are concerned with saving resources. Compressive sensing (CS) [1, 2, 3, 4, 5, 6, 7] is a well suited choice for resource limited devices because it enables the sensing and compression of massive data without the complexity burden imposed by conventional schemes. Recent advances in CS reduce the complexity burden even further by random sampling, so that CS schemes have been successfully applied to broader application areas [8, 9, 10]. More specifically, CS reconstructs the exact signals from the compressed codes with overwhelming probability when incoming data can be sparsely represented (i.e., small numbers of components). In other words, the probability of recovery failure (CS failing to reconstruct the exact input signal) can be bounded when the compressed signals have enough length. However, most CS frameworks are built based on the assumption that incoming data can be sparsely represented by a fixed number of components. This assumption does not hold in practice with many dynamic systems where the numbers of components change over time. The assumption also implies that the reconstruction would fail when input signals have more components (denser) than the predefined threshold. This prevents deriving a tight probabilistic model which exploits 1
the number of components in signal recovery. Recently introduced dynamic CS frameworks [11, 12, 13, 14, 15] enable tracking the changes of signal support and amplitude. However, an amortized error analysis on the dynamic CS framework has not been provided by existing CS frameworks yet. In this paper, we present a flexible CS framework that handles signals without the fixed sparsity assumption by incorporating the distribution of signal sparsity. Our framework models the distribution of the (changing) numbers of components of an input signal in a dynamic system. By incorporating the beta distribution, our model presents the signal recovery success more accurately than the success probability analysis in the CS framework. Furthermore, we relax the concept of signal recovery success and present the number of components included in the signal recovery as a varying quantity, for which we propose gamma distribution modeling. Since our new model can handle signals generated from general dynamic systems, it will pave the way to advance in CS frameworks. The rest of this paper is organized as follows. Section 2 reviews the CS frameworks and random sampling. Section 3 provides a new CS framework which incorporates the distribution of dynamic signal sparsity, where we prove our beta distribution modeling is more accurate than existing CS framework. In Section 4, we consider compressible signals where the success or failure of signal recovery does not exist. Section 5 exhibits experimental results, followed by concluding remarks in Section 6.
2
Compressive Sensing and Random Sampling
Compressive sensing, or compressed sampling (CS), is an efficient signal processing framework that incorporates signal acquisition and compression simultaneously [8, 16]. If a signal can be represented by only a few (significant) components with or without the help of a sparsifying basis, CS allows it to be efficiently acquired with a number of samples that is far fewer than the signal dimension and of the same order as the number of components.
2.1
Compressing While Sensing
In CS, a signal is projected onto random vectors whose cardinality is far below the dimension of the signal. Consider a signal x ∈ RN is compactly represented with a sparsifying basis Ψ having just a few components: x = Ψs, where s ∈ RN is the vector of transformed coefficients with a few significant coefficients. Here, Ψ could be a basis that makes x sparse in a transform domain such as the DCT, wavelet transform domains, or even the canonical basis, i.e., the identity matrix I, if x is sparse itself without the help of transform. DefinitionPA signal x is called K-sparse if it is a linear combination of only K N basis vectors such that K i=1 sni ψni , where {n1 , . . ., nK } ⊂ {1, . . ., N }; sni is a coefficient in s; and ψni is a column of Ψ. In practice, some signals may not be exactly K-sparse; rather they can be closely approximated with K basis vectors ignoring many small coefficients close to zero. This type of signal is called compressible [8, 9]. CS projects x onto a random sensing basis Φ ∈ RM ×N as follows (M 0, where sK is the vector s with all but the largest K components set to 0. When an original signal is exactly K-sparse, then s = sK with M = O(K log N ) measurements, which implies that the recovery is exact, i.e., s? = s. The random sensing basis Φ have RIP if (1 − δ)ksk22 ≤ kΦΨsk22 ≤ (1 + δ)ksk22 for small δ ≥ 0, and this condition applies to all K-sparse s. 2 The two bases Φ and Ψ are incoherent when the rows of Φ cannot sparsely represent the columns of Ψ and vice versa. 3 Note that greedy methods are not always fast. 1
3
3
Learning Recovery Success
The success of signal reconstruction in compressive sensing (CS) is not deterministic. For instance, when we say an exact recovery of a K-sparse signal is achievable with overwhelming probability, it implies there is also the chance of recovery not being exact. Most existing CS literature assumes a sufficient number of measurements M so that an exact recovery is almost always achievable [16, 9], which is based on the assumption that the sparsity K is already known or does not exceed a certain bound. However, the signal sparsity in dynamic systems can change over time and an excessive number of measurements can waste resources such as network bandwidth and storage space. Therefore, we propose a new dynamic CS framework for the random sampling of CS and provide a new theoretical analysis on signal recovery.
3.1
Compressive Sensing Framework
In the random sampling of CS, the number of required measurements M = O(K log N ) can be detailed as follows [9]: M ≥ C · K ln(N ) ln(−1 ) (4) for some constant C > 0, where ∈ (0, 1) denotes the probability of an inexact recovery of the K-sparse signal. In particular, the signal recovery succeeds with a probability of at least 1 − if (4) holds. We can then express (4) with regard to the probability of failure , which is given by M IP(s? 6= s | M, N, K) := ≤ exp − . C · ln(N )K
(5)
Thus, the probability of failure (inexact recovery) IP(s? 6= s | M, N, K) is conditional upon M , N , and K. Since we are interested in the dynamic signal sparsity K, we model K as a random variable with M and N as fixed quantities. If we denote an arbitrary probability density function (pdf) of K as fK (k), we can marginalize over fK (k) and find the upper bound of failure probability as follows: Z Z M ? ? IP(s 6= s | M, N ) = IP(s 6= s | M, N, K) · fK (k) dk ≤ exp − fK (k) dk. (6) C · ln(N )K k k The upper bound in (6) can have an analytic solution depending on the form of fK (k). In particular, this is the case when K follows certain distributions such as the inverse Gaussian distribution and 4 Therefore, we can state that a signal recovery succeeds with a probability the gamma distribution. R of at least 1 − k exp(−M/(C · ln(N )K))fK (k)dk, given the distribution of signal sparsity fK (k).
3.2
Beta Distribution Modeling
Unfortunately, the probability of signal recovery failure given in (5) does not hold in practice because there is a discrepancy between the failure probabilities in the CS framework and actual random sampling, as will be explained further in Section 5.1. Thus we have to model the success or failure probability of signal recovery from a new perspective. As an illustrative example, suppose that we examine the success probability by generating many different signed spike (±1) vectors for each signal sparsity and then performing experiments 4
Since K ≥ 0, probability distributions supported on semi-infinite intervals, i.e., (0, ∞), are rational choices.
4
K=20
K=21
K=22
K=23
K=24
K=25
K=26
K=27
K=28
K=29
K=30
35 30
Density
25 20 15 10 5 0 0
0.1
0.2
0.3
0.4 0.5 0.6 Success probability
0.7
0.8
0.9
1
Figure 1: Histograms of success probability for various K’s. Success probability distribution for each K was obtained with 300 different random signed spike vectors for N = 512 and M = 100. A single success probability for each signed spike vector was calculated with 300 experiments. for each signed spike vector.5 Fig. 1 shows histograms of success probability for various signal sparsities. We model the success probability shown in Fig. 1 by the beta distribution with its parameters α and β depending on signal sparsity, i.e., IP(s? = s | M, N, K) ∼ Beta(αK , βK ). We use the beta distribution since it is suitable for modeling the success/failure probability in that the beta distribution is supported on (0, 1) and the conjugate prior for the Benoulli and the binomial distributions. Moreover, in Fig. 1, we can see that the mean of success probability distribution constantly decreases as K increases; the variance of success probability distribution keeps increasing till K = 25 and then keeps decreasing. The maximum-likelihood estimates of the parameters α’s and β’s of the beta distributions in Fig. 1 are found by the least-square method. In particular, we consider the quantity M/(ln(N )K) in (5) because the signal sparsity K should be considered together with the number of measurements M and the signal length N . Thus we model both α and β as functions of this quantity, which is given by M M α(M, N, K) = aα ln + bα and β(M, N, K) = aβ + bβ , (7) ln(N )K ln(N )K where aα > 0, bα > 0; and aβ < 0, bβ > 0 are model parameters for α(·) and β(·) respectively. The accuracy of resulting model functions is presented in Appendix A. Incorporating extreme cases where signal recovery always succeeds or always fails (e.g., K < 20 or K > 30 in Fig. 1), we introduce Kmin and Kmax to denote the minimum and the maximum signal sparsities respectively that yield stochastic probability, not a deterministic result which can be represented as a success or a failure. Then we can model the new pdf of signal recovery success using the mixture of the Dirac delta function and the beta distribution. 5
Detailed settings are explained in Section 5.1.
5
Definition Let IP(s? = s | M, N ) := Π. The pdf of Π given K is given by6 k < Kmin δ(π − 1) fΠ|K (π | k) := Beta(αK , βK ) Kmin ≤ k ≤ Kmax . δ(π) Kmax < k
(8)
Combining this definition with an arbitrary pdf fK (k) of the dynamic signal sparsity K, we can find the success probability distribution marginalized over fK (k) as follows: Z fΠ (π) = fΠ|K (π | k)fK (k) dk k Z Kmin Z Kmax Z ∞ = δ(π − 1)fK (k) dk + Beta(αK , βK ) · fK (k) dk + δ(π)fK (k) dk 0
Kmin
Kmax
Z
Kmax
= δ(π − 1)FK (Kmin ) + δ(π)(1 − FK (Kmax )) +
Beta(αK , βK ) · fK (k) dk,
(9)
Kmin
where FK (·) is the cumulative distribution function (CDF) of K. The two Dirac delta function terms in (9) can be interpreted as probability masses. Since R Kmax Kmin Beta(αK , βK ) · fK (k) dk does not have an analytic solution, we compute the values numerically.
3.3
Modeling Accuracy
Here, we present the main theoretical contribution: the recovery success probability modeled by our beta distribution is tighter than the lower bound of that in the existing CS framework. We show the failure probability in the CS framework is incapable of reflecting the actual failure probability of signal recovery. It is not only that the inequality IP(s? 6= s | M, N, K) ≤ exp(−M/(C · ln(N )K)) cannot provide tight probability of failure, but the inequality itself is inaccurate. The following theorem is proven in Appendix B. Theorem 1 The recovery success probability modeled by the beta distribution is tighter than the lower bound of that given by the CS framework.
4
Learning Recovery Quality
When a signal of interest is not exactly K-sparse, but compressible as discussed in Section 2.1, the recovery of signal in Section 2.3 can be treated from a different perspective [20]. In particular, the inequality (3) is considered differently. If an original signal is compressible, then the quality of a recovered signal is proportional to that of the K most significant pieces of information. We get progressively better results as we compute more measurements M , since M = O(K log N ) [16]. Therefore, Ψs? ∈ RN also makes progress on its quality as M increases.7 From this viewpoint, the success or failure of signal recovery no longer exists. Rather, we can view the number of components included in the signal recovery as a varying quantity. Specifically, if a signal recovery is about to fail with a given K, then K can be lowered to make the recovery 6 7
Beta(αK , βK ) here is used to denote the pdf of the beta distribution. The error bound follows (3) as well if Ψ is an orthogonal matrix, which is usually the case.
6
eventually succeed. Here the number of included components K varies for different recoveries and signals, as analogous to the success probability in Section 3.2 that can be calculated with different recoveries and varies for different signals. In this regard, (3) can be utilized to infer varying K’s over different recoveries and signals. Here our assumption is that the upper bound in (3) is tight such that we solve the following optimization problem: max K subject to ks? − sk2 ≤ C1 · ks − sK k1 . (10) In (10), C1 has to be determined, where the maximum signal sparsity Kmax introduced in Section 3.2 plays a key role to set the upper limit on how large K can be, since K > Kmax is not reasonable. In particular, we can generate a compressible signal si ∈ S such that ksi k1 = C`1 for all i, where S is the set containing many different signals; C`1 > 0 is some constant. For each si , we have a set Si? that contains many different recoveries s?ij . Then C1 can be found as follows: C1 =
min ks?ij − si k2 max ksi − sK k1 i
,
(11)
max denotes the compressible signal si with all but the largest Kmax components set to 0. where sK i
Varying K’s obtained through (10) can be represented by a pdf, which has been empirically shown to follow the gamma distribution [20]. We are interested in the shape of this pdf. The following proposition is proven in Appendix C. Proposition 1 The pdf of K, the number of components included in the signal recovery of a compressible signal, is skewed to the right, i.e., right tailed.
4.1
Error Analysis with Dynamic Signal Sparsity
Since the success or failure of signal recovery does not exist in this framework, we instead investigate the amount of error occurring during the recovery procedure in expected value sense. In particular, the best K-term approximation ks − sK k1 in (3) is known to be bounded as follows [21]: ks − sK k1 ≤
2G , K
(12)
where the constant G can be learned by the power-law decay such that each magnitude of components in s sorted in decreasing order is upper bounded by G/i2 . (i = 1, . . . , N is the sorted index.) Then we can analyze the `2 error E of signal recovery assuming fK (k) = Gamma(κK , θK ), which is given by Z 2G 2C1 G E = C1 · fK (k) dk = B(κK − 1, 1), (13) k θK k where B(·, ·) is the beta function [20]. Here the pdf fK (k) is employed to represent varying K’s. In this framework, the dynamic signal sparsity should be included using an alternative approach. Since the mean of K would vary depending on the changing sparsity of a compressible input signal, the parameters κK and θK would vary as well. Specifically, we empirically found that θK does not change much with the varying mean. Thus we assume θK > 0 is a constant and κK > 1 is a random variable. Substituting κK − 1 for lK , we have the following dynamic error Edyn : Z 2C1 G Edyn = B(lK , 1)fLK (lK ) dlK , (14) lK θK 7
where fLK (lK ) is the pdf of LK that represents the dynamic signal sparsity. Note that the dynamic error Edyn in (14) has analytic solutions when LK follows the inverse Gaussian distribution and the gamma distribution, which are ((2C1 G)/θK ) · (1/µL + 1/λL ) and ((2C1 G)/θK )/((κL − 1)θL ) respectively.8
5
Experimental Results
5.1
Recovery Success
In Section 3, we discussed the discrepancy between the failure probabilities in the CS framework and actual random sampling. In order to show this discrepancy, we artificially generated signed spikes ±1 at random locations in proportion to desired sparsities and densified these spikes using Ψ9 to perform the random sampling. For each signal sparsity K, the actual failure probability can be calculated for different recovery experiments. To this end, we adopted an optimization method (basis pursuit) to solve the optimization problem in (2) [22]. Specifically, the primal-dual algorithm based on the interior point method was employed to solve (2) [17]. Fig. 2a shows that the actual failure probability of signal recovery with varying signal sparsity does not follow the failure probability given in the CS framework. The failure probability in (5) cannot model the actual failure probability of signal recovery, regardless of the value chosen for constant C. This result confirms Lemma 1 and Corollary 2. Moreover, in Section 3.2 we modeled the new pdf of signal recovery success fΠ (π) in (9). We compared this new pdf with the upper bound of failure in (6), given a dynamic signal sparsity K. Specifically, we employed the inverse Gaussian distribution such that fK (k) = IG(30, 200). Fig. 2b exhibits the efficacy of our beta distribution modeling, where the lower bounds of success probability given in the CS framework fail to capture actual success probability in random sampling. This result confirms Theorem 1. We also employed real-world environmental data sets obtained from wireless sensor network deployments for comparison [23]: (i) humidity (%), (ii) luminosity (V), and (iii) temperature (◦ C). Random numbers representing the dynamic signal sparsity K were drawn from the inverse Gaussian distribution and we used this K to randomly choose components sorted in decreasing order; other components were set to zero. Fig. 3 displays the success probability of signal recovery follows the shape of Fig. 2b. Specifically, two probability masses are evident for all of three data types.
5.2
Recovery Quality
When a signal is compressible and not exactly K-sparse, every signal is basically dense. In Section 4, we regarded the number of components included in the signal recovery as a varying quantity. We are interested in the general shape of this quantity in distribution. In order to verify Proposition 1, we performed experiments using real data sets as well as artificially generated random signed spikes. We can identify that these distributions are skewed to the right and naturally represented by the gamma distributions, as presented in Appendix D. Furthermore, we analyzed the error with dynamic signal sparsity in Section 4.1. In order to show the efficacy of Edyn in (14), we compared the solutions of (14) with real data sets. In particular, we used the inverse Gaussian distribution such that fLK (lK ) = IG(30/θK − 1, 200) for 8 9
µL and λL are parameters for the inverse Gaussian distribution; κL and θL for the gamma distribution. We used DCT as the sparsifying basis Ψ throughout experiments.
8
1 0.9
0.45
0.35
0.6
0.3
0.5
0.2
0.3
0.15
0.2
0.1
0.1
0.05 15
0.1879
0.25
0.4
0
0.4254
0.4
0.7 Density
Failure probability ( )
0.8
0.5 Actual probability C=0.5 C=1 C=2
20
25 Signal sparsity K
30
35
0
40
Distribution by (9)
0
0.1
0.2
0.3
0.4 0.5 0.6 Success probability
(a)
0.7
0.8
0.9
1
(b)
Figure 2: (a) Comparison between actual failure probability and failure probabilities given in (5) with varying C’s. Actual failure probability of each signal sparsity K was obtained with 300 experiments for N = 512 and M = 100. (b) Comparison between our new success probability distribution in (9) and the lower bounds of success probability obtained by (6) with varying C’s. The inverse Gaussian distribution was used for fK (k). Two probability masses are shown by vertical arrows, where solid boxes atop the arrows denote their probabilities. Three vertical dashed/dotted lines represent the lower bounds by (6): C = 0.5 at 0.6781; C = 1 at 0.4450; and C = 2 at 0.2596. Here, Kmin = 20 and Kmax = 30. 700
500
600
700 600
400
500
400 300
300
Density
Density
Density
500
200
200
300 200
100
100 0 0
400
0.2
0.4 0.6 Success probability
(a) humidity
0.8
1
0 0
100 0.2
0.4 0.6 Success probability
(b) luminosity
0.8
1
0 0
0.2
0.4 0.6 Success probability
0.8
1
(c) temperature
Figure 3: Histograms of success probability for (a) humidity (%) data, (b) luminosity (V) data, and (c) temperature (◦ C) data. Histogram was obtained with 1,500 random number generations to choose different signals and 100 different experiments for each signal with N = 512 and M = 100. the calculation of (14). Random numbers representing the dynamic signal sparsity were also drawn from the inverse Gaussian distribution IG(30/θK , 200).10 These numbers were used to find closely matching signals in data sets. Then these matching signals were used for experiments.11 Fig. 4 shows that Edyn is a useful estimator for the `2 error of signal recovery. 10 11
Note that this setting makes IE[K] = κK θK = 30, the same configuration as in Section 5.1. Other numbers without matching signals were discarded.
9
300
80
350 300
250 60
250
150
Density
Density
Density
200 40
100
150 100
20 50 0 0
200
50 10
20 L2 error
(a) humidity
30
40
0 0
2
4 L2 error
(b) luminosity
6
8
0 0
10
20
30 L2 error
40
50
60
(c) temperature
Figure 4: Histograms of `2 error for (a) humidity (%) data, (b) luminosity (V) data, and (c) temperature (◦ C) data. Histogram was obtained with 1,500 random number generations for different dynamic signal sparsities of signals with N = 512 and M = 100. The vertical dashed lines represent Edyn at (a) 41.5698, (b) 3.8354, and (c) 33.3732.
6
Conclusion
We have presented a new dynamic CS framework that handles signals without the fixed sparsity assumption. The success probability of signal recovery in random sampling was investigated when the signal sparsity can vary. The success probability analysis in the existing CS framework was shown to be incapable of reflecting actual success probability by both theoretical analysis and experiments. On the contrary, our beta distribution modeling could closely reflect the actual success probability. We also considered signals that cannot be exactly represented with a sparse representation, where we could alternatively view the number of components included in the signal recovery as a varying quantity. This quantity was shown by both theoretical analysis and experiments to follow a right-tailed distribution such as the gamma distribution. We modeled the dynamic signal sparsity for these signals as well and provided an error analysis whose usefulness was demonstrated by experiments.
References [1] W. Bajwa, J. Haupt, A. Sayeed, and R. Nowak, “Compressive wireless sensing,” in Proc. IPSN, 2006, pp. 134–142. [2] S. Ji and L. Carin, “Bayesian compressive sensing and projection optimization,” in Proc. ICML, 2007, pp. 377–384. [3] M. W. Seeger and H. Nickisch, “Compressed sensing and bayesian experimental design,” in Proc. ICML, 2008, pp. 912–919. [4] C. Luo, F. Wu, J. Sun, and C. W. Chen, “Compressive data gathering for large-scale wireless sensor networks,” in Proc. MobiCom, 2009, pp. 145–156. [5] D. Hsu, S. Kakade, J. Langford, and T. Zhang, “Multi-label prediction via compressed sensing.” in Proc. NIPS, 2009, pp. 772–780.
10
[6] M. Lopes, “Estimating unknown sparsity in compressed sensing,” in Proc. ICML, 2013, pp. 217–225. [7] D. Malioutov and K. Varshney, “Exact rule learning via boolean compressed sensing,” in Proc. ICML, 2013, pp. 765–773. [8] R. G. Baraniuk, “Compressive sensing [lecture notes],” IEEE Signal Process. Mag., vol. 24, no. 4, pp. 118–121, Jul. 2007. [9] S. Foucart and H. Rauhut, A Mathematical Introduction to Compressive Sensing. 2013.
Springer,
[10] D. Lee and J. Choi, “Low complexity sensing for big spatio-temporal data,” in Proc. BigData, 2014, pp. 323–328. [11] D. Sejdinovic, C. Andrieu, and R. Piechocki, “Bayesian sequential compressed sensing in sparse dynamical systems,” in Proc. Allerton, 2010, pp. 1730–1736. [12] B. Shahrasbi, A. Talari, and N. Rahnavard, “TC-CSBP: Compressive sensing for timecorrelated data based on belief propagation,” in Proc. CISS, 2011, pp. 1–6. [13] J. Ziniel and P. Schniter, “Dynamic compressive sensing of time-varying signals via approximate message passing,” IEEE Trans. Signal Process., vol. 61, no. 21, pp. 5270–5284, Nov. 2013. [14] S. Ganguli and H. Sompolinsky, “Short-term memory in neuronal networks through dynamical compressed sensing,” in Proc. NIPS, 2010, pp. 667–675. [15] D. M. Malioutov, S. R. Sanghavi, and A. S. Willsky, “Sequential compressed sensing,” IEEE J. Sel. Top. Signal Process., vol. 4, no. 2, pp. 435–444, Apr. 2010. [16] E. J. Cand`es and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag., vol. 25, no. 2, pp. 21–30, Mar. 2008. [17] S. Boyd and L. Vandenberghe, Convex Optimization.
Cambridge University Press, 2004.
[18] Y. C. Pati, R. Rezaiifar, and P. S. Krishnaprasad, “Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition,” in Proc. ACSSC, 1993, pp. 40–44. [19] T. Blumensath and M. E. Davies, “Iterative thresholding for sparse approximations,” J. Fourier Anal. Appl., vol. 14, no. 5-6, pp. 629–654, Dec. 2008. [20] D. Lee and J. Choi, “Learning compressive sensing models for big spatio-temporal data,” in Proc. SDM, 2015, pp. 667–675. [21] R. G. Baraniuk, V. Cevher, M. F. Duarte, and C. Hegde, “Model-based compressive sensing,” IEEE Trans. Inf. Theory, vol. 56, no. 4, pp. 1982–2001, Apr. 2010. [22] S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition by basis pursuit,” SIAM J. Sci. Comput., vol. 20, no. 1, pp. 33–61, Jan. 1998. [23] G. Quer, R. Masiero, G. Pillonetto, M. Rossi, and M. Zorzi, “Sensing, compression, and recovery for WSNs: Sparse signal modeling and monitoring framework,” IEEE Trans. Wireless Commun., vol. 11, no. 10, pp. 3447–3461, Oct. 2012. 11
Appendix A The model functions in (7) are plotted in Fig. 5 with α’s and β’s obtained by the maximumlikelihood estimation for comparison. 250
250
200
200
150
Parameter Model estimation by (7) Parameter Model estimation by (7)
100
50
0 0.5
150
100
50
0.55
0.6
0.65 0.7 M/(ln(N)K)
0.75
0.8
0 0.85
Figure 5: Two beta distribution parameters as functions of M , N , and K estimated by (7). Fig. 6 also presents (sample) means and variances of the success probability distributions shown in Fig. 1, along with means and variances calculated using the model functions in (7). We can see that the beta distribution modeling with parameters in (7) effectively follows the actual success probability distributions. −3
Mean
0.5
0.5
Variance
x 10 1
1
Actual mean Mean estimation Actual variance Variance estimation 0 0.5
0.55
0.6
0.65 0.7 M/(ln(N)K)
0.75
0.8
0 0.85
Figure 6: Mean and variance of success probability distributions in Fig. 1 estimated with the beta distribution modeling α(·) and β(·).
Appendix B: Proof of Theorem 1 The inaccuracy of the failure probability in the CS framework results from the slowly decaying lower bound of success probability, that is, 1 − exp(−M/(C · ln(N )K)). In fact, we can show this lower bound decays slower than a power-law decay by the following lemma.
12
Lemma 1 There exists K0 > 0 such that for all K > K0 , the lower bound of recovery success probability is greater than the value of a power-law-decay function. Proof We need to show the following inequality M 1 − exp − > K −α C · ln(N )K
(15)
holds if K > K0 for some K0 > 0, where α > 0. Adding, subtracting, and taking the power K on both sides yields M −α K (1 − K ) > exp − . (16) C · ln(N ) We now use the binomial approximation on the left-hand side: (1 − K −α )K ≥ 1 − K · K −α . Thus we instead prove the following inequality M 1 − K 1−α > exp − . (17) C · ln(N ) holds if K > K0 for some K0 > 0. If we assume α > 1, then adding, subtracting, and taking the power 1/(1 − α) on both sides of (17) yields 1/(1−α) M 1 − exp − < K. (18) C · ln(N ) Setting K0 = (1 − exp(−M/(C · ln(N ))))1/(1−α) , we can argue that for all K > K0 , the lower bound of recovery success probability is greater than the value of a power-law-decay function. Corollary 2 In the CS framework, there is always a chance of succeeding at signal recovery however large K is. Proof The power-law-decay function K −α in (15) slowly converges to zero as K → ∞: its value is noticeably greater than zero even with large K. As the lower bound of recovery success probability is greater than the value of the power-law-decay function for all K > K0 , we can say there is always a chance of recovery success however large K is. We can now show that our beta distribution modeling provides more accurate success probability. The claim of the CS framework in Corollary 2 is in fact implausible because it says we can even set K > M and there is still a chance of success. We cannot expect signal recovery with a number of measurements M less than K. On the contrary, our beta distribution modeling can yield IP(s? = s | M, N, K) = 0 with a bounded Kmax . In particular, we show that the mean of Beta(αK , βK ), αK /(αK + βK ), converges to zero with a bounded Kmax . If we state αK /(αK + βK ) = 0 using (7), we obtain M aα ln + bα = 0. (19) ln(N )K Representing (19) with respect to K, we have Kmax
M = exp ln(N ) 13
bα aα
< ∞.
(20)
Similarly, we show this mean converges to one (IP(s? = s | M, N, K) = 1) with Kmin which is not so close to zero, whereas the lower bound of the recovery success probability given by the CS framework converges to one only if K is very close to zero. If we state αK /(αK + βK ) = 1 using (7), we have aβ
M + bβ = 0. ln(N )K
(21)
Representing (21) with respect to K, we have Kmin = −(aβ M )/(bβ ln(N )). Plugging Kmin into 1 − exp(−M/(C · ln(N )K)) shows that bβ 1 − exp > 0, (22) C · aβ as C > 0, aβ < 0, and bβ > 0: that is, 1 − exp(−M/(C · ln(N )K)) → 1 if and only if K → 0. Since 0 < Kmin < Kmax < ∞, we can argue that the beta distribution modeling can provide tighter recovery success probability.
Appendix C: Proof of Proposition 1 Since ksi k1 = C`1 for all i, we can conceive the same sequence {sn } of elements (absolute values) in si for all i. Then we have NX −K K ksi − si k1 = sn . (23) n=1
P −K Without loss of generality, we consider the partial sum N n=1 sn in (23) to be an arithmetic series that can be represented by a quadratic function in terms of K. We also assume the inequality constraint in (10) is the equality constraint such that ks? − sk2 = C1 · ks − sK k1 . q If we take the (partial) inverse function of the quadratic function, we have K ∼ Kmax − ks? − sk2 − (min ks?ij − si k2 ). Assuming the distribution of ks? − sk2 is symmetric (zero skewness), this asymptotic relation says ks? − sk2 will be compressed as it becomes large, which in turn makes the pdf of K right tailed. P −K A similar claim can be made if we consider the partial sum N n=1 sn to be a geometric series, ? where K ∼ N − log(ks − sk2 ). In this case, the pdf of K is skewed to the right as well.
Appendix D First, random signed spikes were artificially generated in different magnitudes at random locations and densified to perform random sampling. In particular, we considered an arithmetic sequence of length 50 (2, 4, 6, . . ., 98, 100), whose elements were placed at random locations in each vector. These signals are dense enough to be used for experiments because signal recovery always fails when K > 30 in our case, as shown in Fig. 2a. Fig. 7 displays the histogram of K, the number of components included in each signal recovery, which was obtained using the method explained in Section 4. We can identify that Proposition 1 actually holds here, as this distribution is skewed to the right. Furthermore, we empirically found that the distribution follows the gamma distribution, which is also natural since the gamma distribution has positive skewness, i.e., right tailed. 14
Actual distribution Gamma fitting
0.25
Density
0.2
0.15
0.1
0.05
0
20
22
24
26 K
28
30
32
Figure 7: Distribution of K fitted with a gamma distribution Gamma(242.81, 0.09), using the maximum likelihood estimation. Histogram was obtained with 300 different signals and 300 different experiments for each signal, for N = 512 and M = 100. We also provide results with real-world sensor data sets to verify Proposition 1 and the gamma distribution fitting. Figs. 8, 9, and 10 display the histograms of K and their gamma distribution fittings, where we can see that K follows the gamma distribution as well. Actual distribution Gamma fitting
0.09 0.08 0.07
Density
0.06 0.05 0.04 0.03 0.02 0.01 0
5
10
15
K
20
25
30
Figure 8: Distribution of K fitted with a gamma distribution Gamma(8.43, 1.71) for humidity (%) data. Histogram was obtained with 1,000 experiments over the same signal for N = 512 and M = 100.
15
0.04 Actual distribution Gamma fitting
0.035 0.03
Density
0.025 0.02 0.015 0.01 0.005 0
0
10
20
30 K
40
50
60
Figure 9: Distribution of K fitted with a gamma distribution Gamma(2.55, 7.94) for luminosity (V) data. Histogram was obtained with 1,000 experiments over the same signal for N = 512 and M = 100.
Actual distribution Gamma fitting
0.06
Density
0.05 0.04 0.03 0.02 0.01 0
0
5
10
15
K
20
25
30
35
Figure 10: Distribution of K fitted with a gamma distribution Gamma(5.75, 2.73) for temperature (◦ C) data. Histogram was obtained with 1,000 experiments over the same signal for N = 512 and M = 100.
16
