nikcleju@54: % File: study_analysis_rec_algos nikcleju@54: % Run experiment to prove that our approx. ABS approach really converges to nikcleju@54: % the anaysis solution as lambda -> infty nikcleju@54: % For this, we need to generate signals that, provably, yield bad results nikcleju@54: % with synthesis recovery, and good results with analysis recovery nikcleju@54: % To do this we choose, N >> d, small l, and m close to d nikcleju@54: nikcleju@54: clear all nikcleju@54: close all nikcleju@54: nikcleju@54: % ================================= nikcleju@54: % Set up experiment parameters nikcleju@54: %================================== nikcleju@54: % Which form factor, delta and rho we want nikcleju@54: sigma = 2; nikcleju@54: %delta = 0.995; nikcleju@54: %rho = 0.7; nikcleju@54: %delta = 0.8; nikcleju@54: %rho = 0.15; nikcleju@54: delta = 0.5; nikcleju@54: rho = 0.05; nikcleju@54: nikcleju@54: % Number of vectors to generate each time nikcleju@54: numvects = 10; nikcleju@54: nikcleju@54: % Add noise nikcleju@54: % This is norm(signal)/norm(noise), so power, not energy nikcleju@54: SNRdb = 20; nikcleju@54: %epsextrafactor = 2 nikcleju@54: nikcleju@54: % ================================= nikcleju@54: % Processing the parameters nikcleju@54: %================================== nikcleju@54: nikcleju@54: % Compute noiselevel from db nikcleju@54: noiselevel = 1 / (10^(SNRdb/10)); nikcleju@54: nikcleju@54: % Set up parameter structure nikcleju@54: % and generate data X as well nikcleju@54: d = 50; nikcleju@54: p = round(sigma*d); nikcleju@54: m = round(delta*d); nikcleju@54: l = round(d - rho*m); nikcleju@54: nikcleju@54: % Generate Omega and data based on parameters nikcleju@54: Omega = Generate_Analysis_Operator(d, p); nikcleju@54: %Make Omega more coherent nikcleju@54: [U, S, V] = svd(Omega); nikcleju@54: %Sdnew = diag(S) ./ (1:numel(diag(S)))'; nikcleju@54: Sdnew = diag(S) .* (1:numel(diag(S)))'; % Make D coherent, not Omega! nikcleju@54: Snew = [diag(Sdnew); zeros(size(S,1) - size(S,2), size(S,2))]; nikcleju@54: Omega = U * Snew * V'; nikcleju@54: [x0,y,M,Lambda, realnoise] = Generate_Data_Known_Omega(Omega, d,p,m,l,noiselevel, numvects,'l0'); nikcleju@54: nikcleju@54: datafilename = 'mat/data_SNR20'; nikcleju@54: save(datafilename); nikcleju@54: %load mat/data_SNR20 nikcleju@54: nikcleju@54: % Values of lambda nikcleju@54: %lambdas = sqrt([1e-10 1e-8 1e-6 1e-4 1e-2 1 100 1e4 1e6 1e8 1e10]); nikcleju@54: lambdas = [0 10.^linspace(-5, 4, 10)]; nikcleju@54: %lambdas = 1000 nikcleju@54: %lambdas = sqrt([0 1 10000]); nikcleju@54: nikcleju@54: % Algorithm identifiers nikcleju@54: algonone = []; nikcleju@54: ompk = 1; nikcleju@54: ompeps = 2; nikcleju@54: tst = 3; nikcleju@54: bp = 4; nikcleju@54: gap = 5; nikcleju@54: nesta = 6; nikcleju@54: sl0 = 7; nikcleju@54: yall1 = 8; nikcleju@54: spgl1 = 9; nikcleju@54: nestasynth = 10; nikcleju@54: numallalgos = 10; nikcleju@54: algoname{ompk} = 'ABS OMPk'; nikcleju@54: algoname{ompeps} = 'ABS OMPeps'; nikcleju@54: algoname{tst} = 'ABS TST'; nikcleju@54: algoname{bp} = 'ABS BP'; nikcleju@54: algoname{gap} = 'GAP'; nikcleju@54: algoname{nesta} = 'NESTA Analysis'; nikcleju@54: algoname{sl0} = 'ABS SL0'; nikcleju@54: algoname{yall1} = 'ABS YALL1'; nikcleju@54: algoname{spgl1} = 'ABS SPGL1'; nikcleju@54: algoname{nestasynth} = 'NESTA Synth'; nikcleju@54: nikcleju@54: % What algos to run nikcleju@54: %algos = [gap, ompk, ompeps, tst, bp, sl0, yall1]; nikcleju@54: %algos = [gap, ompk, ompeps, tst, bp]; nikcleju@54: algos = [gap, sl0]; nikcleju@54: numalgos = numel(algos); nikcleju@54: nikcleju@54: % Save mat file? nikcleju@54: do_save_mat = 0; nikcleju@54: matfilename = 'mat/approx_d50_sigma2_SNR20db_Dcoh'; nikcleju@54: nikcleju@54: % Save figures? nikcleju@54: do_save_figs = 0; nikcleju@54: figfilename = 'figs/approx_d50_sigma2_SNR20db_Dcoh'; nikcleju@54: nikcleju@54: % Show progressbar ? (not recommended when running on parallel threads) nikcleju@54: do_progressbar = 0; nikcleju@54: if do_progressbar nikcleju@54: progressbar('Total', 'Current parameters'); nikcleju@54: end nikcleju@54: nikcleju@54: % Init times nikcleju@54: for i = 1:numalgos nikcleju@54: elapsed(algos(i)) = 0; nikcleju@54: end nikcleju@54: nikcleju@54: % Init results structure nikcleju@54: results = []; nikcleju@54: nikcleju@54: nikcleju@54: % ======== nikcleju@54: % Run nikcleju@54: % ======== nikcleju@54: nikcleju@54: % Run GAP and NESTA first nikcleju@54: if any(algos == gap) nikcleju@54: for iy = 1:size(y,2) nikcleju@54: a = gap; nikcleju@54: % Compute epsilon nikcleju@54: %epsilon = epsextrafactor * noiselevel * norm(y(:,iy)); nikcleju@54: epsilon = 1.1 * norm(realnoise(:,iy)); nikcleju@54: % nikcleju@54: gapparams = []; nikcleju@54: gapparams.num_iteration = 1000; nikcleju@54: gapparams.greedy_level = 0.9; nikcleju@54: gapparams.stopping_coefficient_size = 1e-4; nikcleju@54: gapparams.l2solver = 'pseudoinverse'; nikcleju@54: %gapparams.noise_level = noiselevel; nikcleju@54: gapparams.noise_level = epsilon; nikcleju@54: timer(a) = tic; nikcleju@54: xrec{a}(:,iy) = GAP(y(:,iy), M, M', Omega, Omega', gapparams, zeros(d,1)); nikcleju@54: elapsed(a) = elapsed(a) + toc(timer(a)); nikcleju@54: % nikcleju@54: err{a}(iy) = norm(x0(:,iy) - xrec{a}(:,iy)); nikcleju@54: relerr{a}(iy) = err{a}(iy) / norm(x0(:,iy)); nikcleju@54: end nikcleju@54: disp([ algoname{a} ': avg relative error = ' num2str(mean(relerr{a}))]); nikcleju@54: end nikcleju@54: if any(algos == nesta) nikcleju@54: for iy = 1:size(y,2) nikcleju@54: a = nesta; nikcleju@54: % Compute epsilon nikcleju@54: %epsilon = epsextrafactor * noiselevel * norm(y(:,iy)); nikcleju@54: epsilon = 1.1 * norm(realnoise(:,iy)); nikcleju@54: % nikcleju@54: try nikcleju@54: timer(a) = tic; nikcleju@54: %xrec{a}(:,iy) = do_nesta_DemoNonProjector(x0(:,iy), M, Omega', 0); nikcleju@54: nikcleju@54: % Common heuristic: delta = sqrt(m + 2*sqrt(2*m))*sigma ????????? nikcleju@54: [U,S,V] = svd(M,'econ'); nikcleju@54: opts.U = Omega; nikcleju@54: opts.Ut = Omega'; nikcleju@54: opts.USV.U=U; nikcleju@54: opts.USV.S=S; nikcleju@54: opts.USV.V=V; nikcleju@54: opts.TolVar = 1e-5; nikcleju@54: opts.Verbose = 0; nikcleju@54: xrec{a}(:,iy) = NESTA(M, [], y(:,iy), 1e-3, epsilon, opts); nikcleju@54: elapsed(a) = elapsed(a) + toc(timer(a)); nikcleju@54: catch err nikcleju@54: disp('*****ERROR: NESTA throwed error *****'); nikcleju@54: xrec{a}(:,iy) = zeros(size(x0(:,iy))); nikcleju@54: end nikcleju@54: % nikcleju@54: err{a}(iy) = norm(x0(:,iy) - xrec{a}(:,iy)); nikcleju@54: relerr{a}(iy) = err{a}(iy) / norm(x0(:,iy)); nikcleju@54: end nikcleju@54: disp([ algoname{a} ': avg relative error = ' num2str(mean(relerr{a}))]); nikcleju@54: end nikcleju@54: nikcleju@54: for i = 1:numel(algos) nikcleju@54: if algos(i) == gap nikcleju@54: continue; nikcleju@54: end nikcleju@54: prevgamma{algos(i)} = zeros(p, numvects); nikcleju@54: end nikcleju@54: nikcleju@54: % Run ABS algorithms (lambda-dependent) nikcleju@54: for iparam = 1:numel(lambdas) nikcleju@54: nikcleju@54: % Read current parameters nikcleju@54: lambda = lambdas(iparam); nikcleju@54: nikcleju@54: % Init stuff nikcleju@54: for i = 1:numel(algos) nikcleju@54: if algos(i) == gap || algos(i) == nesta nikcleju@54: continue; nikcleju@54: end nikcleju@54: xrec{algos(i)} = zeros(d, numvects); nikcleju@54: end nikcleju@54: % nikcleju@54: for i = 1:numel(algos) nikcleju@54: if algos(i) == gap || algos(i) == nesta nikcleju@54: continue; nikcleju@54: end nikcleju@54: err{algos(i)} = zeros(numvects,1); nikcleju@54: relerr{algos(i)} = zeros(numvects,1); nikcleju@54: end nikcleju@54: nikcleju@54: % For every generated signal do nikcleju@54: for iy = 1:size(y,2) nikcleju@54: nikcleju@54: % Compute epsilon nikcleju@54: %epsilon = epsextrafactor * noiselevel * norm(y(:,iy)); nikcleju@54: epsilon = 1.1 * norm(realnoise(:,iy)); nikcleju@54: nikcleju@54: %-------------------------------- nikcleju@54: % Reconstruct (and measure delay), Compute reconstruction error nikcleju@54: %-------------------------------- nikcleju@54: for i = 1:numel(algos) nikcleju@54: a = algos(i); nikcleju@54: if a == gap || a == nesta nikcleju@54: continue nikcleju@54: end nikcleju@54: if a == ompk nikcleju@54: timer(a) = tic; nikcleju@54: xrec{a}(:,iy) = ABS_OMPk_approx(y(:,iy), Omega, M, p-l, lambda); nikcleju@54: elapsed(a) = elapsed(a) + toc(timer(a)); nikcleju@54: elseif a == ompeps nikcleju@54: timer(a) = tic; nikcleju@54: xrec{a}(:,iy) = ABS_OMPeps_approx(y(:,iy), Omega, M, epsilon, lambda, Omega * x0(:,iy)); nikcleju@54: elapsed(a) = elapsed(a) + toc(timer(a)); nikcleju@54: elseif a == tst nikcleju@54: timer(a) = tic; nikcleju@54: xrec{a}(:,iy) = ABS_TST_approx(y(:,iy), Omega, M, epsilon, lambda, Omega * x0(:,iy)); nikcleju@54: elapsed(a) = elapsed(a) + toc(timer(a)); nikcleju@54: elseif a == bp nikcleju@54: timer(a) = tic; nikcleju@54: [xrec{a}(:,iy), prevgamma{a}(:,iy)] = ABS_BP_approx(y(:,iy), Omega, M, epsilon, lambda, prevgamma{a}(:,iy), Omega * x0(:,iy)); nikcleju@54: elapsed(a) = elapsed(a) + toc(timer(a)); nikcleju@54: elseif a == yall1 nikcleju@54: timer(a) = tic; nikcleju@54: xrec{a}(:,iy) = ABS_YALL1_approx(y(:,iy), Omega, M, epsilon, lambda, Omega * x0(:,iy)); nikcleju@54: elapsed(a) = elapsed(a) + toc(timer(a)); nikcleju@54: % elseif a == gap nikcleju@54: % gapparams = []; nikcleju@54: % gapparams.num_iteration = 40; nikcleju@54: % gapparams.greedy_level = 0.9; nikcleju@54: % gapparams.stopping_coefficient_size = 1e-4; nikcleju@54: % gapparams.l2solver = 'pseudoinverse'; nikcleju@54: % %gapparams.noise_level = noiselevel; nikcleju@54: % gapparams.noise_level = epsilon; nikcleju@54: % timer(a) = tic; nikcleju@54: % xrec{a}(:,iy) = GAP(y(:,iy), M, M', Omega, Omega', gapparams, zeros(d,1)); nikcleju@54: % elapsed(a) = elapsed(a) + toc(timer(a)); nikcleju@54: % elseif a == nesta nikcleju@54: % try nikcleju@54: % timer(a) = tic; nikcleju@54: % %xrec{a}(:,iy) = do_nesta_DemoNonProjector(x0(:,iy), M, Omega', 0); nikcleju@54: % nikcleju@54: % % Common heuristic: delta = sqrt(m + 2*sqrt(2*m))*sigma ????????? nikcleju@54: % [U,S,V] = svd(M,'econ'); nikcleju@54: % opts.U = Omega; nikcleju@54: % opts.Ut = Omega'; nikcleju@54: % opts.USV.U=U; nikcleju@54: % opts.USV.S=S; nikcleju@54: % opts.USV.V=V; nikcleju@54: % opts.TolVar = 1e-5; nikcleju@54: % opts.Verbose = 0; nikcleju@54: % xrec{a}(:,iy) = NESTA(M, [], y(:,iy), 1e-3, epsilon, opts); nikcleju@54: % elapsed(a) = elapsed(a) + toc(timer(a)); nikcleju@54: % catch err nikcleju@54: % disp('*****ERROR: NESTA throwed error *****'); nikcleju@54: % xrec{a}(:,iy) = zeros(size(x0(:,iy))); nikcleju@54: % end nikcleju@54: elseif a == sl0 nikcleju@54: timer(a) = tic; nikcleju@54: xrec{a}(:,iy) = ABS_SL0_approx(y(:,iy), Omega, M, epsilon, lambda); nikcleju@54: elapsed(a) = elapsed(a) + toc(timer(a)); nikcleju@54: elseif a == spgl1 nikcleju@54: timer(a) = tic; nikcleju@54: xrec{a}(:,iy) = ABS_SPGL1_approx(y(:,iy), Omega, M, epsilon, lambda, Omega * xrec{bp}(:,iy)); nikcleju@54: elapsed(a) = elapsed(a) + toc(timer(a)); nikcleju@54: elseif a == nestasynth nikcleju@54: timer(a) = tic; nikcleju@54: xrec{a}(:,iy) = ABS_NESTAsynth_approx(y(:,iy), Omega, M, epsilon, lambda, Omega * xrec{bp}(:,iy)); nikcleju@54: elapsed(a) = elapsed(a) + toc(timer(a)); nikcleju@54: else nikcleju@54: %error('No such algorithm!'); nikcleju@54: end nikcleju@54: % nikcleju@54: % Compare to GAP instead! nikcleju@54: %x0(:,iy) = xrec{gap}(:,iy); nikcleju@54: % nikcleju@54: err{a}(iy) = norm(x0(:,iy) - xrec{a}(:,iy)); nikcleju@54: relerr{a}(iy) = err{a}(iy) / norm(x0(:,iy)); nikcleju@54: end nikcleju@54: nikcleju@54: % Update progressbar nikcleju@54: % if do_progressbar nikcleju@54: % %frac2 = iy/numvects; nikcleju@54: % %frac1 = ((iparam-1) + frac2) / count; nikcleju@54: % if norm(frac2 - 1) < 1e-6 nikcleju@54: % frac2 = 0; nikcleju@54: % end nikcleju@54: % frac2 = frac2 + incr2; nikcleju@54: % frac1 = frac1 + incr1; nikcleju@54: % progressbar(frac1, frac2); nikcleju@54: % end nikcleju@54: end nikcleju@54: nikcleju@54: %-------------------------------- nikcleju@54: % Save results in big stucture & display nikcleju@54: %-------------------------------- nikcleju@54: % Save reconstructed signals nikcleju@54: % Save rel & abs errors nikcleju@54: % Display error nikcleju@54: results(iparam).xrec = xrec; nikcleju@54: results(iparam).err = err; nikcleju@54: results(iparam).relerr = relerr; nikcleju@54: % nikcleju@54: %disp(['Simulation no. ' num2str(iparam)]); nikcleju@54: disp(['Lambda = ' num2str(lambda) ':']); nikcleju@54: for i = 1:numalgos nikcleju@54: a = algos(i); nikcleju@54: if a == gap || a == nesta nikcleju@54: continue nikcleju@54: end nikcleju@54: disp([ algoname{a} ': avg relative error = ' num2str(mean(relerr{a}))]); nikcleju@54: end nikcleju@54: end nikcleju@54: nikcleju@54: % ================================= nikcleju@54: % Save nikcleju@54: % ================================= nikcleju@54: if do_save_mat nikcleju@54: save(matfilename); nikcleju@54: disp(['Saved to ' matfilename]); nikcleju@54: end nikcleju@54: nikcleju@54: % ================================= nikcleju@54: % Plot nikcleju@54: % ================================= nikcleju@54: toplot = zeros(numel(lambdas), numalgos); nikcleju@54: nikcleju@54: relerrs = {results.relerr}; nikcleju@54: for i = 1:numalgos nikcleju@54: for j = 1:numel(lambdas) nikcleju@54: toplot(j,i) = mean(relerrs{j}{algos(i)}); nikcleju@54: end nikcleju@54: end nikcleju@54: nikcleju@54: %h = plot(toplot); nikcleju@54: h = semilogx(lambdas, toplot); nikcleju@54: legend(algoname{algos}) nikcleju@54: xlabel('Lambda') nikcleju@54: ylabel('Average reconstruction error') nikcleju@54: title('Reconstruction error with different algorithms') nikcleju@54: nikcleju@54: if (do_save_figs) nikcleju@54: saveas(gcf, [figfilename '.fig']); nikcleju@54: saveas(gcf, [figfilename '.png']); nikcleju@54: saveTightFigure(gcf, [figfilename '.pdf']); nikcleju@54: end nikcleju@54: