# Statistical Compressive Sensing of Gaussian Mixture Models

A new framework of compressive sensing (CS), namely statistical compressive sensing (SCS), that aims at efficiently sampling a collection of signals that follow a statistical distribution and achieving accurate reconstruction on average, is introduced. For signals following a Gaussian distribution, with Gaussian or Bernoulli sensing matrices of O(k) measurements, considerably smaller than the O(k log(N/k)) required by conventional CS, where N is the signal dimension, and with an optimal decoder implemented with linear filtering, significantly faster than the pursuit decoders applied in conventional CS, the error of SCS is shown tightly upper bounded by a constant times the k-best term approximation error, with overwhelming probability. The failure probability is also significantly smaller than that of conventional CS. Stronger yet simpler results further show that for any sensing matrix, the error of Gaussian SCS is upper bounded by a constant times the k-best term approximation with probability one, and the bound constant can be efficiently calculated. For signals following Gaussian mixture models, SCS with a piecewise linear decoder is introduced and shown to produce for real images better results than conventional CS based on sparse models.

## Authors

• 8 publications
• 72 publications
• ### Statistical Compressed Sensing of Gaussian Mixture Models

A novel framework of compressed sensing, namely statistical compressed s...
01/30/2011 ∙ by Guoshen Yu, et al. ∙ 0

01/25/2012 ∙ by Julio M. Duarte-Carvajalino, et al. ∙ 0

• ### Compressive Sensing via Low-Rank Gaussian Mixture Models

We develop a new compressive sensing (CS) inversion algorithm by utilizi...
08/27/2015 ∙ by Xin Yuan, et al. ∙ 0

• ### Online Adaptive Statistical Compressed Sensing of Gaussian Mixture Models

A framework of online adaptive statistical compressed sensing is introdu...
12/26/2011 ∙ by Julio Duarte-Carvajalino, et al. ∙ 0

• ### How to find real-world applications for compressive sensing

The potential of compressive sensing (CS) has spurred great interest in ...
05/06/2013 ∙ by Leslie N. Smith, et al. ∙ 0

• ### Towards Efficient Compressive Data Collection in the Internet of Things

It is of paramount importance to achieve efficient data collection in th...
05/29/2021 ∙ by Peng Sun, et al. ∙ 0

• ### More chemical detection through less sampling: amplifying chemical signals in hyperspectral data cubes through compressive sensing

Compressive sensing (CS) is a method of sampling which permits some clas...
06/27/2019 ∙ by Henry Kvinge, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## I Introduction

Compressive sensing (CS) attempts to achieve accurate signal reconstruction while sampling signals at a low sampling rate, typically far smaller than the Nyquist/Shannon rate. Let be a signal of interest, a non-adaptive sensing matrix (encoder), consisting of measurements, a measured signal, and a decoder used to reconstruct from . CS develops encoder-decoder pairs such that a small reconstruction error can be achieved.

Reconstructing from is an ill-posed problem whose solution requires some prior information (model) on the signal. Instead of the frequency band-limit signal model assumed in classic Shannon sampling theory, conventional CS adopts a sparse signal model (manifold models have been considered as well [baraniuk2009random, chen2010compressive]), i.e., there exists a dictionary, typically an orthogonal basis , a linear combination of whose columns generates an accurate approximation of the signal, , with the coefficients , , having their amplitude decaying fast after being sorted. For signals following the sparse model, it has been shown that using some random sensing matrices such as Gaussian and Bernoulli matrices with measurements, and an minimization or a greedy matching pursuit decoder promoting sparsity, with high probability CS leads to accurate signal reconstruction. The obtained approximation error is tightly upper bounded by a constant times the -best term approximation error, the minimum error that one may achieve by keeping the largest coefficients in  [candes2006robust, cohen2009compressed, donoho2006compressed].

