Online Wideband Spectrum Sensing Using Sparsity
Lampros Flokas Petros Maragos
cs.IT math.IT





































































Online Wideband Spectrum Sensing Using Sparsity
    Lampros Flokas
        and Petros Maragos, Fellow, IEEE
Lampros Flokas is with the Columbia University in the City of New York, Department of Computer Science, New York, USA, .Petros Maragos is with the National Tech. University of Athens, School of ECE, Greece, Email: 






    December 30, 2023
========================================================================================================================================================================================================================================================================================






Wideband spectrum sensing is an essential part of cognitive radio systems. Exact spectrum estimation is usually inefficient as it requires sampling rates at or above the Nyquist rate. Using prior information on the structure of the signal could allow near exact reconstruction at much lower sampling rates. Sparsity of the sampled signal in the frequency domain is one of the popular priors studied for cognitive radio applications. Reconstruction of signals under sparsity assumptions has been studied rigorously by researchers in the field of Compressed Sensing (CS). CS algorithms that operate on batches of samples are known to be robust but can be computationally costly, making them unsuitable for cheap low power cognitive radio devices that require spectrum sensing in real time. On the other hand, online algorithms that are based on variations of the Least Mean Squares (LMS) algorithm have very simple updates so they are computationally efficient and can easily adapt in real time to changes of the underlying spectrum. In this paper we will present two variations of the LMS algorithm that enforce sparsity in the estimated spectrum given an upper bound on the number of non-zero coefficients. Assuming that the number of non-zero elements in the spectrum is known we show that under conditions the hard threshold operation can only reduce the error of our estimation. We will also show that we can estimate the number of non-zero elements of the spectrum at each iteration based on our online estimations. Finally, we numerically compare our algorithm with other online sparsity-inducing algorithms in the literature.  




signal processing, sparse representations, LMS, cognitive radio.





















§ INTRODUCTION
















Wireless telecommunications spectrum is a limited resource and with the rapid increase of telecommunication applications, static allocation of spectrum for each case is not a viable solution. Additionally, static allocation of spectrum is also not effective as the primary users of the spectrum may use it from time to time and only in some locations. To overcome this limitation cognitive radio devices try to dynamically manage the spectrum by detecting which part of the spectrum is unused by its primary users and temporarily using it for their own needs. 
In order to be effective, these devices would need to check a wide band of frequencies to increase the possibility of finding unused spectrum. If cognitive devices used sampling rates that are equal or above the Nyquist rate, their cost would be prohibitive for most applications. In order to reduce the sampling rate needed as well as the computational effort, we will need to use some prior information on the structure of the received signal. This prior is that the same one that enables the usage of cognitive radio devices in the first place: Primary users do not use their share of the spectrum all the time so the received signal should be sparse in the frequency domain.
The area of compressed sensing (CS) has provided several celebrated algorithms for the reconstruction of undersampled signals with sparse representations <cit.>. Classic algorithms of CS assume a batch setting where the device is assumed to collect a number of observations and operate on them  in an iterative manner. Therefore it is of great importance to provide algorithms that reduce the number of iterations needed in order to reduce the computational burden on the cognitive radio devices and provide real time spectrum estimations. CS based approaches have been adapted by many researchers in the area of spectrum sensing for cognitive radio applications <cit.>.
On the other hand, online algorithms, based on variations of Least Mean Squares introduced by Widrow and Hoff <cit.>, have also been adapted for the CS setting. Algorithms like the ones presented in <cit.> have been shown to estimate sparse signals with faster convergence and smaller steady state errors than methods that do not exploit sparsity. Additionally, they have much simpler updates based on a single sample at a time. This allows them not only to be more computationally efficient but also to be adaptive to the changes of the estimated signal.   
Here we will propose two new variations of the classical LMS algorithm. The first is a variation of the Zero Attracting LMS <cit.> that does not penalize the s algebraically largest coeffients of the estimation where s is an upper bound on the number of non-zero elements in the estimated vector. The second one alternates the standard LMS update with shrinkage using a hard threshold operator. The hard threshold operator will keep the s algebraicly largest components of the estimated vector where s is again an upper bound on the number of non-zero elements in the estimated vector. This algorithm is the online version of the iterative hard thresholding studied in <cit.> and <cit.> and <cit.>. The sparsity of the estimated vector or even an upper bound on it may not be known in advance so we will also propose a way to estimate it in an adaptive manner. Even though we are going to apply the proposed algorithms for the problem of spectrum estimation, they can also be applied in other telecommunications and general machine learning applications where the incoming signals have a known sparse representation.  
The structure of the paper is as follows. In Section 2 we will define the problem of sparse spectrum reconstruction using below Nyquist rate sampling frequencies. In Section 3 we will present the properties of online sparsity aware estimation techniques in the literature and in Section 4 we will introduce our hard thresholding based algorithms. In Section 5 numerical simulations comparing our algorithm with other sparsity aware algorithms are provided. Finally, Section 6 contains concluding remarks and discusses possible directions for future research.


§ PROBLEM STATEMENT

Let 𝐳∈ℝ^N be the full signal that the cognitive radio device would receive if it was sampling it at the Nyquist rate. We would like to undersample 𝐳, taking just M observations from 𝐳 where M<N. Let us call 𝐔 the undersampling matrix whose rows are a subset of the rows of the identity matrix including only the rows where the respective observation of 𝐳 is sampled. Let us call 𝐲=[y_0,y_1,…,y_M-1]^T  ∈ℝ^M the resulting vector. If each observation y_i is corrupted by an additive error term v_i and 𝐯=[v_0,v_1,…,v_M-1]^T ∈ℝ^M Then we obviously have that 

    𝐲=𝐔𝐳 + 𝐯

