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