The present paper introduces a novel framework of CS, namely statistical compressive sensing (SCS). As opposed to conventional CS that deals with one signal at a time, SCS aims at efficiently sampling a collection of signals and having accurate reconstruction on average. Instead of restricting to sparse models, SCS works with general Bayesian models. Assuming that the signals

follow a distribution with probability density function (pdf)

, SCS designs encoder-decoder pairs so that the average error

 Ex∥x−Δ(Φx)∥X=∫∥x−Δ(Φx)∥Xp(x)dx, (1)

where is a norm, is small. As an important example, SCS with Gaussian models is here shown to have improved performance (bounds) relative to conventional CS, the signal construction calculated with an optimal decoder implemented via a fast linear filtering. Moreover, for Gaussian mixture models (GMM), SCS with a piecewise linear decoder is introduced and shown to be very effective.

The motivation of SCS with Gaussian models is twofold. First, controlling the average error over a collection of signals is useful in signal acquisition, not only because one is often interested in acquiring multiple signals in real applications, but also because more effective processing of an individual signal, an image or a sound for example, is usually achieved by dividing the signal in local subparts, patches or short-time windows for instance, and a signal can be regarded as a collection of subpart signals [yu2008audio, yu2010PLE]. In addition, Gaussian mixture models (GMM), which model signals or subpart signals with a collection of Gaussians, have been shown effective in describing real signals, leading to excellent results in image inverse problems [yu2010PLE]

and missing data estimation

[leger2010Matrix].

After reviewing the optimal decoder for Gaussian signals in Section II, a quick numerical check of the good performance of Gaussian SCS is first given in Section III. Section IV analyzes this performance following a similar mathematical approach as the one adopted in conventional CS performance analysis [cohen2009compressed]. This result shows that with the same random matrices as in conventional CS, but with a considerably reduced number of measurements, and with the optimal decoder implemented with linear filtering, significantly faster than the decoders applied in conventional CS, the average error of Gaussian SCS is tightly upper bounded by a constant times the -best term approximation error with overwhelming probability, the failure probability being orders of magnitude smaller than that of conventional CS. Section V further shows stronger yet simpler results: For any sensing matrix, the average error of Gaussian SCS is upper bounded by a constant times the -best term approximation with probability one, and the bound constant can be efficiently calculated. Section VI

introduces, for SCS with GMM, a piecewise decoding scheme based on an efficient maximum a posteriori expectation-maximization (MAP-EM) algorithm following

[yu2010PLE], and shows better results on real images than those obtained with conventional CS.

In the rest of this paper, we will assume without loss of generality Gaussian signal models , with mean zero and diagonal covariance matrix , with

the sorted eigenvalues. We can always center the data with respect to the Gaussian distribution and make a basis change with principal component analysis (PCA)

[yu2010PLE]. For Gaussian and Bernoulli matrices that are known to be universal, analyzing CS in the canonical basis or the PCA basis is equivalent [baraniuk2008simple]. The Gaussian distributions are assumed to be full rank, i.e., , since a degenerated Gaussian can be regarded as a full-rank Gaussian within a reduced dimension.

## Ii Optimal Decoder for Gaussian SCS

It is well-known that the optimal decoders for Gaussian signals are calculated with linear filtering:

###### Theorem 1.

[kay1998fundamentals] Let

be a random vector with prior pdf

, and , , be a sensing matrix. From the measured signal , the optimal decoder that minimizes the mean square error (MSE) as well as the mean absolute error (MAE) where is any mapping from , is obtained with a linear MAP estimator,

 Δ(Φx)=argmaxxp(x|y)=SΦT(ΦSΦT)−1(Φx), (2)

and the resulting error is Gaussian with mean zero and with covariance matrix whose trace yields the MSE of SCS

 Ex[∥x−Δ(Φx)∥22]=Tr(S−SΦT(ΦSΦT)−1ΦS). (3)

