Mercurial > hg > smallbox

comparison toolboxes/alps/ALPSbenchmark_errornorm_vs_iter.m @ 154:0de08f68256b ivand_dev

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
ALPS toolbox - Algebraic Pursuit added to smallbox
author Ivan Damnjanovic lnx <ivan.damnjanovic@eecs.qmul.ac.uk>
date Fri, 12 Aug 2011 11:17:47 +0100
parents
children
comparison
equal deleted inserted replaced
153:af307f247ac7 154:0de08f68256b
1 % =========================================================================
2 % Algebraic Pursuit algorithms - Demo
3 % =========================================================================
4 % Algebraic Pursuit (ALPS) algorithms are accelerated hard thresholding methods
5 % on sparse recovery of linear inverse systems. In particular, let
6 % Phi : M x N real matrix (M < N)
7 % x* : N x 1 K-sparse data vector
8 % n : N x 1 additive noise vector
9 % then, y = Phi x* + n is the undersampled M x 1 measurement vector. ALPS
10 % solve the following minimization problem
11 % minimize ||y - Phi x||^2 subject to x is K-sparse vector.
12 %
13 % Detailed discussion on the algorithm can be found in
14 % [1] "On Accelerated Hard Thresholding Methods for Sparse Approximation", written
15 % by Volkan Cevher, Technical Report, 2011.
16 % =========================================================================
17 % 01/04/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL.
18 % =========================================================================
19
20 clear all
21 close all
22 clc;
23
24 addpath ALPS/;
25
26 % General parameters
27 N = 4000; % Size of sparse vector x
28 M = 800; % Number of measurements
29 K = ceil(M*(0.2)); % Sparsity of x*
30 sigma = 0.0000; % Additive noise variance
31 verbose = 1; % Print information
32
33 % ALPS parameters
34 its = 600; % Maximum number of iterations
35 Ttol= 1e-6; % Convergence tolerance
36 tau = -1; % Momentum term
37 mu = 0; % Step size selection
38 MC = 10; % Number of Monte-Carlo iterations
39 useCG = 0; % Conjugate Gradient
40
41 % Please refer to generate_matrix.m and generate_vector.m files
42 m_ensemble = 'Gaussian'; % Measurement matrix ensemble type
43 x_ensemble = 'Gaussian'; % Sparse vector ensemble type
44
45 max_its = 0;
46
47 for i = 1:MC
48 x = generate_vector(N, K, x_ensemble);
49 Phi = generate_matrix(M, N, m_ensemble);
50 noise = sigma*randn(M,1);
51 y = Phi*x + noise;
52
53 % 0-ALPS(0)
54 disp('=====================================================================');
55 [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 0, 'mode', 0, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu);
56 str = sprintf('Error norm: %f ', norm(x-x_hat));
57 disp(str);
58 num_iter(1, i) = numiter;
59 error(1,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
60 total_time(1, i) = time;
61
62 if (numiter > max_its)
63 max_its = numiter;
64 end;
65
66 % 0-ALPS(2)
67 disp('=====================================================================');
68 [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 0, 'mode', 2, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu);
69 str = sprintf('Error norm: %f ', norm(x-x_hat));
70 disp(str);
71 num_iter(2, i) = numiter;
72 error(2,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
73 total_time(2, i) = time;
74
75 if (numiter > max_its)
76 max_its = numiter;
77 end;
78
79 % 0-ALPS(4)
80 disp('=====================================================================');
81 [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 0, 'mode', 4, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu);
82 str = sprintf('Error norm: %f ', norm(x-x_hat));
83 disp(str);
84 num_iter(3, i) = numiter;
85 error(3,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
86 total_time(3, i) = time;
87
88 if (numiter > max_its)
89 max_its = numiter;
90 end;
91
92 % 0-ALPS(5)
93 disp('=====================================================================');
94 [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 0, 'mode', 5, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu);
95 str = sprintf('Error norm: %f ', norm(x-x_hat));
96 disp(str);
97 num_iter(4, i) = numiter;
98 error(4,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
99 total_time(4, i) = time;
100
101 if (numiter > max_its)
102 max_its = numiter;
103 end;
104
105 % 1-ALPS(2)
106 disp('=====================================================================');
107 [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 1, 'mode', 2, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu);
108 str = sprintf('Error norm: %f ', norm(x-x_hat));
109 disp(str);
110 num_iter(5, i) = numiter;
111 error(5,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
112 total_time(5, i) = time;
113
114 if (numiter > max_its)
115 max_its = numiter;
116 end;
117
118 % 1-ALPS(4)
119 disp('=====================================================================');
120 [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 1, 'mode', 4, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu);
121 str = sprintf('Error norm: %f ', norm(x-x_hat));
122 disp(str);
123 num_iter(6, i) = numiter;
124 error(6,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
125 total_time(6, i) = time;
126
127 if (numiter > max_its)
128 max_its = numiter;
129 end;
130
131 % 1-ALPS(5)
132 disp('=====================================================================');
133 [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 1, 'mode', 5, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu);
134 str = sprintf('Error norm: %f ', norm(x-x_hat));
135 disp(str);
136 num_iter(7, i) = numiter;
137 error(7,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
138 total_time(7,i) = time;
139
140 if (numiter > max_its)
141 max_its = numiter;
142 end;
143
144 % 2-ALPS(2)
145 disp('=====================================================================');
146 [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 2, 'mode', 2, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu);
147 str = sprintf('Error norm: %f ', norm(x-x_hat));
148 disp(str);
149 num_iter(8, i) = numiter;
150 error(8,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
151 total_time(8, i) = time;
152
153 if (numiter > max_its)
154 max_its = numiter;
155 end;
156
157 % 2-ALPS(4)
158 disp('=====================================================================');
159 [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 2, 'mode', 4, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu);
160 str = sprintf('Error norm: %f ', norm(x-x_hat));
161 disp(str);
162 num_iter(9, i) = numiter;
163 error(9,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
164 total_time(9, i) = time;
165
166 if (numiter > max_its)
167 max_its = numiter;
168 end;
169
170 % 2-ALPS(5)
171 disp('=====================================================================');
172 [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 2, 'mode', 5, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu);
173 str = sprintf('Error norm: %f ', norm(x-x_hat));
174 disp(str);
175 num_iter(10, i) = numiter;
176 error(10,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
177 total_time(10, i) = time;
178
179 if (numiter > max_its)
180 max_its = numiter;
181 end;
182
183 % 3-ALPS(2)
184 disp('=====================================================================');
185 [x_hat, numiter, time, x_path] = AlgebraicPursuit(y, Phi, K, 'memory', 3, 'mode', 2, 'tolerance', Ttol, 'ALPSiterations', its, 'verbose', verbose, 'tau', tau, 'useCG', useCG, 'mu', mu);
186 str = sprintf('Error norm: %f ', norm(x-x_hat));
187 disp(str);
188 num_iter(11, i) = numiter;
189 error(11,i,1:size(x_path,2))=sqrt(sum(abs(x*ones(1,size(x_path,2))-x_path).^2));
190 total_time(11, i) = time;
191
192 if (numiter > max_its)
193 max_its = numiter;
194 end;
195 end;
196
197 set(0,'DefaultAxesFontSize',12);
198
199 figure;
200 semilogy(squeeze(median(error(1,:,:),2)),'LineWidth',3); hold on
201 semilogy(squeeze(median(error(2,:,:),2)),'--','LineWidth',3);
202 semilogy(squeeze(median(error(3,:,:),2)),':','LineWidth',3);
203 semilogy(squeeze(median(error(4,:,:),2)),'k','LineWidth',4);
204 semilogy(squeeze(median(error(5,:,:),2)),'k--','LineWidth',4);
205 semilogy(squeeze(median(error(6,:,:),2)),'k:','LineWidth',4);
206 semilogy(squeeze(median(error(7,:,:),2)),'k-.','LineWidth',4);
207 semilogy(squeeze(median(error(8,:,:),2)),'r--','LineWidth',3);
208 semilogy(squeeze(median(error(9,:,:),2)),'m--','LineWidth',3);
209 semilogy(squeeze(median(error(10,:,:),2)),'k--','LineWidth',3);
210 semilogy(squeeze(median(error(11,:,:),2)),'g--','LineWidth',3);
211
212 xlabel('[i]:iteration number')
213 ylabel('||x-x_i||')
214 axis([1 max_its + 10 10^-5 2*10^0])
215
216 legend(strcat('0-IHT(0) [', num2str(mean(num_iter(1,:))),'][', num2str(mean(total_time(1,:))),']'), ...
217 strcat('0-IHT(2) [', num2str(mean(num_iter(2,:))),'][', num2str(mean(total_time(2,:))),']'), ...
218 strcat('0-IHT(4) [', num2str(mean(num_iter(3,:))),'][', num2str(mean(total_time(3,:))),']'), ...
219 strcat('0-IHT(5) [', num2str(mean(num_iter(4,:))),'][', num2str(mean(total_time(4,:))),']'), ...
220 strcat('1-IHT(2) [', num2str(mean(num_iter(5,:))),'][', num2str(mean(total_time(5,:))),']'), ...
221 strcat('1-IHT(4) [', num2str(mean(num_iter(6,:))),'][', num2str(mean(total_time(6,:))),']'), ...
222 strcat('1-IHT(5) [', num2str(mean(num_iter(7,:))),'][', num2str(mean(total_time(7,:))),']'), ...
223 strcat('2-ALPS(2) [', num2str(mean(num_iter(8,:))),'][', num2str(mean(total_time(8,:))),']'), ...
224 strcat('2-ALPS(4) [', num2str(mean(num_iter(9,:))),'][', num2str(mean(total_time(9,:))),']'), ...
225 strcat('2-ALPS(5) [', num2str(mean(num_iter(10,:))),'][', num2str(mean(total_time(10,:))),']'), ...
226 strcat('3-ALPS(2) [', num2str(mean(num_iter(11,:))),'][', num2str(mean(total_time(11,:))),']'));
227 title( strcat('N= ', num2str(N), ', M= ', num2str(M), ', K= ', num2str(K)) );
228 shg;