Of course without any further assumptions the lost information cannot be recovered and important information about the spectrum of 𝐳 cannot be estimated. However, in our case we can assume that the underlying spectrum of the signal is sparse as a large portion of the spectrum  will be left unused by its primary users. Let  𝐰∈ℂ^N be the complex vector representing the Discrete Fourier Transform (DFT) of 𝐳 and Φ be the Inverse Discrete Fourier Transform (IDFT) matrix so that 𝐳=Φ𝐰. Given our assumption on the sparsity of the spectrum of 𝐳, we have that 𝐰 is a sparse vector and therefore we are interested in solving the following problem:

    min𝐰_0: 𝐲-(𝐔Φ)𝐰_2 ≤δ

where the ℓ_0 norm is the count of non-zero elements of the vector i.e 𝐰_0=|support(𝐰)|, where support(𝐰)={ i ∈{0,1,..,N-1}: w_i ≠ 0},  |S| denotes the cardinality of set S and δ is an upper bound on 𝐯_2. In general this problem is NP-hard and therefore computationally intractable. However, researchers in the area of CS have developed several algorithms that recover the solution of problem described by <ref> when the matrix 𝐔Φ satisfies the Restricted Isometry Property and vector 𝐰 is sparse enough. Out of all the algorithms probably the most popular one is Lasso regression. One of its equivalent formulations is

    min𝐰_1: 𝐲-(𝐔Φ)𝐰_2 ≤δ
 
The resulting optimization problem can be solved with standard convex optimization methods. The limiting factor for Lasso and other classical approaches to CS is that they may require several iterations to converge to the optimal solution. This makes them unsuitable for the low power cognitive radio devices that need real time spectrum estimations in potentially highly volatile settings.
In contrast, online estimation algorithms have much simpler update rules that involve one sample at a time and are robust to changes in the estimated signal. In the online setting there is a stream of measurements of 𝐲 and the corresponding rows of Φ that are fed one by one to the online algorithm. There are at least two options when it comes to forming this stream.


  * The first option is to use an online algorithm as a drop in replacement of a batch algorithm. We can collect M out of N samples of the signal and feed them one by one to the online algorithm. Of course the online algorithm may not converge in a single pass over the data so we can augment the stream by feeding the same measurements  to the algorithm multiple times in order to achieve convergence. The advantage of the online algorithms over batch algorithms in this setting is that they have simpler update rules than their batch counterparts and so they could be more easily implementable in low power cognitive radio devices.  

  * The second option is to form a stream by continuously incorporating new measurements. One way to do this is to split the incoming signal in non overlapping windows of length N, randomly sample M measurements in each window and feed the resulting measurements to the online algorithm. The advantage of the online algorithms over batch algorithms in this setting is that they can track the spectrum changes in the signal in real time.

In Section 5 we shall provide experimental results for both settings.   


§ RELATED WORK
 


 §.§ The LMS algorithm

The algorithms proposed in this paper are based on the updates of the LMS algorithm. To better understand the procedure we review the steps of the classical LMS algorithm. Let y(n) be a sequence of observations of the output of a system following the model

    y(n)=𝐰^H𝐱(n)+v(n)

where 𝐰=[w_0,w_1,…,w_N-1]^T ∈ℂ^N  is the parameter vector to be estimated, 𝐱(n)∈ℂ^N is taken from the rows of the Φ^* that correspond to the observed samples and v(n) is the additive observation noise. Let also 𝐰(n) be the estimation we have up to time n for the unknown vector 𝐰 and e(n) be the sample error. Then

    e(n)=y(n)-𝐰^H(n)𝐱(n)

The LMS update rule is recursive and produces a new estimation given the previous one, following the rule

    𝐰(n+1)=𝐰(n)+μ e^*(n)𝐱(n)

where μ is a an appropriately chosen constant. If 𝐑_x=𝔼[𝐱(n)𝐱^H(n)] is the uncentered covariance matrix of 𝐱(n), here assumed constant over time, and λ_max is its maximum eigenvalue then <cit.> shows that LMS will converge in the mean sense if:

    0<μ<2/λ_max

Of course the simple LMS algorithm has the same guarantees for all estimated signals, sparse and dense alike. Using the sparsity assumption can increase the speed of convergence and yield much lower steady state estimation errors than the classical algorithms.


 §.§ Zero Attracting LMS

The Zero Attracting LMS algorithm (ZA-LMS) <cit.> is a modification of the standard LMS algorithm that specializes in sparse system identification.  This algorithm follows the spirit of the equivalence of the ℓ_1 and ℓ_0 regularization problems in the batch case. Therefore the objective minimized at every iteration becomes

    J_ZA(n) = 1/2e(n)^2 +γ𝐰(n)_1

for some parameter γ. Taking the gradient descent update one can adapt the LMS update scheme to the following

    𝐰(n+1)=𝐰(n)+μ e^*(n)𝐱(n)-ρsgn(𝐰(n))

where ρ=μγ and sgn(x) is the component wise sign function defined as 

    sgn(x)=
    x/x ,   x ≠ 0 
    
    0,  otherwise

It is clear that smaller coefficients of the estimated vector are quickly drawn to zero making the vector sparse while larger coefficients remain mostly unaffected for small values of ρ. Thus the update rule converges to sparse vectors.


 §.§ ℓ_0-LMS

ℓ_0-LMS <cit.> takes a different approach to sparse system identification by trying to minimize the objective

    J_ℓ_0(n) = 1/2e(n)^2 +γ𝐰(n)_0

Of course simply doing a gradient descent on the objective directly is not possible and in general the problem is known to be NP-hard. Instead the ℓ_0 norm is approximated by 

    𝐰(n)_0 ≈∑_i=0^N-1(1-e^-βw_i(n))

The parameter β here controls the quality of the approximation of the ℓ_0 norm and as β tends to infinity the formula becomes exact. Taking the gradient on the modified objective leads to the following update rule

    𝐰(n+1)=𝐰(n)+μ e^*(n)𝐱(n)-ρsgn(𝐰(n))e^-β𝐰(n)

where the exponentiation and the sign is applied element-wise. The same observations as in the previous algorithms apply here also. The difference is that the attraction to zero is even weaker for the coefficients that have large magnitudes so we expect that the convergence should be faster in general.


§ NEW ONLINE ALGORITHMS



 §.§ Selective Zero Attracting LMS
 
In the previous two sub sections we saw two regularized objectives of the standard LMS objective. In this paper we will try to solve a constrained version of the LMS objective. We will try to minimize 

    J(n) = 1/2e(n)^2
 
but under the restriction that 

    𝐰(n)_0 ≤ s

where s, a positive integer less than N, is an upper bound on the sparsity of the vector under estimation that we know in advance. Let us define the operator H_s that outputs a vector having zeros in all coefficients except for the ones with the s largest absolute values that remain the same as in the input vector. For example if 𝐱_0=[2,-2,1,0]^T then H_2(𝐱_0)=[2,-2,0,0]^T. In case of ties we can take a conservative approach and allow all tying coefficients to be nonzero in the resulting vector so that H_1(𝐱_0)=[2,-2,0,0]^T. Thus |support(H_s(𝐱))| ≥ s and therefore it is not guaranteed that the output will always be s-sparse. The operator will give as output vectors that are not s-sparse when there are multiple coefficients in the vector that their absolute value is equal to the s largest absolute value in the vector. However, in most cases such ties will be nonexistent and the result will be an s-sparse vector. 
Given the definition of H_s one could easily see a connection with ℓ_0-LMS. Specifically we could relax the objective just like in the previous subsection. Here we will use however different β_i for each coefficient. Let us approximate the ℓ_0 norm as

    𝐰(n)_0 ≈∑_i=0^N-1(1-e^-β_iw_i(n))

Then if we want to make the estimate to converge to an s-sparse vector we can do the following: For the s algebraically largest coefficients we will use β_i=∞ whereas for all the others we will use β_i=0. This can be interpreted as the following penalty

    P_s(𝐱)_i=
    
    0,    i ∈support(H_s(x))
    sgn(x_i),   otherwise

which then leads to the following update rule

    𝐰(n+1)=𝐰(n)+μ e^*(n)𝐱(n)-ρ P_s(𝐰(n))

This is the same concept of the ℓ_1 penalization presented in <cit.> but applied only to the possibly superfluous coefficients given the a priori estimation of sparsity. We shall call this algorithm Selective Zero Attracting LMS.  Based on this fact we can prove a similar convergence result
 
Let us have a zero mean observation noise v(n) independent of 𝐱(n) and given that 𝐱(n) and 𝐰(n) are independent. Let us also assume that 𝔼[𝐱(n)𝐱^H(n)] is constant over time, invertible and equal to 𝐑_x. Then the algorithm described by (<ref>) converges in the mean sense provided that the condition of (<ref>) holds. The limiting vector satisfies the equation

    𝔼[𝐰(∞)]=𝐰-ρ/μ𝐑_x^-1𝔼[P_s(𝐰(∞))]


The proof of the theorem can be found in Appendix <ref> and is similar to the proof for the ZA-LMS. The interested reader can find an in depth analysis of a similar approximation scheme in <cit.>. The difference is in the choice of coefficients that get penalized and those who do not. In the update scheme presented we choose not to penalize the s largest coefficients. In <cit.> the coefficients that do not get penalized are those who are greater than a predefined threshold.
As we can see in Equation (<ref>), the expected value of the the estimation does not converge necessarily to 𝐰. In fact there might be a O(ρ) deviation per coefficient just like in the simple Zero Attracting LMS. However, if 𝐰 is an s sparse vector and the algorithm identifies the support correctly then the bias for the leading s coefficients should be eliminated as the penalty term will be zero for those coefficients, a property that the Zero Attracting LMS does not have. For the rest of the coefficients, unless the estimate for those does not converge exactly to 0 we will still incur the O(ρ) deviation, which should be negligible for small values of ρ.


 §.§ Hard Threshold LMS
 
The contribution of this paper is the study of the properties of the following update scheme

    𝐰(n+1)=H_s(𝐰(n)+μ e^*(n)𝐱(n))

It is easy to see the similarity of our proposed algorithm with the iterative hard thresholding algorithm studied in <cit.>, <cit.> and <cit.>. There, since the algorithm is developed in a batch setting where all the data are known in advance, the relation between the observations 𝐲 and the estimated vector 𝐰 is 𝐲=𝐀𝐰 where 𝐀 is M× N matrix with M<N; thus the problem is undetermined. The update of the iterative hard thresholding under similar assumptions for the sparsity of 𝐰 is

    𝐰(n+1)=H_s(𝐰(n)+μ𝐀^H𝐞(n))