In contrast to conventional CS, for which the minimization or greedy matching pursuit decoders, calculated with iterative procedures, have been shown optimal under certain conditions on and the signal sparsity [candes2006robust, donoho2006compressed], Gaussian SCS enjoys the advantage of having an optimal decoder (2) calculated via a linear filtering for any .

###### Corollary 1.

If a random matrix

is drawn independently to sense each , with all the other conditions as in Theorem 1, the MSE of SCS is

 Ex,Φ∥x−Δ(Φx)∥22=EΦ[Tr(S−SΦT(ΦSΦT)−1ΦS)]. (4)

The MSE of Gaussian SCS has closed-forms (3), (4), however, a mathematical analysis of (3) and (4) seems uneasy due to the complexity of the involved matrix expression. The next section evaluates these values through Monte Carlo simulations, preceding the theoretical bounds later developed.

## Iii Performance of Gaussian SCS – a numerical analysis at first

This section numerically evaluates the MSE of Gaussian SCS, and compares it with the minimal MSE generated by the best -term approximation. Under the Gaussian signal model with sorted eigenvalues , it is well-known that the best -term approximation error,

 σk({x})X=Ex∥x−xK∥X  and  σk({x})22=Ex∥x−xK∥22, (5)

is obtained with a linear projection to the first entries, i.e., , , and otherwise [mallat2008wts]. Note that . This best

-term approximation sensing is impractical with GMM, since the Gaussian assignment of a signal is unknown at the acquisition moment, GMM describing real data much better than a single Gaussian model

[yu2010PLE].

A power decay of the eigenvalues [mallat2008wts],

 λm=m−α,   1≤m≤N, (6)

, is assumed in the Monte Carlo simulation. An independent random Gaussian matrix realization is applied to sense each signal  [cohen2009compressed], and (4) is used to calculate the MSE of SCS.

Figures 1 (a) and (c)-top plot the MSE of SCS and that of the best -term approximation, as well as their ratio as a function of for Gaussian signals with a typical eigenvalue decay . With increasing, both MSEs decrease, their ratio being almost constant at about . The same is plotted in figures 1 (b) and (c)-bottom, with fixed at a typical value of , and varying from to . When increases, the eigenvalues decay faster, and the MSEs for both methods decrease. The ratio increases almost linearly with .

These results indicate a good performance of Gaussian SCS, its MSE is only a small number of times larger than that of the best -term approximation. The next sections provide mathematical analysis of this performance.

## Iv Performance Bounds of Gaussian SCS

Following [cohen2009compressed], this section shows that with Gaussian and Bernoulli random matrices of measurements, considerably smaller than the required by conventional CS, the average error of Gaussian SCS is tightly upper bounded by a constant times the -best term approximation with overwhelming probability, the failure probability being orders of magnitude smaller than that of conventional CS. The proofs of the theorems will not be given due to the space limitations. We consider only the encoder-decoder pairs that preserve , i.e., , satisfied by the optimal in (2) for Gaussian signals , .

### Iv-a From Null Space Property to Instance Optimality

The instance optimality in expectation bounds the average error of SCS with a constant times that of the best -term approximation (5), defining the desired SCS performance:

###### Definition 1.

We say that is instance optimal in expectation of order in , with a constant , if

 Ex,(Φ)∥x−Δ(Φx)∥X≤C0σk({x})X, (7)

the expectation with respect to , and to if one random is drawn independently for each . Similarly, the MSE instance optimality of order is defined as

 Ex,(Φ)∥x−Δ(Φx)∥22≤C0σk({x})22. (8)

The null space property in expectation defined next will be shown equivalent to the instance optimality in expectation.

###### Definition 2.

We say that in has the null space property in expectation of order in , with constant , if

 Ex,(Φ)∥η∥X≤Cσk({η})X,   ∀η=x−Δ(Φx), (9)

the expectation with respect to , and to if one random is drawn independently for each . Note that . Similarly, the MSE null space property of order is defined as

 Ex,(Φ)∥η∥22≤Cσk({η})22,   ∀η=x−Δ(Φx). (10)
###### Theorem 2.

Given an matrix , a norm , and a positive integer , a sufficient condition that there exists a decoder such that the instance optimality in expectation (7) holds with constant , is that the null space property in expectation (9) holds with for this . A necessary condition is the null space property in expectation (9) with . Similar results hold between MSE instance optimality (8) and null space property (10), with the constant in the sufficient condition.

Comparing to conventional CS that requires the null space property to hold with the best -term approximation error [cohen2009compressed], the requirement for Gaussian SCS is relaxed to , thanks to the linearity of the best -term approximation for Gaussian signals.

Theorem 2 proves the existence of the decoder for which the instance optimality in expectation holds for , given the null space property in expectation. However, it does not explain how such decoder is implemented. The following Corollary, a direct consequence of theorems 1 and 2, shows that for Gaussian signals the optimal decoder (2) leads to the instance optimality in expectation.

###### Corollary 2.

For Gaussian signals , if an sensing matrix satisfies the null space property in expectation (9) of order in , with constant , or the MSE null space property (10) of order with constant , then the optimal and linear decoder satisfies the instance optimality in expectation (7) in , or the MSE instance optimality (8).

### Iv-B From RIP to Null Space Property

The Restricted Isometry Property (RIP) of a matrix measures its ability to preserve distances, and is related to the null space property in conventional CS [candes2006near, donoho2006compressed]. The new linear RIP of order restricts the requirement of conventional RIP of order , to a union of -dimensional linear subspaces with consecutive supports:

###### Definition 3.

Let be a positive integer. Let define a linear subspace of functions with support in the first indices in , a linear subspace of functions with support in the next indices, and so on. The functions in the last linear subspace defined this way may have support with less than indices. An matrix is said to have linear RIP of order with constant if

 (1−δ)∥x∥2≤∥Φx∥2≤(1+δ)∥x∥2,   ∀ x∈∪Jj=1Kj. (11)

The linear RIP is a special case of the block RIP [eldar2009robust], with block sparsity one and blocks having consecutive support of the same size.

The following theorem relates the linear RIP (11) of a matrix to its null space property in expectation (9).

###### Theorem 3.

Let be a random vector that follows a certain distribution. Let be an matrix that satisfies the linear RIP of order with , and let be a decoder. Let . Assume further that decays in : , . Then satisfies the null property in expectation of order in  (9), with constant .

For Gaussian signals , with Gaussian or Bernoulli matrices, one realization drawn independently for each , and with the optimal decoder (2), the decay of assumed in in Theorem 3 is verified through Monte Carlo simulations.

### Iv-C From Random Matrices to RIP

The next Theorem shows that Gaussian and Bernoulli matrices satisfy the RIP for one subspace with overwhelming probability.

###### Theorem 4.

[baraniuk2008simple] Let be a random matrix of size drawn according to any distribution that satisfies the concentration inequality

 Pr(∥Φx∥22−∥x∥22|≥ϵ∥x∥22)≤2e−Mc0(δ/2),   ∀ x∈RN, (12)

where , and is a constant depending only on . Then for any set with , we have

 (1−δ)∥x∥2≤∥Φx∥2≤(1+δ)∥x∥2,   ∀ x∈XK, (13)

where is the set of all vectors in that are zero outside of , with probability greater than or equal to Gaussian and Bernoulli matrices satisfy the concentration inequality (12).

The linear RIP of order  (11) requires that (13) holds for subspaces. The next Theorem follows from Theorem 4 by simply multiplying by the probability that the RIP fails to hold for one subspace.

###### Theorem 5.

Suppose that , and are given. Let be a random matrix of size drawn according to any distribution that satisfies the concentration inequality (12). Then there exist constants depending only on such that the linear RIP of order  (11) holds with probability greater than or equal to for with the prescribed and .

