Entropy-Based Learning of Sensing Matrices

TL;DR

This paper introduces an entropy-maximization learning approach to design orthogonal sensing matrices for compressed sensing, improving signal recovery accuracy and efficiency across various signal types.

Contribution

It presents a novel learning scheme that constructs sensing matrices by maximizing measurement entropy, with bounds for unique recovery and superior performance compared to existing methods.

Findings

Enhanced recovery accuracy for synthetic, speech, and image signals.

Reduced number of measurements needed for accurate reconstruction.

Outperforms existing sensing matrix construction methods.

Abstract

This paper proposes a learning method to construct an efficient sensing (measurement) matrix, having orthogonal rows, for compressed sensing of a class of signals. The learning scheme identifies the sensing matrix by maximizing the entropy of measurement vectors. The bounds on the entropy of the measurement vector necessary for the unique recovery of a signal are also proposed. A comparison of the performance of the designed sensing matrix and the sensing matrices constructed using other existing methods is also presented. The simulation results on the recovery of synthetic, speech, and image signals, compressively sensed using the sensing matrix identified, shows an improvement in the accuracy of recovery. The reconstruction quality is better, using less number of measurements, than those measured using sensing matrices identified by other methods.

Tables

Table 1. Table 1: RIP constants for varying values of M 𝑀 M for synthetic signals with sparsity 10.

No. of measurements (M)	Range of $\delta$
9	0.24-0.50
15	0.17-0.47
20	0.14-0.40
25	0.08-0.29
30	0.06-0.26

Equations

$$\mathbf{y}=\mathbf{\Phi }\mathbf{x}=\mathbf{\Phi \Psi }\mathbf{c}=\mathbf{Ac},$$

$$(1-\delta )\Vert \mathbf{c}{\Vert }_2^2\le \Vert \mathbf{Ac}{\Vert }_2^2\le (1+\delta )\Vert \mathbf{c}{\Vert }_2^2,\text{where }0<\delta <1.$$

$$H(\mathbf{x}\mid \mathbf{\Psi })=-\sum _{i=1}^N{p}_i\ln (p_i)=-\sum _{i=1}^N\frac{|c_i{|}^2}{\Vert \mathbf{c}{\Vert }_2^2}\ln \left(\frac{|c_i{|}^2}{\Vert \mathbf{c}{\Vert }_2^2}\right).$$

$$H(\mathbf{y})=-\sum _{i=1}^M\frac{|y_i{|}^2}{\Vert \mathbf{y}{\Vert }_2^2}\ln \left(\frac{|y_i{|}^2}{\Vert \mathbf{y}{\Vert }_2^2}\right),$$

$$n_{th,c}^{\Psi }=\lceil \exp (H(\mathbf{c}))\rceil .$$

$$H(\mathbf{c})\le \ln (K)<H(\mathbf{y})\le \ln (M).$$

$$\exp (H(\mathbf{y}))\le K\Rightarrow M_{eff}\le K+ϵ.$$

$$M_{eff}=2K-\nu .$$

$$M-M_{eff}\to 0\Rightarrow \nu \to 0.$$

$${\hat{\mathbf{Y}}}_j=\text{arg}\underset{{\mathbf{Y}}_j}{\max }H({\mathbf{Y}}_j),$$

$${\hat{\mathbf{Y}}}_j=\text{arg}\underset{{\mathbf{Y}}_j}{\max }\sum _{i=1}^M-\frac{|y_{ij}{|}^2}{\Vert {\mathbf{Y}}_j{\Vert }_2^2}\ln \left(\frac{|y_{ij}{|}^2}{\Vert {\mathbf{Y}}_j{\Vert }_2^2}\right),$$

$$(1-\delta )\le \frac{\Vert {\mathbf{AC}}_j{\Vert }_2^2}{\Vert {\mathbf{C}}_j{\Vert }_2^2}\le (1+\delta ).\text{vspace}-2.5pt$$

$${\hat{\mathbf{Y}}}_j=\text{arg}\underset{{\mathbf{Y}}_j}{\max }\sum _{i=1}^M-\frac{|y_{ij}{|}^2}{\Vert {\mathbf{Y}}_j{\Vert }_2^2}\ln \left(\frac{|y_{ij}{|}^2}{\Vert {\mathbf{Y}}_j{\Vert }_2^2}\right)$$

$$\text{subject to }{\left\{{\left(\frac{\Vert {\mathbf{Y}}_j{\Vert }_2}{\Vert {\mathbf{C}}_j{\Vert }_2}\right)}^2-1\right\}}^2\le {\delta }^2.$$

$${\hat{Y}}_j=\text{arg}\underset{{\mathbf{Y}}_j}{\min }\left\{\sum _{i=1}^M\frac{|y_{ij}{|}^2}{\Vert {\mathbf{Y}}_j{\Vert }_2^2}\ln \left(\frac{|y_{ij}{|}^2}{\Vert {\mathbf{Y}}_j{\Vert }_2^2}\right)+\alpha \left|\left({\left\{{\left(\Vert {\mathbf{Y}}_j{\Vert }_2/\Vert {\mathbf{C}}_j{\Vert }_2\right)}^2-1\right\}}^2\right)-{\delta }^2\right|\right\}.$$

$$\hat{\mathbf{A}}=\text{arg}\underset{\mathbf{A}}{\min }\Vert \hat{\mathbf{Y}}-\mathbf{AC}{\Vert }_F^2.$$

$$\hat{\mathbf{R}}=\text{arg}\underset{\mathbf{R}}{\min }\Vert \mathbf{RP}-\mathbf{Q}{\Vert }_F^2\text{s.t. }{\mathbf{RR}}^T=I,$$

$$\hat{\mathbf{\Phi }}=\hat{\mathbf{A}}{\mathbf{\Psi }}^{-1},$$

$${\hat{Y}}_j=\text{arg}\underset{{\mathbf{Y}}_j}{\min }\sum _{i=1}^M\frac{|y_{ij}{|}^2}{\Vert {\mathbf{Y}}_j{\Vert }_2^2}\ln \left(\frac{|y_{ij}{|}^2}{\Vert {\mathbf{Y}}_j{\Vert }_2^2}\right)+$$

