vsab2.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

function w=VSab2(x,w)
% VSab2     Vector Selection algorithm that Always returns a Better (or the same) w
% This was made to ensure that the weights selected are at least as good as the 
% previous selecte weights. 
% Essentially this use the same algorithm as VSmp2
%
% w=VSab2(x,w);
% w=VSab2(x,S);    
%----------------------------------------------------------------------
% arguments:
%   w    - the weights where at most S of the elements are non-zero, size Kx1
%          note that w is both input and output argument
%   x    - the vector to approximate, size Nx1
%   S    - if second argument is not the weight vector
%          it is number of vectors to select or non-zero weights 
%          (do not mix with S in FindW function, which has a different meaning)
%          S (which is scalar) may be used instead of w as second input argument,
%          then w is found using VSmp2 or VSfomp2
%-----------------------------------------------------------------------------------
% global variables:
%   F    - the dictionary (or F matrix), size NxK
%          the vectors of F must be normalized, F(:,k)'*F(:,k)==1 for all k
%          normalization may be done by: F=F*diag(1./sqrt(diag(F'*F)));
%   FF   - the inner products F'*F, size KxK
%-----------------------------------------------------------------------------------
% Note that we here use a simple complete F matrix, and do not make any assumtions
% regarding structure of F.


%----------------------------------------------------------------------
% Copyright (c) 1999-2001.  Karl Skretting.  All rights reserved.
% Hogskolen in Stavanger (Stavanger University), Signal Processing Group
% Mail:  [e-mail]   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:
% Ver. 1.0  02.10.2001  KS: function made based on VSmp2
% Ver. 1.1  25.11.2002  KS: moved from ..\Frames to ..\FrameTools
%----------------------------------------------------------------------

global F FF

Mfile='VSab2';
Display=0;
[N,K]=size(F);

if nargin<2
   error([Mfile,': wrong number of arguments, see help.']);
end
w=w(:);
x=x(:);
if (length(w)~=K)       % size of w is not what it is expected to be, 
   S=floor(w(1));       % second argument is probably S
   w=zeros(K,1);        
   % The special case when second argument is S, we must find an initial w 
   if (S>=N); S=(N-1); end;
   if (S<0); S=0; end;
   it=2*S;
   wI=[];
   wIinv=1:K;
   r=x;
   k=0;
   while k<it
      c=(r'*F);                % the inner products
      if length(wI)==S
         c(wIinv)=0;           % we do not select any new vectors
      end
      [temp,i]=max(abs(c));i=i(1);
      if w(i)
         w(i)=w(i)+c(i);
      else
         w(i)=c(i);
         wI=find(w);
         wIinv=find(w==0);
      end   
      r=r-c(i)*F(:,i);
      k=k+1;                   % make sure the loop do not run forever
   end
end

%  now we have w
I=find(w);
S=length(I);
w0=w;
r0=x-F*w0;
rr0=r0'*r0;
if Display
   disp([Mfile,': Norm squared of initial error is ',num2str(rr0)]);
end
rr=rr0*2;
xx=x'*x;
SetSomeToZero=1;

if ((S>1) & (S<N) & (rr0>(1e-10*xx)))
   if (rand(1)<0.01) & (K<200) & (N<100)
      % we let it be a little chance that we use the largest coefficients 
      % after Gaussian eliminasion or linear programming
      if (rand(1)<0.9) 
         w3=F\x;                          % Gaussian elimination
      else
         f=ones(2*K,1);
         A=[F,-F];
         LB=zeros(2*K,1);
         UB=ones(2*K,1)*5000;             % this should be large enough
         w1=linprog(f,A,x,A,x,LB,UB);     % minimize 1-norm applying linear programming
         w3=w1(1:K)-w1((1+K):(2*K));      
      end
      [ws,II]=sort(-abs(w3));
      w=zeros(K,1);
      w(II(1:S))=w3(II(1:S));
      I=find(w);
   else   % the most likely case
      if rand(1)<0.1   
         if rand(1)<0.98        % we may also try VSfomp2
            w=VSfomp2(x,S);
         else                  % or occationally partial search
            P=ones(1,S);
            P(1)=ceil(rand(1)^2*13);
            if (S>2); P(2)=ceil(rand(1)^2*11); end;
            if (S>5); P(4)=ceil(rand(1)^2*7); end;
            if (S>6); P(5)=ceil(rand(1)*2); end;
            w=VSps(F,x,S,P);
            SetSomeToZero=0;
         end   
         I=find(w);
      elseif rand(1)<0.4
          it=floor(S*(2+6*rand(1)));
          w=VSmp2(x,S,it);
      end
   end
end
if SetSomeToZero
   for jj=1:1              % do this some times, or only once             
      wa=abs(w(I));
      [ws,II]=sort(wa);
      c=S:(-1):1;
      for j=1:ceil(S*(rand(1)*0.6+0.02))     % set some of the weights to zero
         t=rand(1)*S*(S+1)/2;
         i=1;k=0;
         while 1
            k=k+c(i);
            if (k>t); break; end;
            i=i+1;
         end
         w(I(II(i)))=0;
      end
      %
      it=S;          % number of iterations to do
      wI=find(w);    % the weights already found
      wIinv=find(w==0);
      r=x-F*w;
      k=0;
      while k<it
         c=(r'*F);                % the inner products
         if length(wI)==S
            c(wIinv)=0;           % we do not select any new vectors
         end
         [temp,i]=max(abs(c));i=i(1);
         if w(i)
            w(i)=w(i)+c(i);
         else
            w(i)=c(i);
            wI=find(w);
            wIinv=find(w==0);
         end   
         r=r-c(i)*F(:,i);
         k=k+1;                   % make sure the loop do not run forever
      end
      if (rand(1)<0.35)
         % we let it be a little chance that we project onto these vectors 
         Fi=F(:,wI);    % Fi is here only the selected vectors
         FiFi=Fi'*Fi;
         if rcond(FiFi)>1e-4
            w(wI)=(FiFi\(Fi'))*x;
            r=x-F*w;
         end
      end
      rr=r'*r;
      if Display
         disp([Mfile,': Norm squared of error ',int2str(jj),' is ',num2str(rr)]);
      end
      if (rr<rr0); break; end;
   end
end
% this check is the clue for this function, 
% we only use the new weights if they are better
I=find(w);
if ((rr0<rr) | (length(I)~=S))
   if Display
      disp([Mfile,': no improvement for this weight vector found.']);
   end
   w=w0;
end

return