Comparing with conventional CS, where the null space property requires that the RIP (13) holds for subspaces [baraniuk2008simple, candes2006near, donoho2006compressed], the number of subspaces in the linear RIP (11) is sharply reduced to for Gaussian SCS. In consequence, with the same number of measurements , the probability that a Gaussian or Bernoulli matrix satisfies the linear RIP is substantially higher than that for the conventional RIP. Equivalently, given the same probability that satisfies the linear RIP or the conventional RIP of order , the required number of measurements for the linear RIP is , substantially smaller than the required for the conventional RIP. Similar improvements have been obtained with model-based CS that assumes structured sparsity on the signals [baraniuk2010model].

With the results above, we have showed that for Gaussian signals, with sensing matrices satisfying the linear RIP (11) of order , for example Gaussian or Bernoulli matrices with rows, with overwhelming probability, and with the optimal decoder (2), SCS leads to the instance optimality in expectation of order in  (7), with constant . By the definition of CS, is typically small.

## V Performance Bounds of Gaussian SCS with RIP in Expectation

This section shows that with an RIP in expectation, a matrix isometry property more adapted to SCS, the Gaussian SCS MSE instance optimality (8) of order and constant , holds in the norm with probability one for any matrix. has a closed-form and can be easily computed numerically.

###### Definition 4.

Let be a random vector that follows a certain distribution. Let be an sensing matrix and let be a decoder. Let . in is said to have RIP in expectation in with constant if

 Ex,(Φ)∥ΦηK∥22=cKEx,(Φ)∥ηK∥22,   ∀ η=x−Δ(Φx), (14)

where , is the signal restricted to (, and otherwise), and the expectation is with respect to , and to if one random is drawn independently for each .

The conventional RIP is known to be satisfied only by some random matrices, Gaussian and Bernoulli matrices for example, with high probability. For a given matrix, checking the RIP property is however NP-hard [baraniuk2008simple]. By contrast, the constant of the RIP in expectation (14) can be measured for any matrix via a fast simulation, the quick convergence guaranteed by the concentration of measure [talagrand1996new]. The next proposition, directly following from (3) and (4), further shows that for Gaussian signals, the RIP in expectation has its constant in a closed form.

###### Proposition 1.

Assume , is an sensing matrix and is the optimal decoder (2). Then in satisfies the RIP in expectation in ,

 =cK(EΦ)[Tr(RKSRTK−RKSΦT(ΦSΦT)−1ΦSRTK)], (15)

where is an extraction matrix giving , i.e., , , all the other entries being zero. The expectation with respect to is calculated if one random is drawn independently for each .

The next Theorem shows that the RIP in expectation leads to the MSE null space property holding in equality.

###### Theorem 6.

Let be a random vector that follows a certain distribution, an sensing matrix, and a decoder. Assume and , for some . Assume that in has the RIP in expectation in with constant , and in with constant :

 Ex,(Φ)∥ΦηK∥22Ex,(Φ)∥ηK∥22=aK,   Ex,(Φ)∥ΦηKC∥22Ex,(Φ)∥ηKC∥22=bK,   ∀ η=x−ΔΦx. (16)

Then satisfies

 Ex,(Φ)∥η∥22=C0Ex,(Φ)∥ηKC∥22, (17)

where . In particular, if , then satisfies the MSE null space property of order , which holds with equality,

 Ex,(Φ)∥η∥22=C0σk({η})22. (18)

Let us check the MSE null space property constant of different sensing matrices in SCS for Gaussian signals , assuming the eigenvalues of follow a power decay (6) with , and . Gaussian, Bernoulli and random subsampling matrices of size are considered, and the optimal and linear decoder  (2) is applied to reconstruct the signals. For each matrix distribution, a different random matrix realization is applied to sense each signal . Note that since the random subsampling matrix , each row containing one entry with value 1 at a random position and 0 otherwise, has the maximal coherence with the canonical basis, this matrix is not suitable for directly sensing  [candes2007sparsity], and is replaced by in the simulation, with a DCT basis having low coherence with .

