DynamicSeq_BPDN_demo

 BPDN update with sequential measurements homotopy comparison
 Author: Salman Asif
 Created: February 2009:

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

0001 % BPDN update with sequential measurements homotopy comparison
0002 % Author: 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 = 256;
0018 % number of spikes to put down
0019 T = 30;
0020 % number of observations to make
0021 K = 100;
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 G = randn(K,N)/sqrt(K);
0030 % G = orth(G')';
0031 A = G; %(1:end-1,:);
0032 % b = G(end,:);
0033 
0034 % observations
0035 sigma = .01;
0036 e = randn(K,1)*sigma;
0037 y = A*x + e;
0038 
0039 % regularization parameter
0040 tau = 0.005*max(abs(A'*y)); %l1_ls
0041 if sigma>0
0042     tau = sigma * sqrt(log(N)*2); % BPDN
0043     %tau = max(abs(A'*(A*x-y))); % ideal ???
0044 end
0045 
0046 % Initial update using Lasso
0047 [xp, gamma_x, xp_iter, t1] = BPDN_homotopy_function(A, y, tau,4*N); %Lasso
0048 
0049 chk_e0 = 1;  % 0 or 1: This selects if we want to take U = (A'A+b'b) (1) or (A'A+e0*b'b) (0)...!
0050 % In simulation U = (A'A+b'b) works better
0051 
0052 b = randn(1,N);
0053 ew = randn(1,1)*sigma;
0054 w = b*x+ew;
0055 
0056 % Gram matrix update
0057 AtAgx = A(:,gamma_x)'*A(:,gamma_x)+chk_e0*b(:,gamma_x)'*b(:,gamma_x);
0058 iAtAgx = inv(AtAgx);
0059 pk = A'*(A*xp-y);
0060 
0061 [xp_h, gamma_xh, xp_h_iter, th] =  DynamicSeq_BPDN_function(A, b, AtAgx, iAtAgx, y, w, xp, gamma_x, pk, tau, chk_e0, 4*N);
0062 
0063 
0064 y2 = [y; w];
0065 G = [A; b];
0066 [xp2, gamma_x2, xp2_iter, t2] = BPDN_homotopy_function(G, y2, tau, 4*N); %Lasso
0067 
0068 % [old_homotopy new_homotopy update_homotopy]
0069 disp(' ');
0070 disp('Results for ');
0071 disp('old_homotopy , new_homotopy , update_homotopy')
0072 cputime_comparison = [t1 t2 th]
0073 Iter_comparison = [xp_iter xp2_iter xp_h_iter]
0074 
0075 % figure(2); plot(xp_h-xp2); shg