ivan@154
|
1 function [x_hat, numiter, x_path] = infty_ALPS(y, Phi, K, params)
|
ivan@154
|
2 % =========================================================================
|
ivan@154
|
3 % infty-ALPS(#) algorithm - Beta Version
|
ivan@154
|
4 % =========================================================================
|
ivan@154
|
5 % Algebraic Pursuit (ALPS) algorithm with infty-memory acceleration.
|
ivan@154
|
6 %
|
ivan@154
|
7 % Detailed discussion on the algorithm can be found in
|
ivan@154
|
8 % [1] "On Accelerated Hard Thresholding Methods for Sparse Approximation", written
|
ivan@154
|
9 % by Volkan Cevher, Technical Report, 2011.
|
ivan@154
|
10 % =========================================================================
|
ivan@154
|
11 % INPUT ARGUMENTS:
|
ivan@154
|
12 % y M x 1 undersampled measurement vector.
|
ivan@154
|
13 % Phi M x N regression matrix.
|
ivan@154
|
14 % K Sparsity of underlying vector x* or desired
|
ivan@154
|
15 % sparsity of solution.
|
ivan@154
|
16 % params Structure of parameters. These are:
|
ivan@154
|
17 %
|
ivan@154
|
18 % tol,... Early stopping tolerance. Default value: tol =
|
ivan@154
|
19 % 1-e5
|
ivan@154
|
20 % ALPSiters,... Maximum number of algorithm iterations. Default
|
ivan@154
|
21 % value: 300.
|
ivan@154
|
22 % solveNewtonb,... Value: solveNewtonb = 0. Not used in infty
|
ivan@154
|
23 % methods.
|
ivan@154
|
24 % gradientDescentx,... If gradientDescentx == 1: single gradient
|
ivan@154
|
25 % update of x_{i+1} restricted ot its support with
|
ivan@154
|
26 % line search. Default value: gradientDescentx =
|
ivan@154
|
27 % 1.
|
ivan@154
|
28 % solveNewtonx,... If solveNewtonx == 1: Akin to Hard Thresholding Pursuit
|
ivan@154
|
29 % (c.f. Simon Foucart, "Hard Thresholding Pursuit,"
|
ivan@154
|
30 % preprint, 2010). Default vale: solveNewtonx = 0
|
ivan@154
|
31 % tau,... Variable that controls the momentum in
|
ivan@154
|
32 % non-memoryless case. Ignored in memoryless
|
ivan@154
|
33 % case. Default value: tau = 1/2.
|
ivan@154
|
34 % Special cases:
|
ivan@154
|
35 % - tau = -1: momentum step size is automatically
|
ivan@154
|
36 % optimized in every step.
|
ivan@154
|
37 % - tau as a function handle: user defined
|
ivan@154
|
38 % behavior of tau momentum term.
|
ivan@154
|
39 % mu,... Variable that controls the step size selection.
|
ivan@154
|
40 % When mu = 0, step size is computed adaptively
|
ivan@154
|
41 % per iteration. Default value: mu = 0.
|
ivan@154
|
42 % cg_maxiter,... Maximum iterations for Conjugate-Gradients method.
|
ivan@154
|
43 % cg_tol Tolerance variable for Conjugate-Gradients method.
|
ivan@154
|
44 % =========================================================================
|
ivan@154
|
45 % OUTPUT ARGUMENTS:
|
ivan@154
|
46 % x_hat N x 1 recovered K-sparse vector.
|
ivan@154
|
47 % numiter Number of iterations executed.
|
ivan@154
|
48 % x_path Keeps a series of computed N x 1 K-sparse vectors
|
ivan@154
|
49 % until the end of the iterative process.
|
ivan@154
|
50 % =========================================================================
|
ivan@154
|
51 % 01/04/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL.
|
ivan@154
|
52 % =========================================================================
|
ivan@154
|
53 % cgsolve.m is written by Justin Romberg, Caltech, Oct. 2005.
|
ivan@154
|
54 % Email: jrom@acm.caltech.edu
|
ivan@154
|
55 % =========================================================================
|
ivan@154
|
56 % This work was supported in part by the European Commission under Grant
|
ivan@154
|
57 % MIRG-268398 and DARPA KeCoM program #11-DARPA-1055. VC also would like
|
ivan@154
|
58 % to acknowledge Rice University for his Faculty Fellowship.
|
ivan@154
|
59 % =========================================================================
|
ivan@154
|
60
|
ivan@154
|
61 [~,N] = size(Phi);
|
ivan@154
|
62
|
ivan@154
|
63 %% Initialize transpose of measurement matrix
|
ivan@154
|
64
|
ivan@154
|
65 Phi_t = Phi';
|
ivan@154
|
66
|
ivan@154
|
67 %% Initialize to zero vector
|
ivan@154
|
68 x_cur = zeros(N,1);
|
ivan@154
|
69 y_cur = zeros(N,1);
|
ivan@154
|
70 X_i = [];
|
ivan@154
|
71
|
ivan@154
|
72 x_path = zeros(N, params.ALPSiters);
|
ivan@154
|
73
|
ivan@154
|
74 %% CG params
|
ivan@154
|
75 if (params.solveNewtonx == 1 || params.solveNewtonb == 1)
|
ivan@154
|
76 cg_verbose = 0;
|
ivan@154
|
77 cg_A = Phi_t*Phi;
|
ivan@154
|
78 cg_b = Phi_t*y;
|
ivan@154
|
79 end;
|
ivan@154
|
80
|
ivan@154
|
81 %% Determine momentum step size selection strategy
|
ivan@154
|
82 optimizeTau = 0;
|
ivan@154
|
83 function_tau = 0;
|
ivan@154
|
84
|
ivan@154
|
85 if (isa(params.tau,'float'))
|
ivan@154
|
86 if (params.tau == -1)
|
ivan@154
|
87 optimizeTau = 1;
|
ivan@154
|
88 end;
|
ivan@154
|
89 elseif (isa(params.tau, 'function_handle'))
|
ivan@154
|
90 function_tau = 1;
|
ivan@154
|
91 end;
|
ivan@154
|
92
|
ivan@154
|
93 %% Determine step size selection strategy
|
ivan@154
|
94 function_mu = 0;
|
ivan@154
|
95 adaptive_mu = 0;
|
ivan@154
|
96
|
ivan@154
|
97 if (isa(params.mu,'float'))
|
ivan@154
|
98 function_mu = 0;
|
ivan@154
|
99 if (params.mu == 0)
|
ivan@154
|
100 adaptive_mu = 1;
|
ivan@154
|
101 else
|
ivan@154
|
102 adaptive_mu = 0;
|
ivan@154
|
103 end;
|
ivan@154
|
104 elseif (isa(params.mu,'function_handle'))
|
ivan@154
|
105 function_mu = 1;
|
ivan@154
|
106 end;
|
ivan@154
|
107
|
ivan@154
|
108 %% Help variables
|
ivan@154
|
109 complementary_Xi = ones(N,1);
|
ivan@154
|
110 setXi = zeros(N,1);
|
ivan@154
|
111 setYi = zeros(N,1);
|
ivan@154
|
112
|
ivan@154
|
113 i = 1;
|
ivan@154
|
114 %% infty-ALPS(#)
|
ivan@154
|
115 while (i <= params.ALPSiters)
|
ivan@154
|
116 x_path(:,i) = x_cur;
|
ivan@154
|
117 x_prev = x_cur;
|
ivan@154
|
118
|
ivan@154
|
119 % Compute the residual
|
ivan@154
|
120 if (i == 1)
|
ivan@154
|
121 res = y;
|
ivan@154
|
122 else
|
ivan@154
|
123 Phi_x_cur = Phi(:,X_i)*x_cur(X_i);
|
ivan@154
|
124 res = y - Phi_x_cur;
|
ivan@154
|
125 end;
|
ivan@154
|
126
|
ivan@154
|
127 % Compute the derivative
|
ivan@154
|
128 der = Phi_t*res;
|
ivan@154
|
129
|
ivan@154
|
130 % Determine S_i set via eq. (11)
|
ivan@154
|
131 complementary_Xi(X_i) = 0;
|
ivan@154
|
132 [~, ind_der] = sort(abs(der).*complementary_Xi, 'descend');
|
ivan@154
|
133 complementary_Xi(X_i) = 1;
|
ivan@154
|
134 S_i = [X_i; ind_der(1:K)];
|
ivan@154
|
135
|
ivan@154
|
136 ider = der(S_i);
|
ivan@154
|
137
|
ivan@154
|
138 setder = zeros(N,1);
|
ivan@154
|
139 setder(S_i) = 1;
|
ivan@154
|
140 if (adaptive_mu)
|
ivan@154
|
141 % Step size selection via eq. (12) and eq. (13)
|
ivan@154
|
142 Pder = Phi(:,S_i)*ider;
|
ivan@154
|
143 mu_bar = ider'*ider/(Pder'*Pder);
|
ivan@154
|
144 end;
|
ivan@154
|
145
|
ivan@154
|
146 iy_cur = y_cur.*setYi;
|
ivan@154
|
147 if (~function_tau) % If tau is not a function handle...
|
ivan@154
|
148 if (optimizeTau) % Compute optimized tau
|
ivan@154
|
149
|
ivan@154
|
150 % tau = argmin || u - Phi(x_i + y_i) ||
|
ivan@154
|
151 % = <Phi*y_i, u - Phi(x_i - mu/2 * grad_Si f(xi))>/||Phi*y_i||^2
|
ivan@154
|
152
|
ivan@154
|
153 if (i == 1)
|
ivan@154
|
154 params.tau = 0;
|
ivan@154
|
155 else
|
ivan@154
|
156 % u - Phi*(x_i - mu/2 grad_Si f(xi)) = u - Phi*b
|
ivan@154
|
157 if (adaptive_mu)
|
ivan@154
|
158 b = x_cur(S_i) + mu_bar*ider; % Non-zero elems: S_i
|
ivan@154
|
159 elseif (function_mu)
|
ivan@154
|
160 b = x_cur(S_i) + params.mu(i)*ider;
|
ivan@154
|
161 else b = x_cur(S_i) + params.mu*ider;
|
ivan@154
|
162 end;
|
ivan@154
|
163
|
ivan@154
|
164 y_Phi_b = y - Phi(:,S_i)*b;
|
ivan@154
|
165 Phi_y_prev = Phi(:,Y_i)*y_cur(Y_i); % Phi * y_i
|
ivan@154
|
166 params.tau = y_Phi_b'*Phi_y_prev/(Phi_y_prev'*Phi_y_prev);
|
ivan@154
|
167 end;
|
ivan@154
|
168
|
ivan@154
|
169 if (adaptive_mu)
|
ivan@154
|
170 y_cur = params.tau*iy_cur + mu_bar*der.*setder;
|
ivan@154
|
171 elseif (function_mu)
|
ivan@154
|
172 y_cur = params.tau*iy_cur + params.mu(i)*der.*setder;
|
ivan@154
|
173 else y_cur = params.tau*iy_cur + params.mu*der.*setder;
|
ivan@154
|
174 end;
|
ivan@154
|
175
|
ivan@154
|
176 Y_i = ne(y_cur,0);
|
ivan@154
|
177 setYi = zeros(N,1);
|
ivan@154
|
178 setYi(Y_i) = 1;
|
ivan@154
|
179 else % Tau fixed and user-defined
|
ivan@154
|
180 if (adaptive_mu)
|
ivan@154
|
181 y_cur = params.tau*iy_cur + mu_bar*der.*setder;
|
ivan@154
|
182 elseif (function_mu)
|
ivan@154
|
183 y_cur = params.tau*iy_cur + params.mu(i)*der.*setder;
|
ivan@154
|
184 else y_cur = params.tau*iy_cur + params.mu*der.*setder;
|
ivan@154
|
185 end;
|
ivan@154
|
186
|
ivan@154
|
187 Y_i = ne(y_cur,0);
|
ivan@154
|
188 setYi = zeros(N,1);
|
ivan@154
|
189 setYi(Y_i) = 1;
|
ivan@154
|
190 end;
|
ivan@154
|
191 else
|
ivan@154
|
192 if (adaptive_mu)
|
ivan@154
|
193 y_cur = params.tau(i)*iy_cur + mu_bar*der.*setder;
|
ivan@154
|
194 elseif (function_mu)
|
ivan@154
|
195 y_cur = params.tau(i)*iy_cur + params.mu(i)*der.*setder;
|
ivan@154
|
196 else y_cur = params.tau(i)*iy_cur + params.mu*der.*setder;
|
ivan@154
|
197 end;
|
ivan@154
|
198
|
ivan@154
|
199 Y_i = ne(y_cur,0);
|
ivan@154
|
200 setYi = zeros(N,1);
|
ivan@154
|
201 setYi(Y_i) = 1;
|
ivan@154
|
202 end;
|
ivan@154
|
203
|
ivan@154
|
204 % Hard-threshold b and compute X_{i+1}
|
ivan@154
|
205 set_Xi_Yi = setXi + setYi;
|
ivan@154
|
206 ind_Xi_Yi = find(set_Xi_Yi > 0);
|
ivan@154
|
207 z = x_cur(ind_Xi_Yi) + y_cur(ind_Xi_Yi);
|
ivan@154
|
208 [~, ind_z] = sort(abs(z), 'descend');
|
ivan@154
|
209 x_cur = zeros(N,1);
|
ivan@154
|
210 x_cur(ind_Xi_Yi(ind_z(1:K))) = z(ind_z(1:K));
|
ivan@154
|
211 X_i = ind_Xi_Yi(ind_z(1:K));
|
ivan@154
|
212
|
ivan@154
|
213 setXi = zeros(N,1);
|
ivan@154
|
214 setXi(X_i) = 1;
|
ivan@154
|
215
|
ivan@154
|
216 if (params.gradientDescentx == 1)
|
ivan@154
|
217 % Calculate gradient of estimated vector x_cur
|
ivan@154
|
218 Phi_x_cur = Phi(:,X_i)*x_cur(X_i);
|
ivan@154
|
219 res = y - Phi_x_cur;
|
ivan@154
|
220 der = Phi_t*res;
|
ivan@154
|
221
|
ivan@154
|
222 ider = der(X_i);
|
ivan@154
|
223
|
ivan@154
|
224 if (adaptive_mu)
|
ivan@154
|
225 Pder = Phi(:,X_i)*ider;
|
ivan@154
|
226 mu_bar = ider'*ider/(Pder'*Pder);
|
ivan@154
|
227 x_cur(X_i) = x_cur(X_i) + mu_bar*ider;
|
ivan@154
|
228 elseif (function_mu)
|
ivan@154
|
229 x_cur = x_cur(X_i) + params.mu(i)*ider;
|
ivan@154
|
230 else x_cur = x_cur(X_i) + params.mu*ider;
|
ivan@154
|
231 end;
|
ivan@154
|
232 elseif (params.solveNewtonx == 1)
|
ivan@154
|
233 % Similar to HTP
|
ivan@154
|
234 if (params.useCG == 1)
|
ivan@154
|
235 [v, ~, ~] = cgsolve(cg_A(X_i, X_i), cg_b(X_i), params.cg_tol, params.cg_maxiter, cg_verbose);
|
ivan@154
|
236 else
|
ivan@154
|
237 v = cg_A(X_i,X_i)\cg_b(X_i);
|
ivan@154
|
238 end;
|
ivan@154
|
239 x_cur(X_i) = v;
|
ivan@154
|
240 end;
|
ivan@154
|
241
|
ivan@154
|
242 % Test stopping criterion
|
ivan@154
|
243 if (i > 1) && (norm(x_cur - x_prev) < params.tol*norm(x_cur))
|
ivan@154
|
244 break;
|
ivan@154
|
245 end;
|
ivan@154
|
246 i = i + 1;
|
ivan@154
|
247 end;
|
ivan@154
|
248
|
ivan@154
|
249 x_hat = x_cur;
|
ivan@154
|
250 numiter = i;
|
ivan@154
|
251
|
ivan@154
|
252 if (i > params.ALPSiters)
|
ivan@154
|
253 x_path = x_path(:,1:numiter-1);
|
ivan@154
|
254 else
|
ivan@154
|
255 x_path = x_path(:,1:numiter);
|
ivan@154
|
256 end;
|