$$\alpha \left|\left({\left\{{\left(\Vert {\mathbf{Y}}_j{\Vert }_2/\Vert {\mathbf{C}}_j{\Vert }_2\right)}^2-1\right\}}^2\right)-{\delta }^2\right|.$$

$$SRER=10{\log }_{10}\frac{{\sum }_i{x}_i^2}{{\sum }_i(x_i-\widehat{x_i}{)}^2},\text{vspace}-2.5pt$$

$$\mathbf{y}$$

$$=\mathbf{\Phi \Psi }(\mathbf{c}+{\mathbf{c}}_{\mathbf{e}})=\mathbf{A}(\mathbf{c}+{\mathbf{c}}_{\mathbf{e}}).$$

$$\text{rank}(\tilde{\mathbf{A}})=M_{eff}$$

$$\text{spark}(\tilde{\mathbf{A}})\in [2,M_{eff}+1].$$

$$2K<\text{spark}(\tilde{\mathbf{A}})\le M_{eff}+1$$

$$\Rightarrow 2K<M_{eff}+1$$

$$\Rightarrow 2K\le M_{eff}.$$

$$\Vert \hat{\mathbf{Y}}-\mathbf{AC}{\Vert }_F^2$$

$$=tr(\hat{\mathbf{Y}}{\hat{\mathbf{Y}}}^T)-2tr(\hat{\mathbf{Y}}{\mathbf{C}}^T{\mathbf{A}}^T)+tr({\mathbf{ACC}}^T{\mathbf{A}}^T)$$

$$tr(\hat{\mathbf{Y}}{\mathbf{C}}^T{\mathbf{A}}^T)$$

$$=tr({\mathbf{AU}\mathbf{\Delta }\mathbf{V}}^T)$$

$$tr({\mathbf{AU}\mathbf{\Delta }\mathbf{V}}^T)$$

$$=tr({\mathbf{V}}^T{\mathbf{A}{\mathbf{U}}_{\mathbf{M}}{\mathbf{\Delta }}_{\mathbf{M}}}^T)$$

Full text

1 Introduction

Compressed sensing (CS) aims at capturing signals, sparse in some domain, in a reduced set of measurements. Consider a signal $\bm{x}\in\mathbb{R}^{N}$, having representation $\bm{c}\in\mathbb{R}^{N}$ relative to a basis $\bm{\Psi}=\{\bm{\psi_{i}}\}_{i=1}^{N}$, such that $\bm{x}=\bm{\Psi c}$. The representation of the signal is $K$- sparse if $\|\bm{c}\|_{0}\leq K$ (which is the number of non-zero coefficients in the representation) for $K\ll N$. For the signal $\bm{x}$, the compressive sensing problem is given as

[TABLE]

where the sensing matrix $\bm{\Phi}\in\mathbb{R}^{M\times N}$, $K<M\ll N$, and $\bm{A}=\bm{\Phi\Psi}$.

The $\bm{\Phi}$ matrix should be constructed such that the matrix $\bm{A}$ captures maximum information from the sparse coefficients $\bm{c}$ for a known sparsifying basis $\bm{\Psi}$. The classical choice of the sensing matrix is a random matrix which simplifies the theoretical analysis [1][2]. The scheme for realizing the random projections is the random demodulator proposed in [3] and [4]. In this paper, we address the problem of identifying an efficient sensing matrix.

In [5], Elad introduced a structured sensing matrix along with a method to construct it by reducing the mutual coherence $\mu(\bm{A})$ of the columns of the matrix $\bm{A}$. Subsequently, many techniques were proposed to construct the sensing matrix by reducing $\mu(\bm{A})$ [7] -[10]. In [11] and [12], the sensing matrices were constructed by applying multidimensional scaling (MDS) on a sparsifying dictionary $\bm{\Psi}$. The methods employed in error control coding theory have also been used for the construction of sensing matrices [13] -[15]. Carson et al [16] proposed the construction of a projection matrix by maximizing the statistical Renyi entropy of the projections. It requires the knowledge of the statistical probability distribution of the signals. Baldassarre et al [17] proposed a method for constructing a measurement matrix for signals with structured sparsity by learning from the set of signals. The technique proposed in [17] is a sub-sampling method which aims at capturing maximum energy of the structured sparse signals.

In [18] and [19], Hegde et al proposed a learning method called NuMax to construct a sensing matrix such as to satisfy the restricted isometry property (RIP)[6].

[TABLE]

In this paper, we propose a learning method to construct an efficient measurement matrix (Section 3), from a set of training signals belonging to a class of signals, by maximizing the Shannon entropy of the measurement vectors constrained to achieving RIP. The measurement matrix constructed using the proposed method has orthogonal rows. The probability distribution of the signals and the Shannon entropy used in this article (Section 2) are different from the statistical probability and entropy used in [16], respectively. A relation between the entropy $H(\bm{y})$ of the measurement vector $\bm{y}$, the number of measurements $M$, and the entropy $H(\bm{c})$ of the vector of representation coefficients of the signal (Section 2) is also proposed. The simulation results (Section 4) with synthetic and speech signals suggest that the method proposed is capable of constructing measurement matrices that give improved recovery from a reduced set of measurements. An advantage of the proposed method is that it works even for signals that do not have structured sparsity (Section 4.1.1).

2 Motivation and Problem Formulation

To motivate the entropy based measurement matrix design, we use the definitions of the probability distribution of the representation of a signal and the Shannon entropy of the representation of the signal as proposed in [20] - [27]:

Definition 1**.**

Let $\bm{\Psi}=\{\bm{\psi}_{i}\}_{i=1}^{N}$ be a basis of an $N$-dimensional space. Let $\bm{x}$ be a signal belonging to a class of signals $\bm{X}$ such that $\bm{x}=\sum_{i=1}^{N}{c_{i}\bm{\psi}_{i}}$, where $\bm{c}=[c_{1},c_{2},...c_{N}]^{T}$ is the vector of representation coefficients of $\bm{x}$ relative to $\bm{\Psi}$. The probability distribution of the representation of the signal $\bm{x}$ relative to $\bm{\Psi}$ is $\bm{p}=\{p_{i}\}_{i=1}^{N}$, where $p_{i}=|c_{i}|^{2}/\|\bm{c}\|_{2}^{2}$.