where 𝐞(n)=𝐲-𝐀𝐰(n). It must be noted that the complexity of implementing such an operator is still linear in N as finding the s largest value in a vector does not require sorting it first.
As a result it is clear that the proposed algorithm is closely related to the special case of iterative hard thresholding having M=1. It is also clear that we cannot use the rigorous proofs found in <cit.>, <cit.> and <cit.> to show that the proposed algorithm also converges since for M=1 it is impossible to fulfill the strict properties needed. However, it is still possible to prove some interesting properties of the hard threshold operator. The main contribution of the operator is to let us focus our attention on the support of the estimated vector. If the algorithm does not provide a correct estimation of the support of the estimated vector then this could have a negative effect on the convergence of the algorithm. So one of the key properties that need to be studied is under which conditions is the estimation of the support using the hard threshold operator correct.
 
Let 𝐰=[w_0,w_1,…,w_N-1]^T ∈ℂ^N with 𝐰_0=s and 𝐰̂ be an approximation. Let q=min_w_i ≠ 0w_i. Then if 𝐰-𝐰̂_2^2 < q^2/2 the following will be true 

    support(H_s(𝐰̂))  = support(𝐰)


The proof of the theorem is quite involved and can be found in Appendix <ref>. The essence of the proof however is rather simple. In order to have the minimal error and still incorrectly specify the support of the vector, the error must be concentrated in two coefficients, one that belongs in support(𝐰) and one that does not. The one coefficient that belongs to the correct support must end up having a smaller  magnitude than the one that should not. Since the first coefficient has at least magnitude q in 𝐰 and the other coefficient must have magnitude 0, the minimal error is achieved when both have magnitude q/2 in  𝐰̂ which leads to the bound of the error that we have in the proof.
In order to understand the significance of the theorem we need to see some equivalent bounds having to do with the signal to error ratio that is needed so that the result in relation (<ref>) still holds. The true vector 𝐰 has s nonzero values  each with an absolute value of at least q. Thus 𝐰_2^2≥ sq^2 and hence we need

    SER=𝐰_2^2/𝐰-𝐰̂_2^2> sq^2/q^2/2=2s

Inequality (<ref>) is a necessary condition so that the required conditions of the theorem are true. Even if it is not sufficient it gives us the intuition that for small values of s it will be easier to come up with an estimate 𝐰̂ for which relation (<ref>) is true. On the other hand the conditions of Theorem <ref> are just sufficient for the relation (<ref>) so in practice relation (<ref>) could be true even with much lower signal to error ratios.
To further relax the conditions of our theorem we could allow the estimate to be less sparse. In order to do this we could use H_d instead of H_s with N>d>s>0 where N is the size of the estimated vector. What happens here is a trade off. On the one hand, the result now is less attractive since we have more nonzero coefficients than what is actually needed and that may lead to excessive estimation error that could possibly be avoided. On the other hand, the estimation error of the input to the threshold operator can be greater without risking of loosing an element of support(𝐰) after the application of the operator. The next theorem quantifies the gain in allowable estimation error.
 
Let 𝐰 be a vector in ℂ^N with 𝐰_0=s and 𝐰̂ be an approximation. Let q=min_w_i ≠ 0w_i and d=s+τ with d<N and τ>0 where s, τ, d are integers. Then if 𝐰-𝐰̂_2^2≤ q^2(1-1/τ+2) and 𝐰̂_0 ≥ d, the following will be true

    support(H_d(𝐰̂))  ⊇support(𝐰)
 

The proof of this theorem, found in the Appendix <ref>, is similar to the previous one. The difference in the result comes from the fact that τ+1 coefficients that are not part of support(𝐰) must have significant magnitudes in 𝐰̂ in order to miss a coefficient of support(𝐰).  The analogous inequality of relation (<ref>) for this theorem is

    SER≥s/(1-1/τ+2)

which is less strict as we have expected.
Given the last theorem one can overcome the need to have an initialization that is too close to the vector to be estimated. If we have an estimate that has an error 𝐰-𝐰̂_2^2 at most q^2, we can use the hard threshold operator to reduce its sparsity up to a degree that depends on the error without loosing an important coefficient and thus reducing the error in the process. Of course this is a worst case analysis and the conditions are sufficient but not necessary. Therefore in practice we should be able to to use the update rule of (<ref>) without waiting to converge so close to the solution.


 §.§ Estimating Sparsity

In some applications knowing an upper bound on sparsity, the parameter s in our algorithms, may be an acceptable assumption. For example, in echo cancellation one can assume that there will be a small batch of tens of coefficients that are non-zero. In spectrum estimation we can calibrate s based on prior knowledge about how many primary and secondary users of the spectrum are usually utilizing the spectrum. In general however, we would like our algorithm to adapt in different settings and therefore we need to be able to estimate the parameter s in an online fashion.
To achieve that we will assume that we have a knowledge of lower bound on q, the minimum magnitude of the non-zero coefficient in the estimated vector. One such lower bound could be the minimum magnitude required to consider the corresponding frequency occupied in the cognitive radio application. Let us call this value q^*. One naive way to estimate the sparsity could be to count the number of coefficients in the current estimate 𝐰(n) that have magnitude greater than q^* and use this as an estimate for the sparsity. 
Unfortunately, the current estimation may not be suitable to use for sparsity estimation when the underlying spectrum is changing. For example, let us assume that the number of non zero coefficients increases. To increase our estimation of s based on 𝐰(n) at least one coefficient's magnitude would need to go from 0 to above q^* in a single iteration. Waiting for multiple iterations does not help if hard thresholding is used to remove negligible coefficients. But for such a significant jump to happen in a single iteration one would need either a small q^* or a large μ both of which are undesirable as the first one reduces the accuracy of our sparsity estimate and the second one may make the estimation unstable. 
Instead we will try to approximate the error of our current estimate in order to construct a more suitable vector for the aforementioned procedure. The intuition behind this is that if we track an the error of our current estimate we can then use it to trigger increases in the parameter s when the error increases significantly. Let 𝐰 be once gain the true vector and 𝐰(n) our current estimate. We want to approximate 𝐰̅(n)=𝐰(n)-𝐰. From equation (<ref>) we can get by taking the expectation and assuming that the noise has zero mean that

    𝔼[e^*(n)𝐱(n)]=-𝔼[𝐱(n)𝐱^H(n)]𝔼[𝐰̅(n)]