Figure 2 (a) plots , obtained by simulating (16), with , for different values of . When the number of SCS measurements increases, the reconstruction error of SCS decreases, resulting in a smaller ratio over the best- term approximation error with a fixed . Gaussian and Bernoulli matrices lead to similar values, slightly smaller than that of random subsampling matrices. Figure 2 (b) plots , as a function of , with . Gaussian and Bernoulli matrices lead to similar that varies little with . For random subsampling matrices slowly increases, almost linearly, and is equal to for a typical value , about 20% larger than that of Gaussian and Bernoulli matrices.

From Corollary 2 and Theorem 6, we obtain the next concluding result:

###### Theorem 7.

Assume . Let be an sensing matrix and the optimal and linear decoder (2). Then satisfies the MSE instance optimality of order  (8) with constant , and given in (16).

## Vi SCS with GMM – Experiments on Real Image Data

This section applies SCS with GMM in real image sensing and compares it with conventional CS based on sparse models.

Following [yu2010PLE], an image is decomposed into patches assumed to follow a GMM: There exist Gaussian distributions , and each image patch is independently drawn from one of these Gaussians with an unknown index . SCS samples each patch , with a possibly different for each . The decoding scheme (i) estimates , and (ii) identifies the Gaussian distribution that generates the -th patch, and reconstructs with the optimal and linear decoder (2) associated to the appropriate Gaussian, . The resulting piecewise linear decoder is implemented with a computationally efficient MAP-EM algorithm alternating between (i) and (ii) [yu2010PLE]. The algorithm has fast convergence and the MAP probability of the measured data increases as the algorithm iterates [yu2010PLE] (refer to [yu2010PLE] Sec. 2 for more details).

The dictionary for conventional CS is learned with K-SVD [aharon2006k] from 720,000 image patches, extracted from the entire standard Berkeley segmentation database containing 300 natural images [MartinFTM01]. The decoder is calculated with an minimization.

Figure 3 (a) shows the sensing performance on about 260,000 (sliding) patches, regarded as signals, extracted from the standard image Lena, as shown in Figure 4. The PSNRs generated by SCS and CS using Gaussian and random subsampling sensing matrices, one independent realization for each patch, are plotted as a function of the sampling rate . At the same sampling rate, SCS outperforms SC. The gain increases from about 0.5 dB at very low sampling rates (), learning GMM from the poor-quality measured data being more challenging, to more than 3.5 dB at high sampling rates (). (SC using an “oracle” dictionary learned from the full Lena itself, undoable in practice, improves its performance from 0.2 dB at low sampling rates to 1.3 dB at high sample rates, still lower than SCS.) For both SCS and CS, Gaussian and random subsampling matrices lead to similar PSNRs at low sampling rates (), and at higher sampling rates the former gains about 0.5 dB. Recall that SCS is not just more accurate and significantly faster, but also only uses the compressed image, while CS uses a pre-learned dictionary from a large database.

Figure 3 (b) shows the sensing performance on the whole Lena image. The sensing is performed on non-overlapped patches. Random subsampling matrices, which are diagonal operators (one non-zero entry per row), have the advantage of being able to lead to reconstruction on overlapped patches, averaging the overlapped reconstructed patches further improving the image estimation. The PSNRs generated by SCS using random subsampling matrices and overlapped reconstruction are plotted, in comparison with those obtained using Gaussian sensing matrices and non-overlapped reconstruction. The former improves from about 3.5 to 1.5 dB, at a cost of times computation. As illustrated in Figure 4, the overlapped reconstruction removes the block artifacts and considerably improves the reconstructed image.

Acknowledgements: Work partially supported by NSF, ONR, NGA, ARO, and NSSEFF. The authors thank very much Stéphane Mallat for co-developing the GMM framework for solving inverse problems.