DynamicX_BPDN_demo

 BPDN/LASSO update homotopy comparison
 Author: Muhammad Salman Asif
 Created: February 2009:

-------------------------------------------+
 Copyright (c) 2009.  Muhammad Salman Asif 
-------------------------------------------+

0001 % BPDN/LASSO update homotopy comparison
0002 % Author: Muhammad Salman Asif
0003 % Created: February 2009:
0004 %
0005 %-------------------------------------------+
0006 % Copyright (c) 2009.  Muhammad Salman Asif
0007 %-------------------------------------------+
0008 
0009 clear; clc;
0010 
0011 % % load fixed random states
0012 % load RandomStates
0013 % randn('state',s_randn);
0014 % rand('state',s_rand);
0015 
0016 % signal length
0017 N = 512;
0018 % number of spikes to put down
0019 T = 40;
0020 % number of observations to make
0021 M = 350;
0022 
0023 % random +/- 1 signal
0024 x = zeros(N,1);
0025 q = randperm(N);
0026 x(q(1:T)) = randn(T,1);%sign(randn(T,1)).*randint(T,1,256);
0027 
0028 % measurement matrix
0029 A = randn(M,N)/sqrt(M);
0030 
0031 % observations
0032 sigma = .01;
0033 e = randn(M,1)*sigma;
0034 y = A*x + e;
0035 
0036 % regularization parameter
0037 tau = 0.005*max(abs(A'*y)); %l1_ls
0038 if sigma>0
0039     tau = sigma * sqrt(log(N)*2); % BPDN
0040     %tau = max(abs(A'*(A*x-y))); % ideal ???
0041 end
0042 
0043 % Initial update using BPDN
0044 [xp, gamma_x, xp_iter, t1] = BPDN_homotopy_function(A, y, tau,4*N); %Lasso
0045 
0046 % Model the change in signal
0047 dx_model = 1; % 0 - Perturb only non-zero locations
0048               % 1 - Add some new elements (T_in) and remove some existing
0049               % elements (T_out)
0050 switch dx_model
0051     case  1
0052         % Add some elements (T_in) or remove some existing elements (T_out)
0053         dx  = zeros(N,1);
0054         dx(q(1:T)) = randn*.1;
0055         T_in = round(T/10);
0056         T_out = round(T/10);
0057         qin = randperm(N);
0058 
0059         dx(qin(1:T_in)) = (randn(T_in,1));
0060 
0061         qout = randperm(T);
0062         oldx_indices = q(1:T);
0063         dx(oldx_indices(qout(1:T_out))) = -x(oldx_indices(qout(1:T_out)));
0064 
0065     case  2;
0066         % Perturb the current locations by Gaussian noise
0067         oldx_indices = q(1:T);
0068         dx = zeros(N,1);
0069         dx(oldx_indices) = randn(T,1)*.2;
0070     otherwise
0071         disp('Nooooo');
0072 end
0073 
0074 xt = x+dx;
0075 yt = y+A*dx;
0076 
0077 % Gram matrix update
0078 AtAgx = A(:,gamma_x)'*A(:,gamma_x);
0079 iAtAgx = inv(AtAgx);
0080 pk = A'*(A*xp-y);
0081 
0082 [xp_h, gamma_xh, xp_h_iter, th] =  DynamicX_BPDN_function(A, AtAgx, iAtAgx, y, yt, xp, gamma_x, pk, tau, 4*N);
0083 
0084 [xp2, gamma_x2, xp2_iter, t2] = BPDN_homotopy_function(A, yt, tau, 4*N); %Lasso
0085 
0086 % [old_homotopy new_homotopy update_homotopy]
0087 Cputime_comparison = [t1 t2 th]
0088 Iteration_comparison = [xp_iter xp2_iter xp_h_iter]
0089 
0090 % figure(2); plot(xp_h-xp2); shg