𝐱(n) correspond to rows of Φ^*. Since they are chosen uniformly at random we know that 𝔼[𝐱(n)𝐱^H(n)]=Φ^*Φ^T=𝐈_N where 𝐈_N is the N× N identity matrix. This equality is based on the properties of the IDFT matrix. Therefore the equation becomes

    𝔼[e^*(n)𝐱(n)]=-𝔼[𝐰̅(n)]

Let 𝐞𝐫𝐫(n) be our approximation of 𝐰̅(n). Ideally, we could take a set of new or even past measurements and calculate e for them in every iteration to approximate the right hand side. This however would be wasteful and it would invalidate the online nature of the proposed algorithms. To avoid that we can reuse the estimations of the errors of the previous iterations. However, as our algorithm progresses, errors that were calculated many iterations ago in the past are not representative for our current estimate so they should be down-weighted compared to errors in recent iterations.  To overcome this we can take an exponentially weighted window average of the errors. Let λ∈ (0,1] be the forgetting factor of the window and 𝐛(n)=e^*(n)𝐱(n). Then we can write the following equations

    κ_n+1   = λκ_n +1 
    𝐞𝐫𝐫(n+1)    = ( 1 -1/κ_n+1) 𝐞𝐫𝐫(n) -1/κ_n+1𝐛(n)

where 𝐞𝐫𝐫(0) is all zeros and  κ_0 is zero as well. 
In the end we will get a 𝐰'(n)=𝐰(n)-𝐞𝐫𝐫(n) and we will compare each coefficient magnitude and compare it to the threshold q^*. The number of coefficients that pass this test is the estimate of s. Optionally we can use the coefficients that pass the test as the support of 𝐰(n+1) in the next iteration in order to reduce the computational effort. 
The advantage of using this process instead of operating directly on  𝐰(n) is that we can increase the importance of errors only for sparsity estimation and thus we avoid making our estimate unstable. In general we can even scale the importance of the error correction 

    𝐰'(n)=𝐰(n)-ξ𝐞𝐫𝐫(n)

where ξ is a positive real number. Holding q^* fixed we can increase ξ to make our sparsity estimation more sensitive to the error estimate. 


§ EXPERIMENTATION

 
In this section we will compare the performance of the various algorithms discussed previously. Let us first define the signals on which we will compare these algorithms. The signals of interest are going to be sums of sines affected by additive white Gaussian noise. Specifically the signals of interest here will have the form

    g(n)= ∑_i=1^k A_isin(2π f_i t(n)) + v(n)

where k will be the number of the signals added, f_i is the randomly chosen frequency of each sine, A_i is the amplitude of each sine wave and v(n) is the white zero mean noise. Therefore the spectrum of these signals will be sparse with s=2k non-zero coefficients. The sampling positions t(n) are spread uniformly in a time interval T and the corresponding sampling frequency is equal to the Nyquist frequency. This results in a vector of N samples per time interval T out of which we will sample M of those. Here we will assume for simplicity that A_i=1. 
The first thing we would like to show is that using sparsity aware techniques for spectrum estimation is a necessity when we are undersampling. We will therefore compare the spectrum estimations of the Hard Threshold LMS and the classical LMS. We will use the sum of k=10 sine waves of length N=1000 samples out of which we collect only M=300 samples corrupted by noise so that the SNR is equal to 20db. In order for the algorithms to converge we will make 10 iterations over the data. For the Hard Threshold LMS (HARD-LMS) algorithm we will use s=20 and we will refrain from thresholding during the first pass over the data. For both HARD-LMS and LMS we will use μ=1. The results can be seen in Figure <ref>. As we can clearly see the LMS algorithm does not converge to the true sparse vector that generated the measurements but simply to one of the many non-sparse solutions. In contrast HARD-LMS identified the support of the spectrum perfectly and the error is minimal compared to the one of LMS.
  
Moreover, we would like to show experimentally how the parameter s in our Hard Threshold LMS algorithm influences the speed of convergence as well as the steady state error. We set N to be equal to 1000 and M=200 leading to a one to 5 undersampling ratio. We set k to be equal to 10 and set the noise power so that the SNR of the observed samples is equal to 20db. We collect the M samples and repeat them K=100 times to test the convergence of the algorithms. We repeat the whole process of choosing the different frequencies f_i and the random M samples for R=200 times. The algorithms that are compared are the following: The Hard Threshold LMS algorithm for values s equal to 20 (HARD-20), 40 (HARD-40) and 80 (HARD-80) as well as the version of the Hard Threshold LMS with sparsity estimation (HARD-EST). For the sparsity estimation we use λ=0.99, q^* equal to one tenth of the magnitude of any of the non zero coefficients of the spectrum (all have equal magnitude in this case) and ξ=1. For all those algorithms we refrain from hard thresholding for the first 2M samples so that we get a good enough approximation. Additionally, for all algorithms μ is set equal to 1. We also include the standard LMS (LMS) as a baseline comparison.
The results we get from Figure <ref> are quite interesting. What is evaluated is the relative Mean Square Error (r-MSE). For every run of the algorithm and for every iteration we calculate

    r-MSE=𝐰-𝐰(n)^2/𝐰^2