The entropy of representation conditioned to $\bm{\Psi}$ is the amount of information left in $\bm{x}$ when the representation basis $\bm{\Psi}$ is known [22]. This conditional entropy $H(\bm{x}\mid\bm{\Psi})$ is defined as follows.

Definition 2**.**

The Shannon entropy of the signal $\bm{x}$ with respect to basis $\bm{\Psi}$ is given by

[TABLE]

Remarks: By definition, $0\times\ln(0)=0$ [28]. The terms Shannon entropy and entropy are used interchangeably in this article.

Since the entropy of representation with respect to the basis $\bm{\Psi}$ depends only on the probability distribution of representation $\bm{p}$, which depends on the coefficients of representation $\bm{c}$, we can conclude $H(\bm{x}\mid\bm{\Psi})=H(\bm{c})$, with $H(\bm{c})$ being the entropy of the representation of $\bm{x}$ relative to $\bm{\Psi}$. The entropy of representation of the reduced set of measurements $\bm{y}$ with respect to some basis of the range space of $\bm{A}$ is termed as the entropy of the set of measurements or entropy of the measurement vector. For simplicity, we consider the representation basis of the measurements to be the standard ordered basis of an $M$-dimensional space. Hence the entropy of the measurement vector is

[TABLE]

where $\bm{y}$ is given by (1) and $y_{i}$’s are the representation coefficients relative to the standard ordered basis of the reduced space of measurements which is the range space of $\bm{A}$.

From the definition of the entropy of representation $H(\bm{c})$, it can be seen that the more concentrated the probability distribution of the representation, the lower is the entropy. Since the probability of representation is directly proportional to the magnitude of the representation coefficient, it can be argued that the lower the entropy $H(\bm{c})$, the higher is the compressibility of the representation [23]. To quantify the compressibility of the signal with respect to a basis, we introduce theoretical dimension as the entropy-based measure of sparsity. Given the entropy of representation $H(\bm{c})$, the theoretical dimension $n_{th,c}^{\Psi}$ of the representation of a signal in the basis $\bm{\Psi}$ is given by [21] [23] [24]-[27]

[TABLE]

The theoretical dimension takes values $1\leq n_{th,c}^{\Psi}\leq N$, where $N$ is the total number of basis vectors in the representation basis. As entropy decreases, the theoretical dimension also decreases leading to a compressible representation of the signal. The theoretical dimension specifies the number of basis vectors required to represent a compressible signal without unduly degrading the signal quality. Experimental results of [23] show that the theoretical dimension gives the number of basis vectors required to capture at least $90\%$ of the signal energy. The various criteria that a measure of sparsity must satisfy are discussed in [26]. Meena and Abhilash [27] discuss the various criteria, mentioned in [26], that the theoretical dimension based measure of sparsity satisfies.

For a compressible signal, the $l_{0}$-sparsity is achieved by restricting the representation to those coefficients which carry the maximum energy of the signal, such that the signal quality is not degraded unduly. Hence, the theoretical dimension of the signal representation $\bm{c}$ can be approximated as the $l_{0}$- sparsity $K$ of the representation of the signal, that is $K\approx n_{th,c}^{\Psi}$. According to the definition of theoretical dimension, the lower the entropy, the higher is the sparsity or the lower is the value of $K$. The theoretical dimension of $\bm{y}$ with respect to the standard ordered basis of an $M$ dimensional space is $n_{th,y}^{I}=\lceil\exp(H(\bm{y}))\rceil=M_{eff}$. The quantity $M_{eff}$ is the effective number of coefficients, in the representation of the measurement vector $\bm{y}$ relative to the standard ordered basis, that captures at least $90\%$ of the signal energy. The remaining $M-M_{eff}$ number of coefficients in the representation of $\bm{y}$ carries insignificant amount of information.

According to the theory of CS, for any measurement vector $\bm{y}\in R^{M}$ there exists at most one $K$-sparse signal $\bm{c}\in R^{N}$, such that $\bm{y}=\bm{Ac}$, if and only if $\text{spark}(\bm{A})>2K$ [6], where $\text{spark}(\bm{A})$ is the smallest number of columns of $\bm{A}$ that are linearly dependent. The value of spark($\bm{A}$) lies in the range $[2,M+1]$, where $M$ is the dimension of the range space of $\bm{A}$ [6]. Hence the requirement $M\geq 2K$ holds good. Since $M_{eff}$ gives the effective number of measurements, ideally, the necessary condition for unique recovery would be $M_{eff}\geq 2K$ (See Appendix A). Since we use the entropy of the measurement vector $\bm{y}$ to construct the measurement matrix, we identify the bounds on the entropy of the measurement vector, for a fixed $M$. These bounds are necessary, but not sufficient, for achieving unambiguous measurements for unique recovery. These bounds are presented in Lemma 1.

Lemma 1**.**

If $\bm{y}$ is a measurement vector consisting of a set of $M\geq 2K$ measurements of a non-zero compressible signal, then the entropy $H(\bm{y})$ of the measurement vector, relative to the standard ordered basis of an $M$-dimensional space, satisfies

[TABLE]

where $H(\bm{c})$ is the entropy of the representation of a compressible signal relative to a sparsifying basis and the approximate $l_{0}$ sparsity of the signal is $K\approx\lceil\exp(H(\bm{c}))\rceil$.

Proof.

If $H(\bm{y})\leq\ln(K)$, then $\exp(H(\bm{y}))\leq K$. We know that $M_{eff}=\lceil\exp(H(\bm{y}))\rceil$. Hence $M_{eff}=\exp(H(\bm{y}))+\epsilon$, with $0\leq\epsilon<1$. Thus,

[TABLE]

Therefore, the case $H(\bm{y})\leq\ln(K)$ violates the requirement $M_{eff}\geq 2K$. Hence by contradiction, $\ln(K)<H(\bm{y})$. Since by (5), $H(\bm{c})$ is at the most $\ln(K)$, the lowest bound in (6) is true.

