ivan@137: function gamma = omp2(varargin)
ivan@137: %OMP2 Error-constrained Orthogonal Matching Pursuit.
ivan@137: %  GAMMA = OMP2(D,X,G,EPSILON) solves the optimization problem
ivan@137: %
ivan@137: %       min  |GAMMA|_0     s.t.  |X - D*GAMMA|_2 <= EPSILON
ivan@137: %      gamma
ivan@137: %
ivan@137: %  for each of the signals in X, using Batch Orthogonal Matching Pursuit.
ivan@137: %  Here, D is a dictionary with normalized columns, X is a matrix
ivan@137: %  containing column signals, EPSILON is the error target for each signal,
ivan@137: %  and G is the Gramm matrix D'*D. The output GAMMA is a matrix containing
ivan@137: %  the sparse representations as its columns. 
ivan@137: %
ivan@137: %  GAMMA = OMP2(D,X,[],EPSILON) performs the same operation, but without
ivan@137: %  the matrix G, using OMP-Cholesky. This call produces the same output as
ivan@137: %  Batch-OMP, but is significantly slower. Using this syntax is only
ivan@137: %  recommended when available memory is too small to store G.
ivan@137: %
ivan@137: %  GAMMA = OMP2(DtX,XtX,G,EPSILON) is the fastest implementation of OMP2,
ivan@137: %  but also requires the most memory. Here, DtX stores the projections
ivan@137: %  D'*X, and XtX is a row vector containing the squared norms of the
ivan@137: %  signals, sum(X.*X). In this case Batch-OMP is used, but without having
ivan@137: %  to compute D'*X and XtX in advance, which slightly improves runtime.
ivan@137: %  Note that in general, the call
ivan@137: %
ivan@137: %    GAMMA = OMP2(D'*X, sum(X.*X), G, EPSILON);
ivan@137: %
ivan@137: %  will be faster than the call
ivan@137: %
ivan@137: %    GAMMA = OMP2(D,X,G,EPSILON);
ivan@137: %
ivan@137: %  due to optimized matrix multiplications in Matlab. However, when the
ivan@137: %  entire matrix D'*X cannot be stored in memory, one of the other two
ivan@137: %  versions can be used. Both compute D'*X for just one signal at a time,
ivan@137: %  and thus require much less memory.
ivan@137: %
ivan@137: %  GAMMA = OMP2(...,PARAM1,VAL1,PARAM2,VAL2,...) specifies additional
ivan@137: %  parameters for OMP2. Available parameters are:
ivan@137: %
ivan@137: %    'gammamode' - Specifies the representation mode for GAMMA. Can be
ivan@137: %                  either 'full' or 'sparse', corresponding to a full or
ivan@137: %                  sparse matrix, respectively. By default, GAMMA is
ivan@137: %                  returned as a sparse matrix.
ivan@137: %    'maxatoms' -  Limits the number of atoms in the representation of each
ivan@137: %                  signal. If specified, the number of atoms in each
ivan@137: %                  representation does not exceed this number, even if the
ivan@137: %                  error target is not met. Specifying maxatoms<0 implies
ivan@137: %                  no limit (default).
ivan@137: %    'messages'  - Specifies whether progress messages should be displayed.
ivan@137: %                  When positive, this is the number of seconds between
ivan@137: %                  status prints. When negative, indicates that no messages
ivan@137: %                  should be displayed (this is the default).
ivan@137: %    'checkdict' - Specifies whether dictionary normalization should be
ivan@137: %                  verified. When set to 'on' (default) the dictionary
ivan@137: %                  atoms are verified to be of unit L2-norm. Setting this
ivan@137: %                  parameter to 'off' disables verification and accelerates
ivan@137: %                  function performance. Note that an unnormalized
ivan@137: %                  dictionary will produce invalid results.
ivan@137: %    'profile'   - Can be either 'on' or 'off'. When 'on', profiling
ivan@137: %                  information is displayed at the end of the funciton
ivan@137: %                  execution.
ivan@137: %
ivan@137: %
ivan@137: %  Summary of OMP2 versions:
ivan@137: %
ivan@137: %    version                 |   speed     |   memory
ivan@137: %  -------------------------------------------------------------
ivan@137: %   OMP2(DtX,XtX,G,EPSILON)  |  very fast  |  very large
ivan@137: %   OMP2(D,X,G,EPSILON)      |  fast       |  moderate
ivan@137: %   OMP2(D,X,[],EPSILON)     |  very slow  |  small
ivan@137: %  -------------------------------------------------------------
ivan@137: %
ivan@137: %
ivan@137: %  References:
ivan@137: %  [1] M. Elad, R. Rubinstein, and M. Zibulevsky, "Efficient Implementation
ivan@137: %      of the K-SVD Algorithm using Batch Orthogonal Matching Pursuit",
ivan@137: %      Technical Report - CS, Technion, April 2008.
ivan@137: %
ivan@137: %  See also OMP.
ivan@137: 
ivan@137: 
ivan@137: %  Ron Rubinstein
ivan@137: %  Computer Science Department
ivan@137: %  Technion, Haifa 32000 Israel
ivan@137: %  ronrubin@cs
ivan@137: %
ivan@137: %  April 2009
ivan@137: 
ivan@137: 
ivan@137: % default options
ivan@137: 
ivan@137: sparse_gamma = 1;
ivan@137: msgdelta = -1;
ivan@137: maxatoms = -1;
ivan@137: checkdict = 1;
ivan@137: profile = 0;
ivan@137: 
ivan@137: 
ivan@137: % determine number of parameters
ivan@137: 
ivan@137: paramnum = 1;
ivan@137: while (paramnum<=nargin && ~ischar(varargin{paramnum}))
ivan@137:   paramnum = paramnum+1;
ivan@137: end
ivan@137: paramnum = paramnum-1;
ivan@137: 
ivan@137: 
ivan@137: % parse options
ivan@137: 
ivan@137: for i = paramnum+1:2:length(varargin)
ivan@137:   paramname = varargin{i};
ivan@137:   paramval = varargin{i+1};
ivan@137: 
ivan@137:   switch lower(paramname)
ivan@137: 
ivan@137:     case 'gammamode'
ivan@137:       if (strcmpi(paramval,'sparse'))
ivan@137:         sparse_gamma = 1;
ivan@137:       elseif (strcmpi(paramval,'full'))
ivan@137:         sparse_gamma = 0;
ivan@137:       else
ivan@137:         error('Invalid GAMMA mode');
ivan@137:       end
ivan@137:       
ivan@137:     case 'maxatoms'
ivan@137:       maxatoms = paramval;
ivan@137:       
ivan@137:     case 'messages'
ivan@137:       msgdelta = paramval;
ivan@137: 
ivan@137:     case 'checkdict'
ivan@137:       if (strcmpi(paramval,'on'))
ivan@137:         checkdict = 1;
ivan@137:       elseif (strcmpi(paramval,'off'))
ivan@137:         checkdict = 0;
ivan@137:       else
ivan@137:         error('Invalid checkdict option');
ivan@137:       end
ivan@137: 
ivan@137:     case 'profile'
ivan@137:       if (strcmpi(paramval,'on'))
ivan@137:         profile = 1;
ivan@137:       elseif (strcmpi(paramval,'off'))
ivan@137:         profile = 0;
ivan@137:       else
ivan@137:         error('Invalid profile mode');
ivan@137:       end
ivan@137: 
ivan@137:     otherwise
ivan@137:       error(['Unknown option: ' paramname]);
ivan@137:   end
ivan@137:   
ivan@137: end
ivan@137: 
ivan@137: 
ivan@137: % determine call type
ivan@137: 
ivan@137: if (paramnum==4)
ivan@137:   
ivan@137:   n1 = size(varargin{1},1);
ivan@137:   n2 = size(varargin{2},1);
ivan@137:   n3 = size(varargin{3},1);
ivan@137:   
ivan@137:   if ( (n1>1 && n2==1) || (n1==1 && n2==1 && n3==1) )  %  DtX,XtX,G,EPSILON
ivan@137:     
ivan@137:     DtX = varargin{1};
ivan@137:     XtX = varargin{2};
ivan@137:     G = varargin{3};
ivan@137:     epsilon = varargin{4};
ivan@137:     D = [];
ivan@137:     X = [];
ivan@137:     
ivan@137:   else  % D,X,G,EPSILON
ivan@137:     
ivan@137:     D = varargin{1};
ivan@137:     X = varargin{2};
ivan@137:     G = varargin{3};
ivan@137:     epsilon = varargin{4};
ivan@137:     DtX = [];
ivan@137:     XtX = [];
ivan@137:     
ivan@137:   end
ivan@137:   
ivan@137: else
ivan@137:   error('Invalid number of parameters');
ivan@137: end
ivan@137: 
ivan@137: 
ivan@137: % verify dictionary normalization
ivan@137: 
ivan@137: if (checkdict)
ivan@137:   if (isempty(G))
ivan@137:     atomnorms = sum(D.*D);
ivan@137:   else
ivan@137:     atomnorms = diag(G);
ivan@137:   end
ivan@137:   if (any(abs(atomnorms-1) > 1e-2))
ivan@137:     error('Dictionary columns must be normalized to unit length');
ivan@137:   end
ivan@137: end
ivan@137: 
ivan@137: 
ivan@137: % omp
ivan@137: 
ivan@137: gamma = omp2mexGabor(D,X,DtX,XtX,G,epsilon,sparse_gamma,msgdelta,maxatoms,profile);