and then we take the average performance in db. As we can see selecting s being exactly equal to the true sparsity is not always optimal. The algorithm for s=20 quickly converges to a suboptimal solution with high steady state error. This is because the algorithm has made a wrong estimation of the spectrum's support. In contrast allowing more non-zero coefficients allows the algorithm to include the true support of the spectrum as well as some superfluous coefficients. This allows both s=40 and s=80 to achieve much lower steady state errors. However, increasing the parameter s will tend to significantly decrease the speed of convergence of the algorithm. On the other hand the hard thresholding algorithm with sparsity estimation by making better estimates of the true spectrum and using a conservative magnitude threshold gradually decreases the sparsity estimate in order to smoothly and quickly converge. This aligns with what we proved in the previous section especially with Theorem <ref>. Of course the classical LMS algorithm had no hope of finding the true spectrum as the problem is undetermined and LMS gets stuck in a non sparse spectrum that could give the observed measurements. Since HARD-EST achieved the best performance compared to all other methods we will compare it with other online sparsity aware algorithms from the literature. 
     
Specifically we will also compare with Zero Attracting LMS (ZA-LMS) and Reweighted Zero Attracting LMS (RZA-LMS) from <cit.> as well as with ℓ_0-LMS from <cit.>. We will set the parameter ρ of ZA-LMS and RZA-LMS equal to 0.005 and ϵ=2.25. For the ℓ_0-LMS we will set  β=0.5 and κβ=0.005. We will also include the Selective Zero Attracting LMS (SZA-LMS) that we proposed in this paper using the true sparsity s=20 and ρ=0.005. Additionally, the proposed hard thesholding scheme with sparsity estimation can be combined with other more complicated update rules to further improve performance. So for this experiment we will combine it with the update rule of ℓ_0-LMS using the same parameters to show that we can get improved performance over each method alone. We shall call this algorithm HARD-ℓ_0 and its update rule will be

    𝐰(n+1)=H_s(𝐰(n)+μ e^*(n)𝐱(n)-ρsgn(𝐰(n))e^-β𝐰(n))
 
where s is estimated the same way as for HARD-EST. For HARD-EST we will refrain from hard thesholding for the first M samples and for HARD-ℓ_0 for the first 2M samples. The experimental settings are the same as in the previous experiment.

The results can be seen in Figure <ref> where we show again the r-MSE. Clearly, we can see that all algorithms manage to reach nearly the same level of r-MSE after some iterations with SZA-LMS and ℓ_0-LMS achieving a slightly smaller r-MSE than the other methods. So it makes sense to compare them in terms of speed of convergence. The fastest convergence belongs to SZA-LMS. SZA-LMS has ground truth knowledge of the sparsity of the vector just like the hard thresholding algorithms in the previous experiments but uses it slowly but steadily in order not to reduce coefficients of the true support to 0.  Then we have HARD-ℓ_0 which combines the quickly convergent update rule of ℓ_0-LMS with hard thresholding improving the convergence speed of an already fast algorithm like ℓ_0-LMS.  Then ℓ_0-LMS with an update rule that uses two parameters to tune the amount of zero attraction to use for each coefficient manages to converge faster than the simpler HARD-EST algorithm. HARD-EST then manages to converge faster than RZA-LMS. Finally the simple ZA-LMS algorithm fails to achieve a low steady state error. 
The third experiment that we will present has to with the robustness of the proposed procedures with varying degrees of undersampling. We evaluate the sparsity aware algorithms for M=100, 200, … 1000 samples where 1000 samples corresponds to the full measurement of the signal. In each setting we take 50 instantiations of the random sampling procedure. Then we calculate the steady state r-MSE after K=50 iterations over all the available measurements. The results are shown in Figure <ref>. We compare once again the same algorithms with the same parameters as in the previous experiment. We can clearly see that reducing the number of samples to 100 which corresponds to a 1 over 10 undersampling ratio is prohibitive for all algorithms except maybe SZA-LMS which has ground truth knowledge. However, once we get to 200 samples, which was the case in the previous experiment all algorithms improve their predictions considerably. Adding even more samples leads to better performance although with diminishing returns. One pattern that may seem counter-intuitive is that the hard thresholding algorithms, HARD-EST and HARD-ℓ_0, manage to outperform all other methods by a small margin after 200 samples which is in contrast to what we saw in the previous experiment. The reasoning behind this is that HARD-EST and HARD-ℓ_0 have no misalignment with the ground truth for the coefficients that are exactly zero. In contrast for the other methods these weights oscillate around zero  due to the zero attraction term for the same reasons we analyzed for the case of SZA based on Equation (<ref>). This difference in performance is quite small so it is only observable when HARD-EST and HARD-ℓ_0 have reached their optimal performance.


