1function [s, err_mse, iter_time]=greed_omp_cg(x,A,m,varargin) 2%-*- texinfo -*- 3%@deftypefn {Function} greed_omp_cg 4%@verbatim 5% greed_omp_cg: Orthogonal Matching Pursuit algorithm based on full 6% conjugate gradient solver 7%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 8% Usage 9% [s, err_mse, iter_time]=greed_omp_cg(x,P,m,'option_name','option_value') 10%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 11%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 12% Input 13% Mandatory: 14% x Observation vector to be decomposed 15% P Either: 16% 1) An nxm matrix (n must be dimension of x) 17% 2) A function handle (type "help function_format" 18% for more information) 19% Also requires specification of P_trans option. 20% 3) An object handle (type "help object_format" for 21% more information) 22% m length of s 23% 24% Possible additional options: 25% (specify as many as you want using 'option_name','option_value' pairs) 26% See below for explanation of options: 27%__________________________________________________________________________ 28% option_name AVAILABLE OPTION_VALUES default 29%-------------------------------------------------------------------------- 30% stopCrit M, CORR, MSE, MSE_CHANGE M 31% stopTol NUMBER (SEE BELOW) n/4 32% cgTol NUMBER (SEE BELOW) 1e-6 33% P_trans FUNCTION_HANDLE (SEE BELOW) 34% maxIter POSITIVE INTEGER (SEE BELOW) n 35% verbose TRUE, FALSE false 36% start_val VECTOR OF LENGTH M zeros 37% 38% Available stopping criteria : 39% M - Extracts exactly M = stopTol elements. 40% corr - Stops when maximum correlation between 41% residual and atoms is below stopTol value. 42% mse - Stops when mean squared error of residual 43% is below stopTol value. 44% mse_change - Stops when the change in the mean squared 45% error falls below stopTol value. 46% 47% stopTol: Value for stopping criterion. 48% 49% cgTol: Value for conjugate gradient stoping criterion. 50% 51% P_trans: If P is a function handle, then P_trans has to be specified and 52% must be a function handle. 53% 54% maxIter: Maximum of allowed iterations. 55% 56% verbose: Logical value to allow algorithm progress to be displayed. 57% 58% start_val: Allows algorithms to start from partial solution. 59% 60%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 61% Outputs 62% s Solution vector 63% err_mse Vector containing mse of approximation error for each 64% iteration 65% iter_time Vector containing times for each iteration 66% 67%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 68% Description 69% greed_omp_cg performs a greedy signal decomposition. 70% In each iteration a new element is selected depending on the inner 71% product between the current residual and columns in P. 72% The non-zero elements of s are approximated by orthogonally projecting 73% x onto the selected elements in each iteration. 74% The inverse problem required for orthogonal matching pursuit is solved 75% using a conjugate gradient algorithm in each iteration. This is slow but 76% does not require any additional storage and can therefore be used for 77% large scale problems. 78% 79% See Also 80% greed_qr, greed_omp_chol, greed_omp_cgp, greed_omp_pinv, 81% greed_omp_linsolve, greed_gp, greed_nomp 82% 83% Copyright (c) 2007 Thomas Blumensath 84% 85% The University of Edinburgh 86% Email: thomas.blumensath@ed.ac.uk 87% Comments and bug reports welcome 88% 89% This file is part of sparsity Version 0.1 90% Created: April 2007 91% 92% Part of this toolbox was developed with the support of EPSRC Grant 93% D000246/1 94% 95% Please read COPYRIGHT.m for terms and conditions. 96%@end verbatim 97%@strong{Url}: @url{http://ltfat.github.io/doc/thirdparty/sparsify/private/greed_omp_cg.html} 98%@end deftypefn 99 100% Copyright (C) 2005-2016 Peter L. Soendergaard <peter@sonderport.dk>. 101% This file is part of LTFAT version 2.3.1 102% 103% This program is free software: you can redistribute it and/or modify 104% it under the terms of the GNU General Public License as published by 105% the Free Software Foundation, either version 3 of the License, or 106% (at your option) any later version. 107% 108% This program is distributed in the hope that it will be useful, 109% but WITHOUT ANY WARRANTY; without even the implied warranty of 110% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 111% GNU General Public License for more details. 112% 113% You should have received a copy of the GNU General Public License 114% along with this program. If not, see <http://www.gnu.org/licenses/>. 115 116%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 117% Default values and initialisation 118%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 119 120 121[n1 n2]=size(x); 122if n2 == 1 123 n=n1; 124elseif n1 == 1 125 x=x'; 126 n=n2; 127else 128 display('x must be a vector.'); 129 return 130end 131 132sigsize = x'*x/n; 133initial_given=0; 134err_mse = []; 135iter_time = []; 136STOPCRIT = 'M'; 137STOPTOL = ceil(n/4); 138MAXITER = n; 139verbose = false; 140s_initial = zeros(m,1); 141cg_tol = 1e-6; 142vectnfact = ones(m,1); 143 144if verbose 145 display('Initialising...') 146end 147 148%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 149% Output variables 150%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 151 152switch nargout 153 case 3 154 comp_err=true; 155 comp_time=true; 156 case 2 157 comp_err=true; 158 comp_time=false; 159 case 1 160 comp_err=false; 161 comp_time=false; 162 case 0 163 error('Please assign output variable.') 164 otherwise 165 error('Too many output arguments specified') 166end 167 168%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 169% Look through options 170%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 171 172% Put option into nice format 173Options={}; 174OS=nargin-3; 175c=1; 176for i=1:OS 177 if isa(varargin{i},'cell') 178 CellSize=length(varargin{i}); 179 ThisCell=varargin{i}; 180 for j=1:CellSize 181 Options{c}=ThisCell{j}; 182 c=c+1; 183 end 184 else 185 Options{c}=varargin{i}; 186 c=c+1; 187 end 188end 189OS=length(Options); 190if rem(OS,2) 191 error('Something is wrong with argument name and argument value pairs.') 192end 193 194for i=1:2:OS 195 switch Options{i} 196 case {'stopCrit'} 197 if (strmatch(Options{i+1},{'M'; 'corr'; 'mse'; 'mse_change'},'exact')); 198 STOPCRIT = Options{i+1}; 199 else error('stopCrit must be char string [M, corr, mse, mse_change]. Exiting.'); end 200 case {'stopTol'} 201 if isa(Options{i+1},'numeric') ; STOPTOL = Options{i+1}; 202 else error('stopTol must be number. Exiting.'); end 203 case {'cgTol'} 204 if isa(Options{i+1},'numeric') ; cg_tol = Options{i+1}; 205 else error('cgTol must be number. Exiting.'); end 206 case {'P_trans'} 207 if isa(Options{i+1},'function_handle'); Pt = Options{i+1}; 208 else error('P_trans must be function _handle. Exiting.'); end 209 case {'maxIter'} 210 if isa(Options{i+1},'numeric'); MAXITER = Options{i+1}; 211 else error('maxIter must be a number. Exiting.'); end 212 case {'verbose'} 213 if isa(Options{i+1},'logical'); verbose = Options{i+1}; 214 else error('verbose must be a logical. Exiting.'); end 215 case {'vecNormFac'} 216 if isa(Options{i+1},'numeric')& length(Options{i+1}) == m , vectnfact = Options{i+1}; 217 else error('verbose must be a logical. Exiting.'); end 218 case {'start_val'} 219 if isa(Options{i+1},'numeric') & length(Options{i+1}) == m ; 220 s_initial = Options{i+1}; 221 initial_given=1; 222 else error('start_val must be a vector of length m. Exiting.'); end 223 otherwise 224 error('Unrecognised option. Exiting.') 225 end 226end 227 228 229 230if strcmp(STOPCRIT,'M') 231 maxM=STOPTOL; 232else 233 maxM=MAXITER; 234end 235 236if nargout >=2 237 err_mse = zeros(maxM,1); 238end 239if nargout ==3 240 iter_time = zeros(maxM,1); 241end 242 243 244 245%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 246% Make P and Pt functions 247%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 248 249if isa(A,'float') P =@(z) A*z; Pt =@(z) A'*z; 250elseif isobject(A) P =@(z) A*z; Pt =@(z) A'*z; 251elseif isa(A,'function_handle') 252 try 253 if isa(Pt,'function_handle'); P=A; 254 else error('If P is a function handle, Pt also needs to be a function handle. Exiting.'); end 255 catch error('If P is a function handle, Pt needs to be specified. Exiting.'); end 256else error('P is of unsupported type. Use matrix, function_handle or object. Exiting.'); end 257 258%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 259% Do we start from zero or not? 260%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 261 262if initial_given ==1; 263 IN = find(s_initial); 264 Residual = x-P(s_initial); 265 s = s_initial; 266 oldERR = Residual'*Residual/n; 267 268else 269 IN = []; 270 Residual = x; 271 s = s_initial; 272 sigsize = x'*x/n; 273 oldERR = sigsize; 274end 275 276%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 277% Random Check to see if dictionary is normalised 278%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 279 280% mask=zeros(m,1); 281% mask(ceil(rand*m))=1; 282% nP=norm(P(mask)); 283% if abs(1-nP)>1e-3; 284% display('Dictionary appears not to have unit norm columns.') 285% end 286 287 288%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 289% Main algorithm 290%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 291if verbose 292 display('Main iterations...') 293end 294tic 295t=0; 296DR=Pt(Residual).*vectnfact; 297done = 0; 298iter=1; 299while ~done 300 DR(IN)=0; 301 [v I]=max(abs(DR)); 302 IN=[IN I]; 303 % abs(DR(IN)) 304 % pause 305 [s Residual]=SubsetCG(x,s,P,Pt,IN,cg_tol,verbose); 306 DR=Pt(Residual).*vectnfact; 307 308 309 ERR=Residual'*Residual/n; 310 if comp_err 311 err_mse(iter)=ERR; 312 end 313 314 if comp_time 315 iter_time(iter)=toc; 316 end 317 318%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 319% Are we done yet? 320%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 321 if strcmp(STOPCRIT,'M') 322 if iter >= STOPTOL 323 done =1; 324 elseif verbose && toc-t>10 325 display(sprintf('Iteration %i. --- %i iterations to go',iter ,STOPTOL-iter)) 326 t=toc; 327 end 328 elseif strcmp(STOPCRIT,'mse') 329 if comp_err 330 if err_mse(iter)<STOPTOL; 331 done = 1; 332 elseif verbose && toc-t>10 333 display(sprintf('Iteration %i. --- %i mse',iter ,err_mse(iter))) 334 t=toc; 335 end 336 else 337 if ERR<STOPTOL; 338 done = 1; 339 elseif verbose && toc-t>10 340 display(sprintf('Iteration %i. --- %i mse',iter ,ERR)) 341 t=toc; 342 end 343 end 344 elseif strcmp(STOPCRIT,'mse_change') && iter >=2 345 if comp_err && iter >=2 346 if ((err_mse(iter-1)-err_mse(iter))/sigsize <STOPTOL); 347 done = 1; 348 elseif verbose && toc-t>10 349 display(sprintf('Iteration %i. --- %i mse change',iter ,(err_mse(iter-1)-err_mse(iter))/sigsize )) 350 t=toc; 351 end 352 else 353 if ((oldERR - ERR)/sigsize < STOPTOL); 354 done = 1; 355 elseif verbose && toc-t>10 356 display(sprintf('Iteration %i. --- %i mse change',iter ,(oldERR - ERR)/sigsize)) 357 t=toc; 358 end 359 end 360 elseif strcmp(STOPCRIT,'corr') 361 if max(abs(DR)) < STOPTOL; 362 done = 1; 363 elseif verbose && toc-t>10 364 display(sprintf('Iteration %i. --- %i corr',iter ,max(abs(DR)))) 365 t=toc; 366 end 367 end 368 369 % Also stop if residual gets too small or maxIter reached 370 if comp_err 371 if err_mse(iter)<1e-16 372 display('Stopping. Exact signal representation found!') 373 done=1; 374 end 375 else 376 377 378 if iter>1 379 if ERR<1e-16 380 display('Stopping. Exact signal representation found!') 381 done=1; 382 end 383 end 384 end 385 386 if iter >= MAXITER 387 display('Stopping. Maximum number of iterations reached!') 388 done = 1; 389 end 390 391%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 392% If not done, take another round 393%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 394 395 if ~done 396 iter=iter+1; 397 oldERR=ERR; 398 end 399end 400 401%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 402% Only return as many elements as iterations 403%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 404 405if nargout >=2 406 err_mse = err_mse(1:iter); 407end 408if nargout ==3 409 iter_time = iter_time(1:iter); 410end 411 412if verbose 413 display('Done') 414end 415 416% Change history 417% 418% 8 of Februray: Algo does no longer stop if dictionary is not normaliesd. 419