0
$\begingroup$

The model is expressed as, $$y(n) = \sum_{i=0}^{p-1} r(i) x(n-i) + v(n) \tag{1}$$ where $\mathbf{r} = [r_1,r_2,\ldots,r_p]^T$ is the sparse channel coefficients of length $p$, $\mathbf{x} = [x_1,x_2,...,x_{n-p}]^T$ is the one dimensional input and v is additive White Gaussian noise of zero mean. y = filter(.) command is used to model the above equation and thus creating an FIR filter or a moving average (MA) model. The order of the MA model is p=3.

So, $y = [y(1),y(2),....,y(100)]$ is a vector of 100 elements. I am generating noise of variance 0.1 I want to estimate the sparse channel coefficients using LASSO. As there are p channel coefficients, I should get p estimates.

According to the equation of LASSO, ||rx - y||_2^2 + lambda * ||r||_1 I am estimating the sparse coefficients, r. As the true coefficient array contains p elements, I should get p estimated elements. I am not quite sure if this is the way to do. I have not found any example on LASSO being applied to univariate time series model such as ARMA. I don't know how to estimate the sparse coefficients using the appropriate algorithm and need help.

The first part of the Equation : $||rx - y||_2^2$ is a least squares formulation which I can solve using Least Squares Approach. In order to implement LS, I have to arrange the input in terms of regressors. However, if the coefficients, $\mathbf{r}$ are sparse then I should use LASSO approach. I have tried using Matlab's LASSO function. For LASSO, I rearranged the input data $x$ in terms of regressors, but I don't know if this the correct approach.

I need help. Is there an approach to include the sparsity term in the LS?

Please find below the code for LASSO using Matlab function. As a toy example I am just assuming model order to be of lag 3 but I know that LASSO can be applied efficiently to a large model. I can test for larger order MA model having lag > 3.


% Code for LASSO estimation technique for 
%MA system, L = 3 is the order,  

%Generate input
 x = -5:.1:5;

r = [1    0.0   0.0];% L elements of the channel coefficients     
%Data preparation into regressors    
X1 = [ ones(length(x),1) x' x']; %first column treated as all ones since    x_1=1

y = filter(r,1,x); % Generate the MA model
[r_hat_lasso, FitInfo] = lasso(X1, y, 'alpha', 1, 'Lambda', 1, 'Standardize', 1);

OUTPUT :

The estimates returned are r_hat_lasso = 0, 0.657002829714982, 0

Question : This differs very much from the actual r. Is this wrong?

$\endgroup$
2
$\begingroup$

I am not entirely sure what matlab's LASSO routine does so I started with Ordinary Least Squares (OLS) and worked backwards. From an OLS perspective X1 as you have it won't work. You've got a regressor that is as all ones, but your parameter inputs for you example data (r) doesn't contain an offset. Essentially your model doesn't fit you're data well, which is why you get a weird answer.

Second, you have repeated regressors which is going to confuse things even more. I've rewritten your code to shift your regressors to line up the three inputs that would result from MA using p = 3.

You can perform regular old OLS on the data X1 and y to get the exact input r. If I use matlab's LASSO function I can get the proper components of r but they are off by a scalar.

    % Code for LASSO estimation technique for 
    %MA system, L = 3 is the order,  

    %Generate input
    x = -5:.1:5;
    xm = circshift(x,[0,1]);
    xm2 = circshift(x,[0,2]);

    r = 0.5*[1.0    1.0   0.0]'; % L elements of the channel coefficients     

    %Data preparation into regressors    
    X1 = [ x' xm' xm2'];

    y = X1*r;
    y1 = filter(r,1,x); % Generate the MA model, so I can plot and check.

    r_hat = X1'*X1\(X1'*y);

    [r_hat_lasso, FitInfo] = lasso(X1, y, 'alpha', 1, 'Lambda', 1, 'Standardize', 1);

Output:

    r_hat =
    0.5000
    0.5000
    0.0000

    r_hat_lasso =
    0.3235
    0.3231
         0
$\endgroup$
  • $\begingroup$ Thank you for your reply. But I think your code has a typo resulting in the error Undefined function or variable 'xp'. in the line X1 = [x' xm' xp']. The corrected should be in the line X1 = [x' xm' xm2'] Can you please explain how to use the cirshift() command?If the order is say p=5, which means that r contains 5 elements, then would I be doing ` xm = circshift(x,[0,1]); xm2 = circshift(x,[0,2]); xm3 = circshift(x,[0,3]); ` xm4 = circshift(x,[0,4]);` X1 = [x' xm' xm2' xm3' xm4'] ? $\endgroup$ – SKM May 1 '17 at 23:12
  • $\begingroup$ Lastly, is it possible to see what the value of \lambda is selected by Matlab if I don't specify the value in the function signature? It seems LS is performing better than LASSO when there are less number of non-zero coefficients. The reason for developing LASSO was to handle such sparse casess. However, it seems here that LS is better. This is quite confusing. Do you have thoughts on why LS estimation performance is better for estimating zero coefficients as well? $\endgroup$ – SKM May 1 '17 at 23:30
  • $\begingroup$ @SKM Yes, that is how you use circshift. The matlab help can probably tell you whether to use xm3 = circshift(x,[0,3]); or xm3 = circshift(x,[0,-3]); to go forward or backward. I can never remember. $\endgroup$ – Chad Sexington May 2 '17 at 1:43
  • $\begingroup$ @SKM The values for Lamba are in the FitInfo. If you add in noise OLS doesn't do any better than lasso. They both pick out the zeros but they get the coefficients wrong -- as one would expect with noise. This exchange talks about why one would want to use lasso over OLS: stats.stackexchange.com/questions/82466/… $\endgroup$ – Chad Sexington May 2 '17 at 1:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.