Given the number of measurements $M$, the entropy $H(\bm{y})$ of the measurement $\bm{y}$ is at the most $\ln(M)$ (by (4)); hence the upper bound in (6) also holds. Consider two signals $\bm{x_{1}}$ and $\bm{x_{2}}$ $(\bm{x_{1}}\neq\bm{x_{2}})$, with measurements $\bm{y_{1}}$ and $\bm{y_{2}}$ such that $H(\bm{y_{1}})=H(\bm{y_{2}})=\ln(M)$, then by Definitions 1 and 2, the probability distribution of the $j$-th measurement vector $\bm{y_{j}}$ with respect to the standard ordered basis of a vector space of dimension $M$ is $p_{ij}=\frac{\mid y_{ij}\mid^{2}}{\|\bm{y_{j}}\|_{2}^{2}}=1/M$, for $i=1,2...M$. The equality of the probability distribution does not imply equality of the measurement vectors. That is, $\frac{\mid y_{i1}\mid^{2}}{\|\bm{y_{1}}\|_{2}^{2}}=\frac{\mid y_{i2}\mid^{2}}{\|\bm{y_{2}}\|_{2}^{2}}$ does not imply $y_{i1}=y_{i2}$. Hence the equality $H(\bm{y})=\ln(M)$ does not affect the uniqueness of the measurements. ∎

Let $\mathcal{N}(\bm{A})$ represent the null space of the matrix $\bm{A}$. If $\bm{c}\in\mathcal{N}(\bm{A})$ (where $\bm{c}$ is the $K$-sparse approximation of the representation of a non-zero compressible signal), then the vector of the reduced set of measurements is $\bm{y}=\bm{Ac}=\bm{\theta}$, where $\bm{\theta}\in\mathbb{R}^{M}$ is the zero vector. Since the probability distribution of the measurements is defined as $p_{i}=\left\{\frac{|y_{i}|^{2}}{\|\bm{y}\|_{2}^{2}}\right\}_{i=1}^{M}$, the entropy $H(\bm{y})$ of $\bm{y}$ is undefined for $\bm{c}\in\mathcal{N}(\bm{A})$, as $p_{i}$’s are undefined. Hence, $\bm{c}\in\mathcal{N}(\bm{A})$ is undesirable. The theory of CS also discusses that the measurement operator $\bm{A}$ can uniquely represent all non-zero $K$-sparse signals if and only if no non-zero $2K$-sparse signal lies in $\mathcal{N}(\bm{A})$ [6]. The following Lemma shows that maximization of the entropy $H(\bm{y})$ of measurements implies that no non-zero $K$-sparse signal lies in $\mathcal{N}(\bm{A})$. Ideally, if no non-zero $K$-sparse signal falls in $\mathcal{N}(\bm{A})$, then no non-zero $2K$-sparse signal will belong to $\mathcal{N}(\bm{A})$.

Lemma 2**.**

If $\bm{A}$ is a measurement matrix having rank $M$ and $\bm{y}\in\mathbb{R}^{M}$ is a vector of measurements of the $K$-sparse approximation of a non-zero compressible signal with $M\geq 2K$, then the maximization of the entropy $H(\bm{y})$ of $\bm{y}$ implies that no $K$-sparse signal falls in the null space $\mathcal{N}(\bm{A})$ of $\bm{A}$.

Proof.

Since the rank of $\bm{A}$ is $M$, $\text{spark}(\bm{A})\in[2,M+1]$. Since $M\geq 2K$, $\text{spark}(\bm{A})>2K$, implying the existence of at most one $K$-sparse signal $\bm{c}\in\mathbb{R}^{N}$ such that $\bm{y}=\bm{Ac}$.

Let $\bm{c}$ represent the vector of sparse representation coefficients of a non-zero compressible signal approximated to be $K$-sparse. If $\bm{c}\in\mathcal{N}(\bm{A})$, then the vector of the reduced set of measurements is $\bm{y}=\bm{Ac}=\bm{\theta}$. If $\|\bm{y}\|_{0}=1$, then $\bm{c}\not\in\mathcal{N}(\bm{A})$ and $H(\bm{y})=0$. But, by Lemma 1, for $K=1$, $0<H(\bm{y})\leq\ln(M)$. Hence, $H(\bm{y})=0$ indicates that the measurements are incomplete. Therefore, to capture maximum information of the signal into $M$ measurements, the entropy $H(\bm{y})$ of $\bm{y}$ should be maximized. An entropy maximized non-zero measurement $\bm{y}$ implies that $\bm{c}\not\in\mathcal{N}(\bm{A})$. ∎

Lemma 2 maintains that maximization of the entropy $H(\bm{y})$ guarantees the null space property of $\bm{A}$. However, it does not claim that the maximization of entropy $H(\bm{y})$ would make the set of measurements complete. But, for unique recovery the set of measurements should capture maximum information contained in the signal. In particular, if the number of measurements is as small as possible, then the compression gain is the highest. Lemma 3 establishes how maximization of the entropy $H(\bm{y})$ enables unique recovery with the least possible number of measurements.

Lemma 3**.**

The maximization of the entropy $H(\bm{y})$ of the set of measurements of a compressible signal leads to unique recovery with a reduced number of measurements $M$ close to $2K$.

Proof.

$M_{eff}$ is the number of non-zero coefficients of the measurement vector in the standard ordered basis, that capture at least $90\%$ of the energy. By definition $M_{eff}\leq M$. Let $M_{eff}=M-\nu$ with $0\leq\nu$. As mentioned earlier, the necessary condition for unique recovery is $M\geq 2K$ [6]. If $M=2K$ then,

[TABLE]

Since $M_{eff}$ gives the effective number of measurements, ideally the necessary condition for unique recovery would be $M_{eff}\geq 2K$. This is not satisfied by (8). Hence the ideal lower bound on $M$ would be $2K+\nu$.

Since $M_{eff}=\lceil\exp(H(\bm{y}))\rceil$, increasing $H(\bm{y})$ leads to increase in $M_{eff}$, thus

[TABLE]

