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