In the fourth experiment we are going to validate that HARD-EST is capable of tracking changing sparsity patterns in real time. In this experiment, the incoming signal will change over time. At first, the incoming signal consists of 10 sine waves just like before. The pattern of N=1000 samples is repeated 150 times. Then the incoming signal is augmented with another 10 sine waves of different frequency. Then the new pattern is repeated for 150 times. The incoming signal is split in non overlapping windows of N samples and we randomly sample M=200 measurements corrupted by additive noise in each window. The SNR is 20db.  
To help HARD-EST perform well in scenarios where the incoming signal is changing abruptly we are going to change the algorithms configuration. We are going to set λ=0.98, q^* equal to one hundredth of the magnitude of any of the non zero coefficients of the spectrum (again all the coefficients have equal magnitude) and ξ=20.  The smaller λ allows us to forget previous error estimates more quickly whereas the combination of the smaller q^* and the higher ξ allows us to adapt more quickly to changes in the sparsity of the spectrum. We can see the results in Figure <ref>. The algorithm converges very close to the true spectrum using the new samples it gets in every window. When the change in the spectrum happens the estimation is initially far away from the new spectrum. Then, the increased error estimates trigger the increase of the estimated sparsity from 20 towards 40 non zero coefficients. This allows the estimate to adapt to the less sparse spectrum and eventually converge to it. The r-MSE in the end is higher than before the spectrum change but this is to expected since the spectrum now is less sparse.
To understand the effect of the parameter ξ in Equation <ref> and the significance of our sparsity estimation procedure  we add an additional algorithm HARD-EST-SIMPLE which is the same as HARD-EST with the only difference being that for HARD-EST-SIMPLE we set ξ=0. In this setting the sparsity estimation is successful in the first half of the simulation yielding the same approximation error with HARD-EST. However, in the second half while HARD-EST manages to increase its estimate s to 40, HARD-EST-SIMPLE does not manage to adapt resulting in s being equal to 20 in the second half as well and in an r-MSE of -3 db. Therefore, it is clear that, when the underlying patterns of sparsity are changing, setting a positive value for ξ is crucial for the adaptation of the sparsity estimation.




§ CONCLUSION

In this paper we studied the problem of online sparse spectrum estimation for cognitive radio applications using sub-Nyquist sampling rates. To solve the problem, we analyzed the properties of two algorithms that try to minimize the squared error at each iteration while maintaining the ℓ_0 norm of the estimated vector under a predefined threshold s. Moreover, we analyzed the convergence properties of the Selective Zero Attracting LMS as well as the properties of the Hard Thresholding operator. Specifically, we proved that if our current estimation is close enough to the solution we can use the Hard Thresholding operator to reduce the error without risking to loose an important coefficient of the spectrum especially when we allow the operator to use more non-zero coefficients. Additionally, we proposed a way to estimate in an adaptive way the parameter s so that the estimation can gradually become sparser without misspecifying the support of the estimated spectrum. Further, in the experimentation section we analyzed the importance of the parameter s for the steady state error as well as the speed of convergence. Then we compared our results with other online sparsity aware algorithms in the literature. We also showed that the two proposed algorithms have robust performance even when the sampling rate is low and that they can produce even better estimates when the number of samples increases. Finally, we showed experimentally that the proposed sparsity estimation technique is robust to signal changes.
Of course spectrum estimation for cognitive radio applications is only one of the many possible applications of the proposed algorithms. Obviously an a priori estimation of the sparsity of the estimated vector may not be available in all use cases, even though we showed that this estimate must not be exact in order to actually take benefit. However, there are other use cases where the algorithms proposed here could make a difference. The standard LMS algorithm has been used in many adaptive machine learning tasks like neural network training and others as discussed in <cit.> so taking advantage of sparsity could be advantageous. For example, in the case of training a perceptron with an abundance of available features one could begin training with all the features but then proceed to use one of the proposed algorithms to impose feature selection through sparsity. By increasing the imposed sparsity one can then train several classifiers and then compare them using criteria like the Bayesian information criterion.



§ PROOF OF THEOREM <REF>
 

Let us define 𝐰̅(n) as the difference between the estimation 𝐰(n) and the true vector 𝐰. Subtracting 𝐰 from both sides of the equation (<ref>) gives

    𝐰̅(n+1)   =𝐰(n)+μ e^*(n)𝐱(n)-𝐰 -ρ P_s(𝐰(n))+v(n)
       =𝐰̅(n)+μ e^*(n) 𝐱(n) -ρ P_s(𝐰(n))+v(n)

After some calculations, which are the same as in the case of the classical LMS, we have that

    e^*(n)𝐱(n)=-𝐱(n)𝐱^H(n)𝐰̅(n)+v(n)𝐱(n)

Taking the mean under the independence assumptions made and given that the observation noise mean is zero will yield

    𝔼[e^*(n)𝐱(n)]=-𝐑_x𝔼[𝐰̅(n)]

Then from equation (<ref>) we obtain 

    𝔼[𝐰̅(n+1)]=(𝐈_N-μ𝐑_x)𝔼[𝐰̅(n)]-ρ𝔼[P_s(𝐰(n))]

where 𝐈_N is the N× N identity matrix. Given the bound in (<ref>) the algebraically largest eigenvalue of 𝐈_N-μ𝐑_x is less than one. Further the term induced by the penalty is bounded by the vectors -ρ1 and ρ1 where 1 is the vector of ℝ^N whose every element is one. Thus we can conclude that the  𝔼[𝐰̅(n)] converges and as a result so does  𝔼[𝐰(n)]. Therefore the algorithm provided by equation (<ref>) converges. The limiting vector cannot be found in a closed form but is guaranteed to be the solution of equation (<ref>).