From (8) and (9) we see that increasing $H(\bm{y})$ results in reducing $\nu$ and hence $M_{eff}\approx M$. Hence, as $\nu$ decreases, the lower bound $2K+\nu$ on $M$ comes close to $2K$. Thus, maximization of $H(\bm{y})$ leads to unique recovery with a reduced number of measurements $M$ close to $2K$. ∎

Based on Lemmas 1, 2 and 3, it can be concluded that the matrix $\bm{A}$ that maximizes $H(\bm{y})$ would capture maximum information from the vector of coefficients $\bm{c}$ into a reduced set of measurements $\bm{y}$. It also implies that the smallest required number of measurements $M$ could be as small as $2K$. Hence, we propose a learning scheme for identifying $\bm{A}$, and thus $\bm{\Phi}$ for a class of signals such that $\bm{\Phi}$ maximizes the entropy $H(\bm{y})$ of the measurement vector $\bm{y}$.

3 Entropy Maximizing Sensing (EMS) Matrix Design

We propose a two-stage learning procedure to identify an efficient measurement matrix $\bm{\Phi}$ for a class of compressible signals. The learning method is motivated by the two stage dictionary/transform learning algorithms [23][29] [30]. The dictionary learning also finds application in Blind Compressive Sensing (BCS) [31], where the measurement matrix $\bm{A}=\bm{\Phi\Psi}$ (with $\bm{\Psi}$ unknown) is learned using the dictionary learning approach. In this paper, we do not consider the BCS framework.

The first stage of the proposed algorithm finds a set of entropy maximized measurements of the signals in the training set $\bm{X}$. The second stage tries to learn a $\bm{\Phi}$ based on the desired measurements, sparsifying basis, and the set of training signals. The two stages of the algorithm are alternately performed for a fixed number of iterations.

3.1 Stage I

Consider the matrix of $N$ dimensional training signals $\bm{X}\in\mathbb{R}^{N\times L}$, where each column is a signal from the training set and $L$ is the number of signals in the training set. Let $\bm{\Psi}\in\mathbb{R}^{N\times N}$ be a known orthonormal sparsifying basis (with the basis vectors $\{\bm{\psi}_{i}\}_{i=1}^{N}$ arranged as its columns) such that $\bm{X}=\bm{\Psi C}$, where $\bm{C}\in\mathbb{R}^{N\times L}$ is the matrix of the representation coefficients of the signals in $\bm{X}$ with respect to $\bm{\Psi}$. That is, $\bm{C}_{p}$ (the p-th column of $\bm{C}$) is the representation of the p-th signal $\bm{X}_{p}$ with respect to $\bm{\Psi}$. Given $\bm{\Psi}$ and an initial measurement matrix $\bm{\Phi}\in\mathbb{R}^{M\times N}$, we need to find the reduced set of measurements $\bm{Y}\in\mathbb{R}^{M\times L}$, of the training signals in $\bm{X}$, where $\bm{Y}=\bm{\Phi X}=\bm{\Phi\Psi C}=\bm{AC}$, such that the entropy of the set of measurements of each signal in $\bm{X}$ is maximized. Since for the $p$-th signal $\bm{X_{p}}$, the measurement vector $\bm{Y_{p}}=\bm{AC_{p}}$ is independent of the measurement of the $q$-th signal $\bm{Y_{q}}=\bm{AC_{q}}$ $(p\neq q)$, we can update the measurements considering all the signals separately. Thus, the problem can be formulated as

[TABLE]

where $\bm{Y}_{j}$ is the $j$-th column of $\bm{Y}$, and $\widehat{\bm{Y}}_{j}$, the $j$-th column of $\widehat{\bm{Y}}$, is the measurement vector of the $j$-th signal $\bm{X}_{j}$ having maximized entropy. Using Definition 3, we can rewrite the problem as

[TABLE]

where $\bm{Y}_{j}=\{y_{ij}\}_{i=1}^{M}$. The RIP in (2) can be restated as

[TABLE]

Thus, we can reformulate the problem (11), to satisfy the RIP, as

[TABLE]

The constrained problem in (13) can be made unconstrained by using the penalty method. To incorporate the penalty, we convert the maximization problem to a minimization problem. Hence the stage I solution is given as

[TABLE]

Since the absolute value function of the penalty term is non-differentiable at zero, we approximate it using the relaxation $|\bm{z}|\approx|\bm{z}|^{\zeta}=\sqrt{\bm{z}^{*}\bm{z}+\zeta}$ where $\zeta$ is taken to be $10^{-15}$ [32]. The value of $\delta$ is chosen depending on the desired RIP constant. The bounds on the value of $\delta$ for signals compressible in a dictionary are discussed in [33]. We have arbitrarily chosen $\delta=0.1$ for the experiments in the paper. The problem (14) can be solved using any optimization algorithm with the initial vector as $\bm{Y}_{j}=\bm{\Phi\Psi C}_{j}$.

3.2 Stage II

Stage I gives the desired measurements $\widehat{\bm{Y}}$ that maximizes the entropy of the measurement vectors. In stage II, we identify the matrix $\bm{A}$ and hence $\bm{\Phi}$ that would lead to $\widehat{\bm{Y}}$. The problem can be formulated as

[TABLE]

where the matrix $\bm{A}$ has orthogonal rows and $\|.\|_{F}$ represents the Frobenius norm. To solve the problem in (15), we use the orthogonal Procrustes method. The orthogonal Procrustes problem [30] [34] [35] is:

[TABLE]

where $\widehat{\bm{R}}$ is an orthogonal square matrix to be found such that $\widehat{\bm{R}}$ acts on a matrix $\bm{P}$ to result in the matrix $\bm{Q}$. Considering the singular value decomposition of $\bm{PQ}^{T}$ as $\widetilde{\bm{U}}\widetilde{\bm{\Delta}}\widetilde{\bm{V}}^{T}$, the solution to (16) is $\bm{R}=\widetilde{\bm{V}}\widetilde{\bm{U}}^{T}$.

