idamnjanovic@1: function [A]=SMALL_chol(Dict,X, m, maxNumCoef, errorGoal, varargin) 
ivan@128: %% Implementation of OMP with Cholesky factorisation
ivan@128: % Sparse coding of a group of signals based on a given 
ivan@128: % dictionary and specified number of atoms to use. 
ivan@128: % input arguments:  Dict - the dictionary
ivan@128: %                   X - the signals to represent
ivan@128: %                   m - number of atoms in Dictionary
ivan@128: %                   errorGoal - the maximal allowed representation error for
ivan@128: %                   each signal.
ivan@128: %
ivan@128: % optional:         if Dict is function handle then Transpose Dictionary
ivan@128: %                   handle needs to be specified.
ivan@128: %
ivan@128: % output arguments: A - sparse coefficient matrix.
ivan@128: %
ivan@128: 
idamnjanovic@22: %
idamnjanovic@22: %   Centre for Digital Music, Queen Mary, University of London.
idamnjanovic@22: %   This file copyright 2009 Ivan Damnjanovic.
idamnjanovic@22: %
idamnjanovic@22: %   This program is free software; you can redistribute it and/or
idamnjanovic@22: %   modify it under the terms of the GNU General Public License as
idamnjanovic@22: %   published by the Free Software Foundation; either version 2 of the
idamnjanovic@22: %   License, or (at your option) any later version.  See the file
idamnjanovic@22: %   COPYING included with this distribution for more information.
idamnjanovic@22: %   
idamnjanovic@1: %%
idamnjanovic@1: % This Dictionary check is based on Thomas Blumensath work in sparsify 0_4 greedy solvers 
idamnjanovic@1: explicitD=0;
idamnjanovic@1: if     isa(Dict,'float')      
idamnjanovic@1:             D =@(z) Dict*z;  
idamnjanovic@1:             Dt =@(z) Dict'*z;
idamnjanovic@63:             explicitD=0;
idamnjanovic@1: elseif isobject(Dict) 
idamnjanovic@1:             D =@(z) Dict*z;  
idamnjanovic@1:             Dt =@(z) Dict'*z;
idamnjanovic@1: elseif isa(Dict,'function_handle') 
idamnjanovic@1:     try
idamnjanovic@1:         DictT=varargin{1};
idamnjanovic@1:         if isa(DictT,'function_handle'); 
idamnjanovic@74:                 D=Dict; 
idamnjanovic@74:                 Dt=DictT;
idamnjanovic@1:         else
idamnjanovic@1:                 error('If Dictionary is a function handle,Transpose Dictionary also needs to be a function handle. ');
idamnjanovic@1:         end
idamnjanovic@1:     catch
idamnjanovic@1:         error('If Dictionary is a function handle, Transpose Dictionary  needs to be specified. Exiting.');
idamnjanovic@1:     end
idamnjanovic@1: else
idamnjanovic@1:     error('Dictionary is of unsupported type. Use explicit matrix, function_handle or object. Exiting.');
idamnjanovic@1: end
idamnjanovic@1: %%
idamnjanovic@1: [n,P]=size(X);
idamnjanovic@1: 
idamnjanovic@1: 
idamnjanovic@1: 
idamnjanovic@1: global opts opts_tr machPrec
idamnjanovic@1: opts.UT = true; 
idamnjanovic@1: opts_tr.UT = true; opts_tr.TRANSA = true;
idamnjanovic@1: machPrec = 1e-5;
idamnjanovic@1: 
idamnjanovic@1: A = sparse(m,size(X,2));
idamnjanovic@1: for k=1:1:P,
idamnjanovic@1:     
idamnjanovic@1:     R_I = [];
idamnjanovic@1:     x=X(:,k);
idamnjanovic@1:     residual=x;
idamnjanovic@1: 	indx = [];
idamnjanovic@1: 	a = zeros(m,1);
idamnjanovic@1: 	currResNorm = norm(residual);
idamnjanovic@63:     errorGoal1=errorGoal*currResNorm;
idamnjanovic@1: 	j = 0;
idamnjanovic@1:     while currResNorm>errorGoal & j < maxNumCoef,
idamnjanovic@1: 		j = j+1;
idamnjanovic@1:         dir=Dt(residual);
idamnjanovic@1:         
idamnjanovic@1:         [tmp__, pos]=max(abs(dir));
idamnjanovic@1:         
idamnjanovic@1:         [R_I, flag] = updateChol(R_I, n, m, D, explicitD, indx, pos, Dt);
idamnjanovic@1:         
idamnjanovic@1:        
idamnjanovic@1:             indx(j)=pos;
idamnjanovic@1:             dx=zeros(m,1);
idamnjanovic@1:         
idamnjanovic@1:             z = linsolve(R_I,dir(indx),opts_tr);
idamnjanovic@1:         
idamnjanovic@1:             dx(indx) = linsolve(R_I,z,opts);
idamnjanovic@1:             a(indx) = a(indx) + dx(indx);
idamnjanovic@1:                 
idamnjanovic@1:             residual=x-D(a);
idamnjanovic@1:             currResNorm = norm(residual);
idamnjanovic@1:         
idamnjanovic@1:     
idamnjanovic@1:    end;
idamnjanovic@1:    if (~isempty(indx))
idamnjanovic@1:        A(indx,k)=a(indx);
idamnjanovic@1:    end
idamnjanovic@1: end;
idamnjanovic@1: return;
idamnjanovic@1: 
idamnjanovic@1: 
idamnjanovic@1: function [R, flag] = updateChol(R, n, N, A, explicitA, activeSet, newIndex, varargin)
idamnjanovic@1: 
idamnjanovic@1: % updateChol: Updates the Cholesky factor R of the matrix 
idamnjanovic@1: % A(:,activeSet)'*A(:,activeSet) by adding A(:,newIndex)
idamnjanovic@1: % If the candidate column is in the span of the existing 
idamnjanovic@1: % active set, R is not updated, and flag is set to 1.
idamnjanovic@1: 
idamnjanovic@1: global opts_tr machPrec
idamnjanovic@1: flag = 0;
idamnjanovic@1: 
idamnjanovic@1: if (explicitA)
idamnjanovic@1:     newVec = A(:,newIndex);
idamnjanovic@1: else
idamnjanovic@1:     At=varargin{1};
idamnjanovic@1:     e = zeros(N,1);
idamnjanovic@1:     e(newIndex) = 1;
idamnjanovic@1:     newVec = A(e);%feval(A,1,n,N,e,1:N,N); 
idamnjanovic@1: end
idamnjanovic@1: 
idamnjanovic@1: if isempty(activeSet),
idamnjanovic@1:     R = sqrt(sum(newVec.^2));
idamnjanovic@1: else
idamnjanovic@1:     if (explicitA)
idamnjanovic@1:         p = linsolve(R,A(:,activeSet)'*A(:,newIndex),opts_tr);
idamnjanovic@1:     else
idamnjanovic@1:         AnewVec = At(newVec);%feval(A,2,n,length(activeSet),newVec,activeSet,N);
idamnjanovic@1:         p = linsolve(R,AnewVec(activeSet),opts_tr);
idamnjanovic@1:     end
idamnjanovic@1:     q = sum(newVec.^2) - sum(p.^2);
idamnjanovic@1:     if (q <= machPrec) % Collinear vector
idamnjanovic@1:         flag = 1;
idamnjanovic@1:     else
idamnjanovic@1:         R = [R p; zeros(1, size(R,2)) sqrt(q)];
idamnjanovic@1:     end
idamnjanovic@1: end