Comparison of various solvers for iterative reweighting and adaptive reweighting
Solves the following basis pursuit denoising (BPDN) problem
min_x \Sum \w_i |x_i| + 1/2*||y-Ax||_2^2
while updating the weights w_i
Written by: Salman Asif, Georgia Tech
Email: sasif@gatech.edu
Created: June 16, 2011
Reference:
"Fast and accurate algorithms for re-weighted L1 norm minimization," by
M. Salman Asif and Justin Romberg
To reproduce experiments in the paper, use this scrip with
rseed = 2012;
rand('state',rseed);
randn('state',rseed);
before running each simulation in
script_simulation_wtBPDN
The paper has results for Blocks and HeaviSine at SNR = 40
Run publish_results.m in the misc folder to generate pdfs presented int the paper.
NOTE: While preparing the code for release, I realized that I used the average
of SER in db when I compiled results for publication, which is probably
not the correct way to compute the average.
The average SER in the paper was computed as mean(-10log10(normalized error))
The standard way to compute average is -10log10(mean(normalized squared error))
where mean is over all the experiments.
To reproduce results in the paper, use publish_results.m
Otherwise use publish_results_new.m
This changes the values of average SER for all the solvers. However, the
comparison among the solvers remains almost same.0001 % Comparison of various solvers for iterative reweighting and adaptive reweighting 0002 % 0003 % Solves the following basis pursuit denoising (BPDN) problem 0004 % min_x \Sum \w_i |x_i| + 1/2*||y-Ax||_2^2 0005 % 0006 % while updating the weights w_i 0007 % 0008 % Written by: Salman Asif, Georgia Tech 0009 % Email: sasif@gatech.edu 0010 % Created: June 16, 2011 0011 % 0012 % Reference: 0013 % "Fast and accurate algorithms for re-weighted L1 norm minimization," by 0014 % M. Salman Asif and Justin Romberg 0015 % 0016 % To reproduce experiments in the paper, use this scrip with 0017 % 0018 % rseed = 2012; 0019 % rand('state',rseed); 0020 % randn('state',rseed); 0021 % 0022 % before running each simulation in 0023 % script_simulation_wtBPDN 0024 % 0025 % The paper has results for Blocks and HeaviSine at SNR = 40 0026 % Run publish_results.m in the misc folder to generate pdfs presented int the paper. 0027 % 0028 % NOTE: While preparing the code for release, I realized that I used the average 0029 % of SER in db when I compiled results for publication, which is probably 0030 % not the correct way to compute the average. 0031 % 0032 % The average SER in the paper was computed as mean(-10log10(normalized error)) 0033 % The standard way to compute average is -10log10(mean(normalized squared error)) 0034 % where mean is over all the experiments. 0035 % 0036 % To reproduce results in the paper, use publish_results.m 0037 % Otherwise use publish_results_new.m 0038 % 0039 % This changes the values of average SER for all the solvers. However, the 0040 % comparison among the solvers remains almost same. 0041 0042 0043 % function job_wtBPDN_WAVE(mT, sT, snr, rwt_mode, lam_mode) 0044 % mT = 1, sT = 2, snr = 3, rwt_mode = 5, lam_mode = 1; 0045 0046 clear 0047 0048 EXP_LIST = [1 1 2 5 1; 1 1 3 5 1; 1 1 4 5 1; 1 2 2 5 1; 1 2 3 5 1; 1 2 4 5 1]; 0049 maxNumCompThreads(1); 0050 % numCores = str2double(getenv('NUMBER_OF_PROCESSORS')); 0051 % mpSize = 1;% size(EXP_LIST,1); 0052 % if matlabpool('size') ~= mpSize 0053 % if matlabpool('size')~=0 0054 % matlabpool close; 0055 % else 0056 % matlabpool('open', mpSize); 0057 % end 0058 % end 0059 0060 %% Setup path 0061 mname = mfilename; 0062 mpath = mfilename('fullpath'); 0063 mdir = mpath(1:end-length(mname)); 0064 cd(mdir); 0065 0066 addpath ../Pursuits_Homotopy/ 0067 addpath ../utils/ 0068 addpath ../utils/utils_Wavelet 0069 addpath ../solvers/ 0070 addpath src/ 0071 0072 %% Simulation parameters 0073 % reweighted setup 0074 rwt = 5; % number of reweighting iterations 0075 rwt_adp = 0; % number of reweighting iterations after adaptive reweighting 0076 0077 % rank-1 update mode 0078 delx_mode = 'mil'; % mil or qr 0079 0080 % simulation setup 0081 maxsim = 100; 0082 0083 mType_list = {'randn','orth'}; 0084 sType_list = {'blocks','HeaviSine','pcwPoly'}; 0085 SNR_list = [20:10:40 inf]; 0086 0087 lambda_list = [0, 1e-1, 1e-2, 1e-4]; 0088 % multiplication factor for tau = lambda *\|A'y\|_\infty 0089 % 0 --> sigma*log(N) 0090 0091 N_list = [256 512 1024]; 0092 M_ratio = [2:.5:4]; 0093 0094 for pf = 1:length(EXP_LIST) 0095 0096 mT = EXP_LIST(pf,1); 0097 sT = EXP_LIST(pf,2); 0098 snr = EXP_LIST(pf,3); 0099 rwt_mode = EXP_LIST(pf,4); 0100 lam_mode = EXP_LIST(pf,5); 0101 0102 EXP_stack = {}; 0103 estack = 1; 0104 EXP_stack{1,1} = 'mType'; 0105 EXP_stack{1,2} = 'sType'; 0106 EXP_stack{1,3} = 'SNR'; 0107 EXP_stack{1,4} = 'rwt_mode'; 0108 EXP_stack{1,5} = 'lambda_mode'; 0109 EXP_stack{1,6} = 'N'; 0110 EXP_stack{1,7} = 'M'; 0111 EXP_stack{1,8} = 'T'; 0112 EXP_stack{1,9} = 'str0'; 0113 EXP_stack{1,10} = 'str2'; 0114 EXP_stack{1,11} = 'avg SIM_stack'; 0115 EXP_stack{1,12} = 'SIM_stack'; 0116 EXP_stack{1,13} = 'SIM_memory'; 0117 0118 % for mT = 1:length(mType_list); 0119 % for sT = 1:length(sType_list); 0120 % for snr = 1:length(SNR_list); 0121 % for lam_mode = 1:length(lambda_list) 0122 0123 for nl = 1:length(N_list) 0124 for mr = 1:length(M_ratio) 0125 0126 N = N_list(nl); % signal length 0127 M = round(N/M_ratio(mr)); % no. of measurements 0128 T = N; % sparsity level 0129 0130 sType = char(sType_list{sT}); 0131 mType = char(mType_list{mT}); 0132 SNR = SNR_list(snr); 0133 0134 lambda = lambda_list(lam_mode); 0135 0136 str0 = sprintf('mType-%s, sType-%s, SNR = %d, (N,M,T) = %d, %d, %d, rwt_mode-%d, lambda%3.4g.', mType, sType, SNR, N, M, T, rwt_mode, lambda); 0137 disp(str0); 0138 0139 % filename_save = sprintf('results_comparison_ALL/comparison_wtBPDN_mType%s_sType%s_SNR%d_rwt_mode%d_lambda%3.4g_N%d_M%d_T%d.mat',mType,sType,SNR,rwt_mode,lambda,N,M,T); 0140 % filename_save = sprintf('results_comparison_ALL/comparison_wtBPDN_mT-%s_sT-%s_SNR%d_rwt_mode%d_lam_mode%d.mat',mType,sType,SNR,rwt_mode,lam_mode); 0141 0142 filename_save = sprintf('results_comparison_ALL/comparison_wtBPDN_mT-%s_sT-%s_SNR%d-reproduce-Trwt.mat',mType,sType,SNR); 0143 % filename_save = sprintf('results_comparison_ALL/comparison_wtBPDN_mT-%s_sT-%s_SNR%d-reproduce-rwt.mat',mType,sType,SNR); 0144 0145 % load RandomStates 0146 % 0147 % Experiments in the earlier draft... 0148 % rseed = 0; 0149 % rseed = sum(100*clock); 0150 % RandStream.setDefaultStream(RandStream('mt19937ar','seed',rseed)); 0151 % 0152 % Experiments that reproduce results in the paper... 0153 rseed = 2012; 0154 rand('state',rseed); 0155 randn('state',rseed); 0156 0157 % rng(rseed); 0158 % rseed = rng(2012); 0159 0160 %% Run simulation 0161 script_simulation_wtBPDN 0162 0163 %% save results 0164 estack = estack+1; 0165 0166 EXP_stack{estack,1} = mType; 0167 EXP_stack{estack,2} = sType; 0168 EXP_stack{estack,3} = SNR; 0169 EXP_stack{estack,4} = rwt_mode; 0170 EXP_stack{estack,5} = lambda; 0171 EXP_stack{estack,6} = N; 0172 EXP_stack{estack,7} = M; 0173 EXP_stack{estack,8} = T; 0174 EXP_stack{estack,9} = str0; 0175 EXP_stack{estack,10} = str2; 0176 EXP_stack{estack,11} = mean(cell2mat(SIM_stack),1); 0177 EXP_stack{estack,12} = SIM_stack; 0178 EXP_stack{estack,13} = SIM_memory; 0179 0180 eval(sprintf('save %s', filename_save)); 0181 end 0182 end 0183 end