script_simulation_wtBPDN

 partial script for running different solvers

0001 % partial script for running different solvers
0002 
0003 %% Main loop
0004 SIM_stack = cell(maxsim,1);
0005 SIM_memory = cell(maxsim,1);
0006 
0007 for sim = 1:maxsim
0008     %% Generate a random signal
0009     in = []; in.type = sType; in.T = T; in.randgen = 1; in.take_fwt = 0;
0010     switch sType
0011         case 'blocks'
0012             in.take_fwt = 1;
0013             in.wType = 'haar';
0014         case 'heavisine'
0015             in.take_fwt = 1;
0016             in.wType = 'daub4';
0017         case 'pcwpoly';
0018             in.take_fwt = 1;
0019             in.wType = 'daub8';
0020     end
0021     x = genSignal(N,in);
0022     % gamma_orig = find(abs(x)>0);
0023     [val ind] = sort(abs(x),'descend');
0024     ind_pos = ind(val>0);
0025     gamma_orig = ind_pos(1:min(length(ind_pos),M-1));
0026     
0027     % measurement matrix
0028     in = []; in.type = mType;
0029     A = genAmat(M,N,in);
0030     
0031     % measurements
0032     sigma = sqrt(norm(A*x)^2/10^(SNR/10)/M);
0033     %sigma = .05;
0034     e = randn(M,1)*sigma;
0035     y = A*x+e;
0036     
0037     % % orthogonalize rows of A even if A itself is not orthogonal
0038     %     [Q, R] = qr(A',0);
0039     %     A = Q'; y = R' \ y;
0040     
0041     % parameter selection
0042     if lambda > 0
0043         tau = lambda*max(abs(A'*y));
0044     else
0045         % tau = sigma*sqrt(log(N));
0046         tau = max(1e-4*max(abs(A'*y)),sigma*sqrt(log(N)));        
0047     end
0048     % tau = max(lambda*max(abs(A'*y)),sigma*sqrt(log(N)));
0049     
0050     maxiter = 2*N;
0051     
0052     %% Adaptive support & weight selection
0053     in = [];
0054     in.A = A; in.y = y; in.x = x;
0055     in.x_init = zeros(N,1); in.max_rwt = rwt_adp;
0056     in.tau = tau;
0057     in.rwt_mode = rwt_mode;
0058     in.delx_mode = delx_mode;
0059     in.debias = 0;
0060     in.verbose = 0;
0061     if in.verbose
0062         fprintf('sType-%s, mType-%s, (N,M,T)-%d,%d,%d ',sType,mType,N,M,T);    
0063     end
0064     out = script_rwtBPDN_adaptive(in);
0065     xh_adp = out.x_out;
0066     iter_adp = out.iter;
0067     gamma_adp = out.gamma;
0068     W_new = out.W_new;
0069     xh_adp_init = out.x_init;
0070     err_adp = out.err;
0071     time_adp = out.time;
0072     supp_diff_adp = out.supp_diff;
0073     
0074     %% Iterative reweighted BPDN (IRW-H)
0075     in = [];
0076     in.A = A; in.y = y; in.x = x;
0077     in.x_init = zeros(N,1); in.max_rwt = rwt;
0078     in.tau = tau;
0079     in.rwt_mode = rwt_mode;
0080     in.delx_mode = delx_mode;
0081     out = script_rwtBPDN_iterative(in);
0082     xh_rwt = out.x_out;
0083     xh_rwt_init = out.x_init;
0084     iter_rwt = out.iter;
0085     gamma_rwt = out.gamma;
0086     err_rwt = out.err;
0087     time_rwt = out.time;
0088     supp_diff_rwt = out.supp_diff;
0089     
0090     %% Oracle rwt BPDN
0091     alpha = 5; beta = 10;
0092     W_new = ones(N,1); W_new(gamma_orig) = 1/alpha./(beta*abs(x(gamma_orig)));
0093     % epsilon = 0.1; W_new = 1/alpha./(beta*(abs(x))+epsilon);
0094     
0095     % To check the accuracy of rwtBPDN_adaptive set
0096     % W_new = out.W_new/tau;
0097     
0098     in = [];
0099     in.tau = tau;
0100     in.maxiter = maxiter;
0101     in.x_orig = x;
0102     in.record = 1;
0103     in.delx_mode = 'mil';
0104     %     in.Te = rwt_itr;
0105     AW = A*diag(1./W_new);
0106     tic
0107     out_new = BPDN_homotopy_function(AW, y, in); %BPDN
0108     % time_orac = out_new.time;
0109     time_orac = toc;
0110     xh_orac = out_new.x_out.*(1./W_new);
0111     gamma_orac = out_new.gamma;
0112     iter_orac = out_new.iter;
0113     
0114     %% oracle LS
0115     x_LS = zeros(N,1);
0116     x_LS(gamma_orig) = A(:,gamma_orig)\y;
0117     
0118     %% OMP
0119     in = [];
0120     in.Te = round(1.2*T);
0121     out = OMP_function(y,A,in);
0122     x_omp = out.x_out;
0123     iter_omp = out.iter;
0124     gamma_omp = out.gamma;
0125     
0126     %% FPC_AS
0127     %     opts = [];
0128     %     opts.gtol = 1e-6;
0129     %     opts.xtol = 1e-6;
0130     %     opts.fullMu = false;
0131     %     opts.x0 = zeros(N,1);
0132     %     opts.eta = 4;
0133     %     opts.mxitr = N;
0134     %     % opts.PrintOutput = 0;
0135     %     % opts.f_value_tol = tolA;
0136     %     % opts.PrintOptions = 0;
0137     %     % opts.record = -1;
0138     %     % opts.scale_A = 1;
0139     %     W_new = ones(N,1);
0140     %     AW = A;
0141     %     x_fpc = [];
0142     %     iter_fpc = [];
0143     %     err_fpc = [];
0144     %     supp_diff_fpc = [];
0145     %     for rwt_itr = 1:rwt+1
0146     %         % [x_FPC, out_FPC] = FPC_AS(N,AW,y,tau,[],opts); % call FPC_AS
0147     %         out_FPC = fpc(N,AW,y,1/tau,[],opts);
0148     %         out_FPC.nProdA = out_FPC.itr;
0149     %         out_FPC.nProdAt = out_FPC.itr;
0150     %         x_FPC = out_FPC.x;
0151     %
0152     %         x_fpc = x_FPC./W_new;
0153     %
0154     %         [alpha beta epsilon] = weight_param(rwt_mode,rwt_itr,x_fpc,M);
0155     %         W_new = 1/alpha./(beta*abs(x_fpc)+epsilon);
0156     %         opts.x0 = x_fpc.*W_new;
0157     %         AW = A*diag(1./W_new);
0158     %         iter_fpc = [iter_fpc (out_FPC.nProdA+out_FPC.nProdAt)/2];
0159     %         err_fpc = [err_fpc norm(x-x_fpc)/norm(x)];
0160     %         supp_diff_fpc = [supp_diff_fpc length(setxor(gamma_orig,find(x_fpc)))];
0161     %     end
0162     
0163     %% SPGL1
0164     % call SPGL1
0165     iter_spgl1 = []; err_spgl1 = []; supp_diff_spgl1 = []; time_spgl1 = [];
0166     if ~exist('spgSetParms','file'); error('Solver SPGL1 is not found.'); end
0167     
0168     W_new = ones(N,1);
0169     % delta = norm(e);
0170     delta = sqrt(M)*sigma;
0171     x_spgl1 = [];
0172     
0173     for rwt_itr = 1:rwt+1
0174         spg_opts = spgSetParms('verbosity',0,'weights',W_new,'iterations',maxiter);
0175         [x_spgl1,r,g,info] = spgl1(A,y,0,delta,x_spgl1,spg_opts);
0176         
0177         rerr = norm(x_spgl1-x)/norm(x);
0178         iter_spgl1 = [iter_spgl1  (info.nProdA+info.nProdAt)/2];
0179         err_spgl1 = [err_spgl1 rerr];
0180         supp_diff_spgl1 = [supp_diff_spgl1 length(setxor(gamma_orig,find(x_spgl1)))];
0181         time_spgl1 = [time_spgl1 info.timeTotal];
0182         [alpha beta epsilon] = weight_param(rwt_mode,rwt_itr,x_spgl1,M);
0183         W_new = 1/alpha./(beta*abs(x_spgl1)+epsilon);
0184     end
0185     
0186     %% YALL1
0187     % set options
0188     digit = 6; if sigma > 0; digit = 4; end
0189     opts = [];
0190     opts.tol = 10^(-digit);
0191     W_new = ones(N,1);
0192     opts.weights = W_new;
0193     opts.print = 0;
0194     opts.maxit = maxiter/2;
0195     x_yall1 = [];
0196     
0197     iter_yall1 = [];
0198     err_yall1 = [];
0199     time_yall1 = [];
0200     supp_diff_yall1 = [];
0201     for rwt_itr = 1:rwt+1
0202         opts.nu = 0; opts.rho = tau; 
0203         tic;
0204         [x_yall1,Out] = yall1(A,y,opts);
0205         % time_yall1 = [time_yall1 Out.cputime];
0206         time_yall1 = [time_yall1 toc];
0207         
0208         iter_yall1 = [iter_yall1 (Out.cntA+Out.cntAt)/2];
0209         rerr = norm(x-x_yall1)/norm(x);
0210         err_yall1 = [err_yall1 rerr];
0211         
0212         supp_diff_yall1 = [supp_diff_yall1 length(setxor(gamma_orig,find(x_yall1)))];
0213         
0214         [alpha beta epsilon] = weight_param(rwt_mode,rwt_itr,x_yall1,M);
0215         W_new = 1/alpha./(beta*abs(x_yall1)+epsilon);
0216         opts.x0 = x_yall1;
0217         opts.weights = W_new;
0218     end
0219     
0220     %% SpaRSA
0221     W_new = 1;
0222     % AW = A;
0223     x_sparsa = 0;
0224     iter_sparsa = [];
0225     err_sparsa = [];
0226     time_sparsa = [];
0227     supp_diff_sparsa = [];
0228     for rwt_itr = 1:rwt+1
0229         psi_function = @(x,tau) soft(x,tau*W_new);
0230         phi_function = @(x) sum(abs(W_new.*x));
0231         tic;
0232         [x_SpaRSA,x_debias_SpaRSA_m,obj_SpaRSA_m_cont,...
0233             times_SpaRSA_m_cont,debias_start_SpaRSA_m,mse_SpaRSA_m,taus_SpaRSA_m, numA, numAt]= ...
0234             SpaRSA_adpW(y,A,tau,...
0235             'Monotone',0,...
0236             'adp_wt',0,...
0237             'W_new',W_new,...
0238             'Debias',0,...
0239             'Initialization',x_sparsa,...
0240             'StopCriterion',2,...            
0241             'ToleranceA',1e-4,...
0242             'psi',psi_function,...
0243             'phi',phi_function,...
0244             'Safeguard',1,...
0245             'MaxiterA',maxiter,...
0246             'Verbose',0,...
0247             'True_x',x,...
0248             'Continuation',1,...
0249             'Continuationsteps',-1);
0250         x_sparsa = x_SpaRSA;
0251         
0252         time_sparsa = [time_sparsa toc];
0253         % time_sparsa = [time_sparsa times_SpaRSA_m_cont(end)];
0254 
0255         iter_sparsa = [iter_sparsa (numA+numAt)/2];
0256         rerr = norm(x-x_sparsa)/norm(x);
0257         err_sparsa = [err_sparsa rerr];
0258         
0259         
0260         if rerr > 100
0261             fprintf('rerr big... iter # %d in %s\n',sim,str0);
0262         end
0263         
0264         supp_diff_sparsa = [supp_diff_sparsa length(setxor(gamma_orig,find(x_sparsa)))];
0265         
0266         [alpha beta epsilon] = weight_param(rwt_mode,rwt_itr,x_sparsa,M);
0267         W_new = 1/alpha./(beta*abs(x_sparsa)+epsilon);
0268     end
0269     
0270     %% NESTA
0271     %     opts = [];
0272     %     opts.Verbose = 0;
0273     %     opts.tolvar = 1e-8;
0274     %     delta = 0;
0275     %     muf = 1e-8;
0276     %     W_new = ones(N,1);
0277     %
0278     %     iter_nesta = [];
0279     %     err_nesta = [];
0280     %     supp_diff_nesta = [];
0281     %     x_nesta = [];
0282     %
0283     %     for rwt_itr = 1:rwt+1
0284     %         opts.U = spdiags(W_new,0,N,N);
0285     %         opts.Ut = opts.U;
0286     %         opts.normU = max(W_new);
0287     %         opts.xplug = x_nesta;  % use old solution as starting value
0288     %
0289     %         % constrained
0290     %         [x_nesta,niter,resid,outData,optsOut] = NESTA(A,[],y,muf,delta,opts);
0291     %         % unconstrained
0292     %         La = norm( A*A' );
0293     %         % [x_nesta,niter,resid,outData,optsOut] = NESTA_UP(A,[],y,tau,La,muf,opts);
0294     %
0295     %         iter_nesta = [iter_nesta niter*3];
0296     %         rerr = norm(x-x_nesta)/norm(x);
0297     %         err_nesta = [err_nesta rerr];
0298     %         supp_diff_nesta = [supp_diff_nesta length(setxor(gamma_orig,find(x_nesta)))];
0299     %
0300     %         [alpha beta epsilon] = weight_param(rwt_mode,rwt_itr,x_nesta,M);
0301     %         W_new = 1/alpha./(beta*abs(x_nesta)+epsilon);
0302     %     end
0303     
0304     %% SALSA
0305     %     W_new = ones(N,1);
0306     %     x_salsa = 0;
0307     %     mu = 1/tau;
0308     %     invLS = @(x) inv(A'*A+mu*eye(N))*x;
0309     %
0310     %     iter_salsa = [];
0311     %     err_salsa = [];
0312     %     supp_diff_salsa = [];
0313     %     for rwt_itr = 1:rwt+1
0314     %         psi_function = @(x,tau) soft(x,tau*W_new);
0315     %         phi_function = @(x) sum(abs(W_new.*x));
0316     %         [x_salsa, numA, numAt, objective, distance,  times, mses] = ...
0317     %             SALSA_v2(y, A, tau,...
0318     %             'Mu', 1e-3, ...
0319     %             'True_x', x, ...
0320     %             'psi',psi_function,...
0321     %             'phi',phi_function,...
0322     %             'INITIALIZATION', x_salsa, ...
0323     %             'ToleranceA', 1e-6,...
0324     %             'MAXITERA', maxiter, ...
0325     %             'VERBOSE', 0);
0326     %
0327     %         iter_salsa = [iter_salsa (numA+numAt)/2];
0328     %         rerr = norm(x-x_salsa)/norm(x)
0329     %         err_salsa = [err_salsa rerr];
0330     %         supp_diff_salsa = [supp_diff_salsa length(setxor(gamma_orig,find(x_salsa)))];
0331     %
0332     %         [alpha beta epsilon] = weight_param(rwt_mode,rwt_itr,x_salsa,M);
0333     %         W_new = 1/alpha./(beta*abs(x_salsa)+epsilon);
0334     %     end
0335     
0336     %% Save results...
0337     
0338     % fprintf('sim iter. %d: \n',sim);
0339     % SIM_stack = [SIM_stack; sim, rwt_itr, iter_adp, iter_rwt_ALL, norm(x-xh_adp), norm(x-xh_rwt),length(setxor(gamma_adp,gamma_rwt)),length(setdiff(gamma_adp,gamma_orig)),length(setdiff(gamma_rwt,gamma_orig))];
0340     SIM_stack{sim} = [sim, tau, norm(x-x_LS)/norm(x), norm(x-xh_orac)/norm(x), iter_orac, ...
0341         norm(x-xh_adp)/norm(x), sum(iter_adp,2), sum(time_adp,2), ...
0342         norm(x-xh_rwt)/norm(x), sum(iter_rwt,2), sum(time_rwt,2), ...
0343         norm(x-x_yall1)/norm(x), sum(iter_yall1,2), sum(time_yall1,2), ...
0344         norm(x-x_sparsa)/norm(x), sum(iter_sparsa,2), sum(time_sparsa,2), ...
0345         norm(x-x_spgl1)/norm(x), sum(iter_spgl1,2), sum(time_spgl1,2), ...
0346         norm(x-x_omp)/norm(x), iter_omp];
0347     % norm(x-x_fpc)/norm(x), sum(iter_fpc,2), ...
0348     % sum(iter_nesta,2), norm(x-x_nesta)/norm(x), ...
0349     % sum(iter_salsa,2), norm(x-x_salsa)/norm(x), ...
0350     
0351     % adaptive rwt
0352     exp = 1;
0353     SIM_memory{sim}{exp,1} = 'adp';
0354     SIM_memory{sim}{exp,2} = iter_adp;
0355     SIM_memory{sim}{exp,3} = err_adp;
0356     SIM_memory{sim}{exp,4} = supp_diff_adp;
0357     SIM_memory{sim}{exp,5} = time_adp;
0358     % iterative rwt
0359     exp = exp+1;
0360     SIM_memory{sim}{exp,1} = 'rwt';
0361     SIM_memory{sim}{exp,2} = iter_rwt;
0362     SIM_memory{sim}{exp,3} = err_rwt;
0363     SIM_memory{sim}{exp,4} = supp_diff_rwt;
0364     SIM_memory{sim}{exp,5} = time_rwt;
0365     % yall1
0366     exp = exp+1;
0367     SIM_memory{sim}{exp,1} = 'yall1';
0368     SIM_memory{sim}{exp,2} = iter_yall1;
0369     SIM_memory{sim}{exp,3} = err_yall1;
0370     SIM_memory{sim}{exp,4} = supp_diff_yall1;
0371     SIM_memory{sim}{exp,5} = time_yall1;
0372     % sparsa
0373     exp = exp+1;
0374     SIM_memory{sim}{exp,1} = 'sparsa';
0375     SIM_memory{sim}{exp,2} = iter_sparsa;
0376     SIM_memory{sim}{exp,3} = err_sparsa;
0377     SIM_memory{sim}{exp,4} = supp_diff_sparsa;
0378     SIM_memory{sim}{exp,5} = time_sparsa;
0379     % spgl1
0380     exp = exp+1;
0381     SIM_memory{sim}{exp,1} = 'spgl1';
0382     SIM_memory{sim}{exp,2} = iter_spgl1;
0383     SIM_memory{sim}{exp,3} = err_spgl1;
0384     SIM_memory{sim}{exp,4} = supp_diff_spgl1;
0385     SIM_memory{sim}{exp,5} = time_spgl1;
0386     %     % fpc
0387     %     exp = exp+1;
0388     %     SIM_memory{sim}{exp,1} = 'fpc';
0389     %     SIM_memory{sim}{exp,2} = iter_fpc;
0390     %     SIM_memory{sim}{exp,3} = err_fpc;
0391     %     SIM_memory{sim}{exp,4} = supp_diff_fpc;
0392     %     % nesta
0393     %     exp = exp+1;
0394     %     SIM_memory{sim}{6,1} = 'nesta';
0395     %     SIM_memory{sim}{6,2} = iter_nesta;
0396     %     SIM_memory{sim}{6,3} = err_nesta;
0397     %     SIM_memory{sim}{6,4} = supp_diff_nesta;
0398     %     % salsa
0399     %     exp = exp+1;
0400     %     SIM_memory{sim}{6,1} = 'salsa';
0401     %     SIM_memory{sim}{6,2} = iter_salsa;
0402     %     SIM_memory{sim}{6,3} = err_salsa;
0403     %     SIM_memory{sim}{6,4} = supp_diff_salsa;
0404     
0405     % mS = mean(cell2mat(SIM_stack),1);
0406     % fprintf('maxsim %d. oracLS=%3.4g. (err,iter,time): oracWT-%3.4g,%3.4g; adp-%3.4g,%3.4g,%3.4g; rwt-%3.4g,%3.4g,%3.4g; yall1-%3.4g,%3.4g,%3.4g; sparsa-%3.4g,%3.4g,%3.4g; omp-%3.4g,%3.4g; fpc-%3.4g,%3.4g; spgl1-%3.4g,%3.4g. \n', maxsim, mS(2:end));
0407 end
0408 
0409 mS = mean(cell2mat(SIM_stack),1);
0410 str1 = 'maxsim %d. tau = %3.4g, oracLS=%3.4g. (err,iter,time): oracWT-%3.4g,%3.4g; adp-%3.4g,%3.4g,%3.4g; rwt-%3.4g,%3.4g,%3.4g; yall1-%3.4g,%3.4g,%3.4g; sparsa-%3.4g,%3.4g,%3.4g; spgl1-%3.4g,%3.4g,%3.4g; omp-%3.4g,%3.4g.';
0411 str2 = sprintf([str1,' \n'], maxsim, mS(2:end));
0412 disp(str2);
0413