function [Dictionary,output] = SGK(...
Data,... % an nXN matrix that contins N signals (Y), each of dimension n.
param)
% -----------------------------------------------------------------------
% Author: Sujit Kumar Sahoo, School of Electrical and Electronic
% Engineering, Nanyang Technological University, Singapore.
% =========================================================================
% A SEQUENTIAL GENERALIZATION OF K-MEANS (SGK) algorithm
% =========================================================================
% The SGK algorithm finds a dictionary for linear representation of
% signals. Given a set of signals, it searches for the best dictionary that
% can sparsely represent each signal. Detailed discussion on the algorithm
% can be found in "Dictionary Training for Sparse Representation as
% Generalization of K-means Clustering", written by S.K. Sahoo and A. Makur
% and to be appeared in the IEEE Signal Processing Letter, 2013.
% =========================================================================
% INPUT ARGUMENTS:
% Data an nXN matrix that contins N signals (Y), each of dimension n.
% param structure that includes all required
% parameters for the SGK execution.
% Required fields are:
% K, ... the number of dictionary elements to train
% numIteration,... number of iterations to perform.
% errorFlag... if =0, a fix number of coefficients is
% used for representation of each signal. If so, param.m must be
% specified as the number of representing atom. if =1, arbitrary number
% of atoms represent each signal, until a specific representation error
% is reached. If so, param.errorGoal must be specified as the allowed
% error.
% preserveDCAtom... if =1 then the first atom in the dictionary
% is set to be constant, and does not ever change. This
% might be useful for working with natural
% images (in this case, only param.K-1
% atoms are trained).
% (optional, see errorFlag) L,... % maximum coefficients to use in OMP coefficient calculations.
% (optional, see errorFlag) errorGoal, ... % allowed representation error in representing each signal.
% InitializationMethod,... mehtod to initialize the dictionary, can
% be one of the following arguments:
% * 'DataElements' (initialization by the signals themselves), or:
% * 'GivenMatrix' (initialization by a given matrix param.initialDictionary).
% (optional, see InitializationMethod) initialDictionary,... % if the initialization method
% is 'GivenMatrix', this is the matrix that will be used.
% (optional) TrueDictionary, ... % if specified, in each
% iteration the difference between this dictionary and the trained one
% is measured and displayed.
% displayProgress, ... if =1 progress information is displyed. If param.errorFlag==0,
% the average repersentation error (RMSE) is displayed, while if
% param.errorFlag==1, the average number of required coefficients for
% representation of each signal is displayed.
% =========================================================================
% OUTPUT ARGUMENTS:
% Dictionary The extracted dictionary of size nX(param.K).
% output Struct that contains information about the current run. It may include the following fields:
% CoefMatrix The final coefficients matrix (it should hold that Data equals approximately Dictionary*output.CoefMatrix.
% ratio If the true dictionary was defined (in
% synthetic experiments), this parameter holds a vector of length
% param.numIteration that includes the detection ratios in each
% iteration).
% totalerr The total representation error after each
% iteration (defined only if
% param.displayProgress=1 and
% param.errorFlag = 0)
% numCoef A vector of length param.numIteration that
% include the average number of coefficients required for representation
% of each signal (in each iteration) (defined only if
% param.displayProgress=1 and
% param.errorFlag = 1)
% time Average time taken for one iteration of
% the dictionary update algorithm.
% =========================================================================
if (~isfield(param,'displayProgress'))
param.displayProgress = 0;
end
totalerr = 99999*ones(param.numIteration,1);
if (isfield(param,'errorFlag')==0)
param.errorFlag = 0;
end
if (isfield(param,'TrueDictionary'))
displayErrorWithTrueDictionary = 1;
ErrorBetweenDictionaries = zeros(param.numIteration+1,1);
ratio = zeros(param.numIteration+1,1);
else
displayErrorWithTrueDictionary = 0;
ratio = 0;
end
if (param.preserveDCAtom>0)
FixedDictionaryElement = ones(size(Data,1),1)/sqrt(size(Data,1));
else
FixedDictionaryElement = [];
end
% coefficient calculation method is OMP with fixed number of coefficients
Dictionary = zeros(size(Data,1),param.K);
if (size(Data,2) < param.K)
disp('Size of data is smaller than the dictionary size. Trivial solution...');
Dictionary = Data(:,1:size(Data,2));
return;
elseif (strcmp(param.InitializationMethod,'DataElements'))
Dictionary(:,1:param.K-param.preserveDCAtom) = Data(:,1:param.K-param.preserveDCAtom);
elseif (strcmp(param.InitializationMethod,'GivenMatrix'))
Dictionary = param.initialDictionary;
end
% reduce the components in Dictionary that are spanned by the fixed
% elements
if (param.preserveDCAtom)
tmpMat = FixedDictionaryElement \ Dictionary;
Dictionary = Dictionary - FixedDictionaryElement*tmpMat;
Dictionary = Dictionary(:,sum(Dictionary.^2)>10^-2);
end
%normalize the dictionary.
Dictionary = Dictionary*diag(1./sqrt(sum(Dictionary.*Dictionary)));
Dictionary = Dictionary.*repmat(sign(Dictionary(1,:)),size(Dictionary,1),1); % multiply in the sign of the first element.
[n,N] = size(Data);
% the SGK Approximation algorithm starts here.
time=0;
for iterNum = 1:param.numIteration
D = [FixedDictionaryElement,Dictionary];
aNorms = sqrt(sum(D.^2));
D = D./repmat(aNorms,[n,1]);
G = D'*D; epsilon = sqrt(size(D,1))*param.errorGoal;
% sparse coding [You can replace below with your own OMP algorithm]
%----------------------------------------------------------------------
if (param.errorFlag==0)
%%% Faster computation of OMP
CoefMatrix = OMPa(D,Data,param.m,0,0);
else
%%% Faster computation of OMP
param.m = 1;
CoefMatrix = OMPa(D,Data,param.m,inf,param.errorGoal^2);
end
%----------------------------------------------------------------------
CoefMatrix = CoefMatrix./repmat(aNorms',[1,N]);
%***********************************************************************
tic;
replacedVectorCounter = 0;
rPerm = randperm(size(Dictionary,2));
% rPerm = 1:size(Dictionary,2);
for j = rPerm
[better
- 1
- 2
前往页