Since the matrix $\widehat{\bm{A}}$ in (15) is not a square matrix, the solution to the Orthogonal Procrustes Problem has to be altered. If the singular value decomposition of $\bm{C}\widehat{\bm{Y}}^{T}$ is given by $\bm{U\Delta V}^{T}$, where $\bm{U}\in\mathbb{R}^{N\times N}$ and $\bm{V}\in\mathbb{R}^{M\times M}$, we propose that the solution to the problem (15) is given by $\widehat{\bm{A}}=\bm{VU_{M}}^{T}$, where $\bm{U_{M}}$ contains the first $M$ columns of $\bm{U}$ which correspond to the largest $M$ singular values of $\bm{C\widehat{Y}}^{T}$ (see Appendix B). The desired measurement matrix $\widehat{\bm{\Phi}}$ can be obtained from $\widehat{\bm{A}}$ as

[TABLE]

which holds good because $\bm{\Psi}$ is a well conditioned matrix of the representation basis for the class of signals considered.

Fig. 1 (a) and Fig. 1(b) show the variation of the average entropy per signal as the iteration progresses for synthetic and speech signals, respectively; the variation saturates after a finite number of iterations indicating the convergence of the algorithm empirically. The figures were generated by setting $\alpha=1$ and $\delta=0.1$ in the algorithm.

The value of $\alpha$ in stage I decides the rate of convergence and the value of the maximum entropy to which the convergence occurs. Large values of $\alpha$ would mean that the RIP constraint is strictly followed in stage I and the maximum entropy $H(\widehat{\bm{Y}}_{j})$ obtained in stage I of each iteration would be less. Since the stage II finds the matrix $\bm{A}$ such that the measurements $\bm{AC}_{j}$ is as close as possible to $\widehat{\bm{Y}}_{j}$, the final entropy $H(\bm{\widehat{A}C}_{j})$ at each iteration would be in small increments. Hence, the algorithm would converge slowly to a lower value of maximum entropy. For smaller values of $\alpha$, the RIP constraint will not be followed strictly in stage I and the maximum entropy obtained in each iteration would be close to $\ln(M)$. In such a case, the algorithm will converge fast to a higher entropy value. The discussion is validated in the Fig. 2, which shows the saturation of average entropy per signal for synthetic signals for different values of $\alpha$. The experimental setup in this paper uses $\alpha=1$.

3.3 Discussion

Consider the representation $\bm{C}_{j}\in\bm{C}$ of the $j$-th signal $\bm{X_{j}}\in\bm{X}$ with respect to the basis $\bm{\Psi}$ which is measured by $\bm{A}$ to get the measurement vector $\bm{Y}_{j}$. Let $\bm{\Lambda}$ be the set of indices corresponding to the coefficients in $\bm{C}_{j}$ having high magnitudes, which capture at least $90\%$ energy of the signal. Let $\bm{\Lambda^{c}}$ be the complement set containing indices corresponding to the coefficients in $\bm{C}_{j}$ having negligible magnitude ($<10\%$ signal energy). In a strictly sparse case, $\bm{\Lambda}$ will contain the indices corresponding to the non-zero coefficients and $\bm{\Lambda^{c}}$ will contain that corresponding to the zeros.

The algorithm identifies a matrix $\bm{A}$ such that the measurement vector of every signal in the training set has maximum entropy. The $j$-th measurement vector $\bm{Y}_{j}$ attains maximum entropy when its probability distribution of representation tends to be uniform; that is $|y_{ij}|\approx|y_{kj}|$ for $i\neq k$

Maximization of the entropy can lead to $|y_{ij}|\approx|y_{kj}|$ for $i\neq k$ with either all the $y_{ij}$’s being small or all the $y_{ij}$’s being large. This wide separation of values may occur if the rows of $\bm{A}$ scale $\bm{C}_{j}$ unduly. In the proposed algorithm, since $\bm{A}$ is generated using the orthogonal Procrustes method, the rows of $\bm{A}$ are orthonormal. Hence $\bm{A}$ does not drastically scale the coefficients in $\bm{C}_{j}$.

Since $\bm{A}$ does not unduly scale the coefficients in $\bm{C}_{j}$, the case with all the $y_{ij}$’s being small will occur when $\bm{A}$ captures the information from the coefficients present in $\bm{\Lambda^{c}}$, which do not contain important information of the signal. Hence this is an undesired case. The RIP penalty term ensures that the energy of the measurement vector $\bm{Y}_{j}$ does not deviate largely from the signal energy. Hence the case of all the $y_{ij}$’s being small is eliminated.

Hence, the entropy maximized measurements $\widehat{\bm{Y}}_{j}$ will capture maximum information from the coefficients of $\bm{C}_{j}$ corresponding to the indices in $\bm{\Lambda}$.

Remark: The entropy based measure of sparsity used in this paper depends on the $2$-norm of the signal. An entropy based sparsity measure depending on $p$-norm ($p\leq 1$) is proposed in [36]. Rao and Kreutz-Delgado [37] discuss the various criteria that the $p$-norm dependent entropy based measures of sparsity satisfy, for different values of $p$. The experimental results of [38] show that the freedom to choose suitable values of $p$ helps in fully exploiting the sparsity-promoting potential of the entropy function. Hence, the performance of the proposed EMS algorithm may improve on choosing an appropriate sparsity promoting value of $p$ instead of the $2$-norm. We leave the detailed discussion and analysis for future study.

4 Results

This section discusses the performance of the EMS matrix constructed using the proposed algorithm applied to a class of synthetic, speech, and image signals. The performance evaluation was done for noise-free signals and noisy signals. To measure the performance, the recovery methods used were the $l_{1}$ minimization or the Basis pursuit (BP) [1], the Entropy matching pursuit (EMP) [21], and the Orthogonal matching pursuit (OMP) [39] algorithms. The performance evaluation is done by calculating the signal to reconstruction error (SRER) which is given by

[TABLE]

where $\bm{x}$ is the original signal and $\hat{\bm{x}}$ is the recovered signal. Performance comparison is done for the measurements obtained from noise-free signals and noisy signals.

The performance of the EMS matrix constructed is compared with that of the random projection matrix, the optimized measurement matrix proposed by Elad [5], and the projection matrix generated by the NuMax [18][19]. The measurement matrix construction using Elad’s method and the NuMax algorithm were studied experimentally, using the softwares available in [40] and [41], respectively. The values of the parameters used for the construction of Elad’s measurement matrix are as mentioned in [5]. The training set used for the NuMax algorithm is the same as that used to train the EMS matrix.