§ PROOF OF THEOREM <REF>
 

The proof will be completed in three distinct cases.
(i) First, we assume that H_s(𝐰̂)_0<s which can be true only if 𝐰̂_0<s. We can easily see that, since  𝐰_0=s, there is at least one coefficient index i such that ŵ_i=0 and w_i≠ 0, which from the hypothesis also means that w_i≥ q. As a result we have that 
    𝐰-𝐰̂_2^2 ≥w_i-ŵ_i^2= w_i^2 ≥ q^2
 which contradicts the hypothesis; so this case is impossible.
(ii) Now we have that H_s(𝐰̂)_0 = s. Let us assume that relation (<ref>) does not hold. Then since the two sets have the same number of nonzero elements, it is clear that there is a coefficient index ℓ∈support(𝐰) but ℓ∉support(H_s(𝐰̂)) and a coefficient index k so that k ∈support(H_s(𝐰̂)) but k ∉support(𝐰). We directly know that w_k=0 and that w_ℓ≥ q. We can also deduce that ŵ_k > ŵ_ℓ since k belongs in support(H_s(𝐰̂)) but ℓ does not. Then, for the error norm we have 

    𝐰-𝐰̂_2^2 ≥w_k-ŵ_k^2+w_ℓ-ŵ_ℓ^2

Since w_k-ŵ_k^2 = ŵ_k^2 > ŵ_ℓ^2, it follows that 

    𝐰-𝐰̂_2^2 > 2ŵ_ℓ^2-w^*_ℓŵ_ℓ -w_ℓŵ^*_ℓ + w_ℓ^2
 
Therefore we can also write that

    𝐰-𝐰̂_2^2 >min_ŵ_ℓ∈ℂ2ŵ_ℓ^2-w^*_ℓŵ_ℓ -w_ℓŵ^*_ℓ + w_ℓ^2
 
The minimum value of the RHS is attained for ŵ_ℓ=w^*_ℓ/2 and equals w_ℓ^2/2; hence 

    𝐰-𝐰̂_2^2 >w_ℓ^2/2≥q^2/2
 
This once again contradicts the hypothesis and so relation (<ref>) is true in this case.
(iii) Finally, we assume that H_s(𝐰̂)_0 > s. This can happen only if there are ties for the s largest absolute values in 𝐰̂. Let us denote as B the set of tying coefficients, A= support(H_s(𝐰̂)) ∖ B and finally C=(support(H_s(𝐰̂))^c. It is evident that A≤ s-1. We shall prove that this case is impossible. There are two subcases:
(a) B ∩support(𝐰)= ∅. Since A≤ s-1 and w_0=s, support(𝐰) must have an element in common with C. Let us call that element ℓ. Let us also take an element k from B.  Then just like in the second case ŵ_k > ŵ_ℓ since k belongs in support(H_s(𝐰̂)) but ℓ does not. Following the rest of the steps in case (ii) we reach a contradiction.
(b) B ∩support(𝐰)≠∅. Let ℓ a common element of the two sets. Since H_s(𝐰̂)_0 >𝐰_0 there is an element k so that k ∈support(H_s(𝐰̂)) but k ∉support(𝐰). Since ℓ is one of the indexes tying for the last spot, we have ŵ_k≥ŵ_ℓ. Following the steps of case (ii) yields 𝐰-𝐰̂_2^2 ≥w_ℓ^2/2≥q^2/2 and therefore we get a contradiction.



§ PROOF OF THE THEOREM <REF>
 

Let us assume that relation (<ref>) does not hold. Just like in the proof of Theorem <ref> it is clear that there is a coefficient index so that ℓ∈support(𝐰) but ℓ∉support(H_d(𝐰̂)). This time however the set support(H_d(𝐰̂)) has at least d=s+τ elements but support(𝐰) has at most s-1 elements that could exist in  support(H_d(𝐰̂)). As a result we are sure that there are at least τ+1 indexes k_i so that k_i ∈support(H_s(𝐰̂)) but k_i ∉support(𝐰). Once again we know that w_k_i=0 and that w_ℓ≥ q and we can deduce that ŵ_k_i > ŵ_ℓ since k_i exists in support(H_d(𝐰̂)) but ℓ does not.
Like in the the proof of Theorem <ref> we can deduce about the error norm that 
    𝐰-𝐰̂_2^2 ≥∑_i=1^τ+1w_k_i-ŵ_k_i^2+w_ℓ-ŵ_ℓ^2

We bound the first term just like in the previous proof so that it becomes 
    ∑_i=1^τ+1w_k_i-ŵ_k_i^2 = ∑_i=1^τ+1ŵ_k_i^2≥ (τ+1) ŵ_ℓ^2
 Thus, we end up 
    𝐰-𝐰̂_2^2 > (τ+2)ŵ_ℓ^2-2w_ℓŵ_ℓ + w_ℓ^2

Taking the minimum on the right side with respect to ŵ_ℓ will lead once again to finding the minimum value of a quadratic function. The minimum is found for ŵ_ℓ= w_ℓ/τ+2 and equals to w_ℓ^2(1-1/τ+2); hence

    𝐰-𝐰̂_2^2 >w_ℓ^2(1-1/τ+2) ≥ q^2(1-1/τ+2)

which once again contradicts the hypothesis so the proof is completed.




§ ACKNOWLEDGMENT

We wish to thank the anonymous reviewers whose constructive comments helped us improve this paper.

IEEEtran

