Iterative method for simultaneous sparse approximation
Sahar Sadrizadeh Shahrzad Kiani Mahdi Boloursaz Farokh Marvasti
cs.IT math.IT







Journal of  Class Files, Vol. 14, No. 8, August 2015
S. Sadrizadeh and S. Kiani : Simultaneous Sparse Approximation

Iterative method for simultaneous sparse approximation
    Sahar Sadrizadeh,
        Shahrzad Kiani,
        Mahdi Boloursaz, Student Members, IEEE,
        and Farokh Marvasti, Senior Member, IEEE
S. Sadrizadeh, S. Kiani, M. Boloursaz and F. Marvasti are with the Advanced Communication Research
Institute (ACRI), Electrical Engineering Department, Sharif University of Technology,
Tehran, Iran (email: sadrizadeh_s@ee.sharif.edu; shkianid@gmail.com; marvasti@sharif.edu).
    December 30, 2023
==================================================================================================================================================================================================================================================================================================================================================================================================================================





This paper studies the problem of Simultaneous Sparse Approximation (SSA). This problem arises in many applications which work with multiple signals maintaining some degree of dependency such as radar and sensor networks. In this paper, we introduce a new method towards joint recovery of several independent sparse signals with the same support.  We provide an analytical discussion on the convergence of our method called Simultaneous Iterative Method with Adaptive Thresholding (SIMAT). Additionally, we compare our method with other group-sparse reconstruction techniques, i.e., Simultaneous Orthogonal Matching Pursuit (SOMP), and Block Iterative Method with Adaptive Thresholding (BIMAT) through numerical experiments.  The simulation results demonstrate that SIMAT outperforms these algorithms in terms of the metrics Signal to Noise Ratio (SNR) and Success Rate (SR). Moreover, SIMAT is considerably less complicated than BIMAT, which makes it feasible for practical applications such as implementation in MIMO radar systems.




Simultaneous Sparse Approximation; Iterative Method; Adaptive Thresholding; Joint Recovery.





§ INTRODUCTION

Sparse
 signal processing has recently been exploited in various fields of communication, due to fact that sparse signals can be approximated by only a few nonzero coefficients and hence sub-Nyquist sampling and Compressed Sensing (CS) <cit.>. The general CS problem is formulated as follows:

|l|
  𝐱_0
  
  𝐲=𝐀𝐱+𝐯

where x is the main sparse signal, y is the measurement vector, A is the sensing matrix, and v is the additive noise vector. Two main models are considered in CS for reconstruction of sparse signals. Models with one measurement vector are referred to as Single Measurement Vector (SMV) models, while the other models with at least two measurement vectors are called Multiple Measurement Vector (MMV) models.


    The problem investigated in MMV models, known as SSA, aims to jointly recover sparse representation of the measurement vectors. The SSA applications may be encountered in various fields such as sensor networks <cit.>, Electroencephalography and Magnetoencephalography (EEG and MEG) <cit.>, source localization <cit.>, and distributed MIMO radar systems <cit.>.


    <cit.> investigates the theory of MMV models. Some algorithms have been developed by extending the general SMV model into the MMV model to solve the SSA problems. Orthogonal Matching Pursuit (OMP)<cit.> as a greedy algorithm is one of the very first algorithms used for sparse recovery. At each iteration of this algorithm, the best local improvement to the current approximations is found in hope of obtaining a good overall solution. The extension of the OMP algorithm to The MMV paradigm, Simultaneous OMP or SOMP, has been developed in <cit.>.


    The Iterative Method with Adaptive Thresholding (IMAT) algorithm was originally proposed for sparse signal reconstruction from missing samples <cit.>. The Block Iterative Method with Adaptive Thresholding (BIMAT) <cit.> as an extension of IMAT is employed for block sparse recovery for distributed MIMO radar systems.

    In this paper, we propose SIMAT as a generalization of IMAT for simultaneous reconstruction of jointly sparse signals from their missing samples.




 §.§ Paper Overview 

    The rest of this paper is structured as follows. In Section <ref>, we first provide the description of SSA model. Then the proposed method is introduced and its convergence is analyzed. Numerical experiments of our method in comparison with the SOMP algorithm are presented in Section <ref>. SIMAT is then demonstrated as a simple decoding algorithm for MIMO radar systems, and its performance is compared with BIMAT by means of simulation. Finally, the paper is concluded in Section <ref>.




 §.§ Notations

