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