4.1 Noise-free Signals

4.1.1 Synthetic signals

A set of signals, sparse in the Discrete Cosine Transform (DCT) basis with arbitrary support, was generated to test the performance of the EMS matrix. Two experiments were performed with the synthetic signals. For the first experiment, a set of $200$ signals of dimension $64$ and sparsity $10$ were generated. The measurement matrices for different values of $M$ were generated using the proposed algorithm for this class of synthetic signals. The BP algorithm was used to recover the signals from the reduced set of measurements. The average SRER plot for varying $M$ is shown in Fig 3(a).

The second experiment was to find the recovery performance when the number of measurements $M$ is fixed and the sparsity is varied. Signals of dimension $64$ and sparsity varying from $K=1$ to $30$ were constructed. A set of $200$ signals for each $K$ value was generated and concatenated to generate the training set (the training set hence contained 6000 signals). The EMS matrix was generated with $M=10$, $M=20$ and $M=30$ using these training signals. The variation of SRER against $K$ is plotted in Fig. 3(b)-(d).

The SRER plots show that the EMS matrix gives higher values of SRER even with less number of measurements, for strictly sparse signals, than attained using sensing matrices constructed using other methods. Fig. 3(b) shows that for a signal with $l_{0}$-sparsity $K=5$ and $M=10=2K$, the SRER is close to $13$dB when the measurements were made using the EMS matrix. Whereas the SRER is less than $5$dB when other measurement matrices were used for sensing. Similarly, Fig 3(c) shows that for a signal with $K=10$ the SRER is close to $26$dB when the EMS matrix with $M=2K=20$ measurements was used, which is approximately $8$dB greater than that of NuMax which gives the next best performance. Fig. 3(d) shows an improvement of about $50$dB over other measurement matrices for signal with $K=15$ and $M=30$. These observations confirm the claim in Lemma 3.

The constant $\delta$ in Eqn. (14) defines the upper bound on the radius of the recovery error sphere. It fixes the stability of a recovery algorithm which works on the measurement vector obtained using the $\bm{\Phi}$ matrix identified. The smaller the value of $\delta$ used in Eqn. (14), the better will be the recovery as validated in Fig. 4 which shows that as the value of delta increases, the SRER decreases.

The actual RIP constants obtained when the EMS matrix acts on a signal need not be the same as the one used in the algorithm. The values of $\delta$ obtained for the synthetic signals with sparsity 10 is given in the Table 1. The table gives the range of $\delta$ values obtained for the 200 signals.

Table 1 shows that though the obtained values of $\delta$ are not exactly the same as the value of $\delta$ assigned in the algorithm, the measurement matrix does satisfy the RIP. Further, as the number of measurements increases, the value of $\delta$ comes closer to that of the value assigned for $\delta$.

4.1.2 Speech signals

A set of 2450 speech signals of dimension $64$ from the database in the Linguistic Data Consortium for Indian Languages (LDC-IL) [42], sampled at 8kHz, was used as the training set. The measurement matrices for a set of $M$ values were constructed by applying the proposed learning algorithm to these training signals with the DCT as the sparsifying basis. In Fig. 5, a comparison of the SRER against $M$, for the test signals (signals belonging to the class but outside the training set) sampled using sensing matrices identified through various methods, is shown. In Fig. 5 the SRER curves of Elad’s measurement matrix and random measurement matrix are very close. The results in [5] show that the performance of the measurement matrix decreases as the $l_{0}$-sparsity $K$ increases. For speech signals the value of $K$ is not fixed and there may be signals in the test set with high values of $K$. This condition accounts for the low performance of the Elad’s measurement matrix.

We have compared the performance of the EMS matrix with the sensing matrix generated by applying MDS on a biorthogonal wavelet basis. For comparison, we identified the EMS matrix with the biorthogonal basis as the sparsifying basis $\bm{\Psi}$. The SRER plots are shown in Fig. 6.

The plots in Fig. 5 and Fig. 6 indicate that the sensing matrix generated leads to improved performance with the number of measurements smaller than that required by other measurement matrices. The improvement is more pronounced when EMP, an entropy based recovery algorithm, is used for recovery.

4.1.3 Image signals

The 64-dimensional test signals for images were constructed by considering non-overlapping $8\times 8$ blocks from various images. The number of image blocks used for training were 7725. The sparsifying basis used was the DCT basis. The proposed algorithm was used to construct the measurement matrix for this class of image signals. The SRER plots in Fig. 7 show the measurement matrix constructed is able to recover signals with less number of measurements than required by other measurement matrices.

To observe the perceptual quality of the signals, Figures 8 and 9 show the images measured using the EMS matrix of various $M$ values (the $M$ value is the number of measurements in each $8\times 8$ block) and reconstructed using BP and EMP algorithms, respectively. Fig. 8(b)-(d) give the images reconstructed when the measurement rates ($M/N$) are 0.1, 0.2 and 0.3, respectively. The Peak Signal to Noise Ratio (PSNR) values show that the recovery error is small for measurement rate close to 0.3. Fig. 10 gives a comparison of the reconstructed images when the images were measured using various measurement matrices with $M=30$, for each $8\times 8$ block, and BP was used as the recovery algorithm. It can be seen that the perceptual quality of the reconstructed image when the measurement was taken using the EMS matrix is better than that of the reconstructed images when the measurements were taken using other sensing matrices.

Fig. 11 shows images measured using the EMS matrix with $M=30$, for each $8\times 8$ block, and reconstructed using the BP algorithm. The PSNR values show that the reconstructed images retain most of the information contained in the original image and the information loss is minimal. The figure establishes that the EMS matrix is capable of capturing maximum information from any image signal.

4.2 Noisy Signals

The measurement of noisy signal $\bm{x_{e}}=\bm{x}+\bm{e}$, where $\bm{e}$ is the additive white noise, can be expressed as

[TABLE]

where, $\bm{c}$ is the representation of the signal $\bm{x}$ in the basis $\bm{\Psi}$ and $\bm{c_{e}}$ is the representation of the noise $\bm{e}$ in the basis $\bm{\Psi}$. We know that the representation $\bm{c}$ of the signal $\bm{x}$ in $\bm{\Psi}$ is sparse but the representation $\bm{c_{e}}$ of the noise $\bm{e}$ is dense in $\bm{\Psi}$. The matrix $\bm{A}$ constructed using the proposed method captures maximum information from the sparse set of coefficients. Ideally, the measurement matrix $\bm{A}$ will thus capture only the noisy components having support in $\bm{\Lambda}$ and reject the noisy component having support in $\bm{\Lambda^{c}}$. Also, the measurement matrix is constructed such that the RIP is satisfied. The RIP ensures that the measurement matrix-recovery algorithm pair is stable. Hence, stable recovery is ensured when the EMS matrix is paired with any recovery algorithm.

To validate the performance in the presence of noise, signals contaminated with white Gaussian noise were measured using the sensing matrices and recovered using Basis Pursuit De-Noising (BPDN), EMP and OMP algorithms. The SRER plots of signals sparse in DCT basis, with sparsity 10, for input SNR 3dB (Fig. 12) show that the performance of the EMS matrix paired with any recovery algorithm is stable.

For noisy signals, the total error $(x-\hat{x})$ is contributed by the recovery error and the noise. The low SRER with less number of measurements is due to the high recovery error and the noise. As the number of measurements increases, the recovery error decreases and the major contribution to the total error is the noise. The noise reduction property of the recovery methods accounts for the marginal increase in the SRER of the recovered signals, above the input SNR, at large values of $M$. The noise reduction is more pronounced in OMP because the OMP algorithm was run for 10 iterations (since $K=10$), eliminating the noise component in the rest of the coefficients.

5 Conclusion

We have presented a learning method for constructing an efficient sensing matrix for the compressive sensing of a class of signals without assuming structured sparsity. The construction of the measurement matrix employs a learning scheme that maximizes the entropy of the measurement vectors of a set of training signals. We have established the bounds on the entropy of measurements necessary for the unique recovery of a signal. We have also proved that maximization of the entropy of measurements leads to a reduction in the number of measurements required for the unique recovery of signals. The sensing matrix designed was used for the compressive measurements of a class of sparse synthetic signals of arbitrary support, speech signals and image signals. The recovery of the signals is significantly better, with less number of measurements, than the recovery from the measurements obtained using other existing sensing matrices.

Appendix

A. Necessary condition for unique recovery is $M_{e}ff\geq 2K$.

Let $\bm{y}$ be the measurement vector corresponding to the sparse representation vector $\bm{c}$ of a compressible signal obtained as $\bm{y}=\bm{Ac}$, where the matrix $\bm{A}$ has orthogonal rows. Consider the case where $\bm{y}$ is strictly sparse with $\|\bm{y}\|_{0}=M_{eff}<M$. Here, $\bm{y}=\bm{Ac}=\tilde{\bm{A}}\bm{c}$, where $\tilde{\bm{A}}$ is formed by replacing the rows of $\bm{A}$ corresponding to the zero entries in $\bm{y}$ with rows of zeros (there would be $M-M_{eff}$ zero rows in $\tilde{\bm{A}}$). Thus,

[TABLE]

According to the theory of CS, the condition for unique recovery is $\text{spark}(\tilde{\bm{A}})>2K$ [6].

[TABLE]

Thus we can say that the condition for unique recovery is $M_{eff}\geq 2K$. The argument can be extended to the case of a compressible $\bm{y}$ where $n_{th,y}^{I}=\lceil\exp(H(\bm{y}))\rceil=M_{eff}$.

B. Solution to Orthogonal Procrustes Problem with rectangular matrix.

If the singular value decomposition of $\bm{C}\widehat{\bm{Y}}^{T}$ is given by $\bm{U\Delta V}^{T}$, where $\bm{U}\in\mathbb{R}^{N\times N}$ and $\bm{V}\in\mathbb{R}^{M\times M}$, we propose that the solution to the problem (15) is given by $\widehat{\bm{A}}=\bm{VU_{M}}^{T}$, where $U_{M}$ contains the first $M$ columns of $U$ which correspond to the largest $M$ singular values of $C\widehat{Y}^{T}$.

Proof.

We know, $\widehat{\bm{Y}}\in\mathbb{R}^{M\times L}$, $\bm{A}\in\mathbb{R}^{M\times N}$ and $\bm{C}\in\mathbb{R}^{N\times L}$. Assuming $\bm{A}$ to have orthogonal rows we get

[TABLE]

where $tr(\bm{A})$ refers to trace of the matrix $\bm{A}$. Minimizing $\|\widehat{\bm{Y}}-\bm{AC}\|_{F}^{2}$ is equivalent to maximizing $tr(\widehat{\bm{Y}}\bm{C}^{T}\bm{A}^{T})$

[TABLE]

where $\bm{C}\widehat{\bm{Y}}^{T}=\bm{U\Delta V}^{T}$ and $\bm{U}\in\mathbb{R}^{N\times N}$, $\bm{\Delta}\in\mathbb{R}^{N\times M}$ and $\bm{V}\in\mathbb{R}^{M\times M}$. Since $\bm{C}\widehat{\bm{Y}}^{T}$ is a rectangular matrix with $M<N$, the $\bm{\Delta}$ matrix is diagonal with $N-M$ rows equal to zero. Thus, $\bm{U\Delta}=\bm{U_{M}\Delta_{M}}^{T}$, where $\bm{U_{M}}$ contains the first $M$ columns of $\bm{U}$ which correspond to the largest $M$ singular values of $\bm{C}\widehat{\bm{Y}}^{T}$, and $\bm{\Delta_{M}}^{T}$ contains the rows of $\bm{\Delta}$ containing the non-zero $M$ singular values of $\bm{C}\widehat{\bm{Y}}^{T}$ ($\bm{\Delta_{M}}^{T}$ is a square diagonal matrix).

[TABLE]

This trace will be maximum when $\bm{V}^{T}\bm{AU_{M}}=\bm{I}$ [35]. Hence, $\bm{A}=\bm{VU_{M}}^{T}$ with orthogonal rows. ∎