Scalar variables, vectors, and matrices are denoted by italic lower-case, boldface lower-case, and boldface upper-case, respectively. The elements of a vector are denoted by subscript, i.e., x_i is the i-th element of the vector 𝐱. |𝐱| calculates the absolute value of each entries of the vector 𝐱. The pseudoinverse of matrix 𝐀 is represented by 𝐀^†.  Finally, the output of the thresholding operator TH( 𝐱,thr) is defined as a diagonal matrix whose diagonal entries are determined as follows:



    TH( 𝐱,thr)_ii = {[ 1.       | x_i| ≥ thr; 0.       | x_i| < thr ].





§ THE PROPOSED METHOD





 §.§ Problem Statement

    In this section we provide the formulation of SSA problems. Assume that 𝐱^1, 𝐱^2, ..., 𝐱^L are L equal-length signals that share the same sparsity support in a specific transform domain. These signals are random-sampled in another domain by independent sampling masks 𝐬^1, 𝐬^2, ..., 𝐬^L. These masks are binary and each element of them is generated independently based on a Bernoulli distribution , i.e., 𝐬^i_n∼ Bernoulli(p), 1≤ i ≤ L, and 0≤ p ≤ 1 is the sampling probability. The random-sampled signals 𝐲^1, 𝐲^2, ..., 𝐲^L are derived as follows:



    y^i_n = s^i_n ×x^i_n



   The problem is to simultaneously reconstruct the original sparse signals from their random-sampled versions by employing the additional information that the signals share a common support in a specific transform domain.




 §.§ Algorithm

One can find the SIMAT algorithm in Algorithm <ref>. Let n , m and k denote the lengths of the original signal, the length of the observed vectors, and the sparsity number, respectively. Moreover, let 𝐱^j_i represent the reconstruction of the j-th signal after i-th iteration.


    It should be noted that the measurement matrices, i.e., [𝐀^1, 𝐀^2, ..., 𝐀^L], can be calculated by multiplying the transformation matrix, which maps the signals to their sparsity domain, and diagonal matrices whose diagonal entries are equal to the elements of  𝐬^1, 𝐬^2, ..., 𝐬^L.


    This algorithm gradually extracts the sparse components of the signals by thresholding the estimated signals iteratively. Each iteration involves two different steps of thresholding and projection. The thresholding step provides an approximation of the common support of the signals by hard-thresholding the summation of the absolute values of the approximated signals. The projection step projects each of the estimated signals onto the convex set defined by the support vector approximated in the previous step.


    In this algorithm λ  is the relaxation parameter and controls the convergence speed. The threshold value is decreased exponentially by βe^ - α( k - 1), where k is the iteration number. The algorithm performance is less dependent on the parameters λ, β,  and α; however, these parameters are optimized empirically to achieve a faster convergence.





 §.§ Analytical Discussion

    In <cit.>, it is proved that the IMAT method converges to the sparset solution of the random sample SMV problem if the threshold value is less than
    ε^2/k, where k is the sparsity number, and ε is the minimum distance between the subspace induced by the random samples and all the subspaces induced by a specific support which do not intersect with the previously mentioned subspace. Now, we show that the probability of finding the support of signals with the SIMAT is more than that of the IMAT. Before proving this statement, we note the following points.

    We assume that the original signals are sparse in the time domain, and each non-zero coefficient has a Gaussian distribution. Additionally, the signals are random-sampled and polluted by additive white Gaussian noise in the frequency domain. Since the signals are estimated by line 16 of Algorithm <ref> in each iteration, each coefficient of the estimated signals has a Gaussian distribution.

    The distribution of the absolute value of  a random variable with Gaussian distribution N(0,σ^2) is Half-Normal with the following Probability Density (PDF) and Cumulative Distribution (CDF):



    P_Y( y;σ) =  √(2)/√(σ ^2×π)exp(  - y^2/2 ×σ ^2)    y ≥ 0



    Q_Y( y;σ) = erf( y/√(σ ^2× 2))



In the thresholding step of each iteration of our algorithm, the absolute values of the coefficients are calculated. Hence, these variables have Half-Normal distribution. We define the variable z as the absolute value of the estimation of a non-zero coefficient of a signal and assume that it has variance σ ^2_1 and mean μ _1. Similarly, we define the variable w associated with a zero-coefficient of a signal and assume that its variance and mean are σ ^2_0 and μ _0, respectively. We also assume that μ _1 > μ _0 since the variance of the additive noise is small. Due to the fact that the mean of L signals is calculated for the SIMAT, the variables z_L and w_L can be defined as the mean of the variables mentioned above.  The variance and the mean of z_L and w_L are σ ^2_1/L, μ _1,
σ ^2_0/L and μ _0, as a result of the independency of signals.


    We prove the superiority of SIMAT over IMAT through Lemma <ref> and Lemma <ref>. As a general idea of these lemmas, according to the law of large numbers, the coefficients of the signal in the SIMAT become nearer to their mean when L, the number of signals, go to infinity. In fact, the none-zero coefficients tend to μ _1 and the zero coefficients tend to μ _0. Therefore, it would be easier to find the non-zero coefficients of the original signals.



If the number of signals L satisfies the following inequality

    L ≥σ ^2_0/( ε  - μ _0)( 1 - erf( ε/√(σ ^2_0 × 2))) > 0,

then P( w ≥ε) ≥ P( w_L≥ε), ∀ε≥ 2 ×μ _0, where P is the probability sign. w and w_L are the variables representing the estimated zero coefficients of a signal in each iteration of IMAT and SIMAT, respectively.



As a result of ε≥ 2 ×μ _0 and positivity of w_L, one can easily show the following equality:


    P( w_L≥ε) =  P( | w_L - μ _0 | ≥ε  - μ _0)


According to the Chebyshev's theorem, we have


    P( | w_L - μ _0| ≥ε  - μ _0) ≤σ ^2_0/L/ε  - μ _0


If we define Q_w( ε) as the CDF of P( w ), we get


    P( w ≥ε) = 1 - Q_w( ε) = 1 - erf( ε/√(σ ^2_0 × 2))


Since 1 - erf( ε/√(σ ^2_0 × 2)) ≥ 0 and ( ε  - μ _0) ≥ 0, the lemma is proved for the number of signals specified by (<ref>).





If the number of signals L satisfies the following inequality

    L ≥σ ^2_1/( μ _1 - ε)( erf( ε/√(σ ^2_1 × 2))) > 0,

then P( z ≥ε) ≤ P( z_L≥ε), ∀  μ _1≥ε≥ 2 ×μ _0, where P is the probability sign, and the estimated non-zero coefficients of a signal in each iteration of IMAT and SIMAT are denoted by z and z_L, respectively.



It is not too difficult to derive the following inequality


    P( z_L≥ε) ≥ 1 -  P( | z_L - μ _1| ≥μ _1 - ε)


According to the Chebyshev's theorem, we have


    P( | z_L - μ _1| ≥μ _1 - ε) ≤σ ^2_1/L/μ _1 - ε⇒
        P( z_L≥ε) ≥ 1 -  σ ^2_1/L/μ _1 - ε


If we define Q_z( ε) as the CDF of P(z), then we get


    P( z ≥ε) = 1 - Q_z( ε) = 1 - erf( ε/√(σ ^2_1 × 2))


Since erf( ε/√(σ ^2_1 × 2)) ≥ 0 and ( μ _1 - ε) ≥ 0 for μ _1≥ε≥ 2 ×μ _0, one can easily prove the lemma for the L indicated in (<ref>)




If the number of signals L satisfies the following inequality

    ∀   μ _1   ≥ε≥ 2 ×μ _0    , 
        L ≥max(σ ^2_1/( μ _1 - ε)( erf( ε/√(σ ^2_1 × 2)))
        ,σ ^2_0/( ε  - μ _0)( 1 - erf( ε/√(σ ^2_0 × 2)))),

then ∀  μ _1≥ε≥ 2 ×μ _0, the probability of finding the support with SIMAT is higher than IMAT.



According to Lemma <ref>, if L satisfies (<ref>), the probability of mistaking a zero coefficient for a support in SIMAT is less than that of IMAT.
According to Lemma <ref>, if L satisfies (<ref>), the probability of finding a none-zero coefficients in SIMAT is higher than IMAT.

Therefore, if L satisfies (<ref>), it is more probable to find the support with the SIMAT than with the IMAT.


One can find the minimum value of L which satisfies (<ref>) as shown in Table <ref>. This table shows that the minimum number of signals L that guarantees superiority of SIMAT over IMAT is very low, and hence there is no need to have a large number of signals to benefit from SIMAT.







§ SIMULATION RESULTS













In each trial of our simulation, we generate L number of K-sparse signals. We choose K components out ofN=256 randomly, and set them to a random number in the interval [-1,1]. Then the noisy signals are random-sampled by a sampling rate of M/N. By the law of algebra, the number of samples needed to specify the sparsity profile of the signals is at least twice the sparsity number, hence K ≤ M/((2× N)). We optimize the parameters of the algorithm in each trial.

Table <ref> compares the average reconstruction SNR (dB) of three algorithms, IMAT, SIMAT and SOMP, for different density numbers and input SNRs. As observed in Table <ref>, the simultaneous reconstruction methods outperform IMAT, especially for higher sparsity numbers.


In the case of not knowing the density rate, SIMAT yields the best results both in low and high density rate, and in noisier channels. In noiseless channels, the signal can be perfectly reconstructed by SIMAT, when the sparsity number is small.


In the case of knowing the density rate, SOMP and SIMAT exhibit similar performances in all the cases. However, based on the results of Table <ref>, SIMAT outperforms SOMP in terms of the complexity measured by the run-time.

 The success rate of IMAT, SIMAT and SOMP algorithms for different density rates is depicted in Fig. <ref>. A reconstruction is considered to be successful if the output SNR is more than 20dB. As seen in this figure, all curves experience a sudden knee-like fall as the density rate increases. This fall is considered as the boundary between successful and unsuccessful reconstruction. The simulation results reveal that for the SIMAT, success rate falls around 20% density rate. While, the knee-like fall happens in 18% and 12% density rate for IMAT and SOMP, respectively. This indicates that SIMAT can successfully reconstruct the signals with higher sparsity number in comparison with SOMP and IMAT algorithms.

Fig. <ref> shows the effect of the number of signals on the performance of SOMP and SIMAT.  As observed in this figure, by increasing the number of signals, we get better results. Additionally,
the SIMAT algorithm can reconstruct the signals better than the SOMP when the number of signals is small.


  The reconstruction SNR for different sampling rates is depicted in Fig. <ref>. According to this figure, the SNR values can be improved by increasing the sampling rate. Although the SOMP algorithm has better results for a small interval in the high sampling rates, SIMAT outperforms SOMP in most cases.

   


   Table <ref> lists the SNR of two reconstruction methods extended from IMAT. Simulation results demonstrate the superiority of SIMAT over BIMAT in terms of the output SNR and the complexity. Therefore, the SIMAT can be used instead of the BIMAT in many applications such as distributed MIMO radar systems.





§ CONCLUSION

In this paper a novel method, SIMAT, was introduced for SSA problems. The proposed method is an extension of IMAT into the MMV models, and the idea of this extension is that the summation of sparse vectors, sharing the same support, in the thesholding step of each iteration can enhance the probability of reconstruction. Indeed, the theoretical analysis with simulation results prove that SIMAT outperforms IMAT with respect to the SNR metric. We compared the proposed method with SOMP, as a well-known algorithm in the MMV models. We can conclude from the conducted numerical experiments that SIMAT is preferable in terms of SNR or SR, specially in noisier channels with low sampling rates and high density rates. Finally, it was observed that SIMAT is superior to BIMAT when complexity and efficiency are important factors.



§ ACKNOWLEDGMENT

The authors would like to thank Ehsan Asadi for helping in the analytical discussion subsection and Babak Barazandeh for introducing this topic to one of the authors.


10

url@samestyle








candes2006robust
E. J. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles:
  Exact signal reconstruction from highly incomplete frequency information,”
  IEEE Transactions on information theory, vol. 52, no. 2, pp. 489–509,
  2006.

donoho2006compressed
D. L. Donoho, “Compressed sensing,” IEEE Transactions on information
  theory, vol. 52, no. 4, pp. 1289–1306, 2006.

duarte2009model
M. F. Duarte, V. Cevher, and R. G. Baraniuk, “Model-based compressive sensing
  for signal ensembles,” in Communication, Control, and Computing, 2009.
  Allerton 2009. 47th Annual Allerton Conference on.1em plus 0.5em
  minus 0.4emIEEE, 2009, pp. 244–250.

gastpar2003source
M. Gastpar and M. Vetterli, “Source-channel communication in sensor
  networks,” in Information Processing in Sensor Networks.1em
  plus 0.5em minus 0.4emSpringer, 2003, pp. 553–553.

gorodnitsky1995neuromagnetic
I. F. Gorodnitsky, J. S. George, and B. D. Rao, “Neuromagnetic source imaging
  with focuss: a recursive weighted minimum norm algorithm,”
  Electroencephalography and clinical Neurophysiology, vol. 95, no. 4,
  pp. 231–251, 1995.

malioutov2005sparse
D. Malioutov, M. Cetin, and A. S. Willsky, “A sparse signal reconstruction
  perspective for source localization with sensor arrays,” IEEE
  transactions on signal processing, vol. 53, no. 8, pp. 3010–3022, 2005.

ibernon2008comparison
R. Ibernon-Fernandez, J.-M. Molina-Garcia-Pardo, and L. Juan-Llacer,
  “Comparison between measurements and simulations of conventional and
  distributed mimo system,” IEEE Antennas and Wireless Propagation
  Letters, vol. 7, pp. 546–549, 2008.

mohammadi2017sampling
E. Mohammadi, A. Fallah, and F. Marvasti, “Sampling and distortion tradeoffs
  for indirect source retrieval.”1em plus 0.5em minus 0.4emIEEE,
  2017, accepted for publication.

tropp2007signal
J. A. Tropp and A. C. Gilbert, “Signal recovery from random measurements via
  orthogonal matching pursuit,” IEEE Transactions on information
  theory, vol. 53, no. 12, pp. 4655–4666, 2007.

determe2017noise
J.-F. Determe, J. Louveaux, L. Jacques, and F. Horlin, “On the noise
  robustness of simultaneous orthogonal matching pursuit,” IEEE
  transactions on signal processing, vol. 65, no. 4, pp. 864–875, 2017.

tropp2006algorithms
J. A. Tropp, A. C. Gilbert, and M. J. Strauss, “Algorithms for simultaneous
  sparse approximation. part i: Greedy pursuit,” Signal Processing,
  vol. 86, no. 3, pp. 572–588, 2006.

rakotomamonjy2011surveying
A. Rakotomamonjy, “Surveying and comparing simultaneous sparse approximation
  (or group-lasso) algorithms,” Signal processing, vol. 91, no. 7, pp.
  1505–1526, 2011.

marvasti2012sparse
F. Marvasti, M. Azghani, P. Imani, P. Pakrouh, S. J. Heydari, A. Golmohammadi,
  A. Kazerouni, and M. Khalili, “Sparse signal processing using iterative
  method with adaptive thresholding (imat),” in Telecommunications
  (ICT), 2012 19th International Conference on.1em plus 0.5em minus
  0.4emIEEE, 2012, pp. 1–6.

marvasti2012unified
F. Marvasti, A. Amini, F. Haddadi, M. Soltanolkotabi, B. H. Khalaj,
  A. Aldroubi, S. Sanei, and J. Chambers, “A unified approach to sparse signal
  processing,” EURASIP journal on advances in signal processing, vol.
  2012, no. 1, p. 44, 2012.

marvasti2012nonuniform
F. Marvasti, Nonuniform sampling: theory and practice.1em plus
  0.5em minus 0.4emSpringer Science & Business Media, 2012.

abtahi2016iterative
A. Abtahi, M. Azghani, J. Tayefi, and F. Marvasti, “Iterative block-sparse
  recovery method for distributed mimo radar,” in Communication and
  Information Theory (IWCIT), 2016 Iran Workshop on.1em plus 0.5em
  minus 0.4emIEEE, 2016, pp. 1–4.

esmaeili2016iterative
A. Esmaeili, E. Asadi, and F. Marvasti, “Iterative null-space projection
  method with adaptive thresholding in sparse signal recovery and matrix
  completion,” arXiv preprint arXiv:1610.00287, 2016.



