How to solve problem of out memory

2 views (last 30 days)
Raoult TEUKAM DABOU
Raoult TEUKAM DABOU on 22 Aug 2018
Edited: OCDER on 22 Aug 2018
function varargout = sparseapprox(X, D, met, varargin)
% sparseapprox Returns coefficients in a sparse approximation of X.
% Several methods for sparse approximation may be used,
% some implemented in this m-file and others depend on external parts.
% For the methods (3) and (4) below, the corresponding packages
% should be installed on your computer and be available from Matlab.
%
% The coefficients or weights, W, are usually (but not always) sparse,
% i.e. number of non-zero coefficients below a limit,
% some methods use the 1-norm but these are not much tested.
%
% The reconstrucion or approximation of X is (D*W).
% The approximation error is R = X - D*W;
% The Signal to Noise Ratio is snr = 10*log10(var(X(:))/var(R(:)));
%
% Use of function:
% ----------------
% W = SPARSEAPPROX(X, D, met, 'option',value, ...)
% W is coefficient matrix, size KxL
% X is data, a matrix of size NxL, (a column vector when L=1)
% D is the dictionary, size NxK
% met is a char string for the different methods, see below
% Additional arguments may be given as pairs: 'option',value, ...
% (or options may be given in one (or more) struct or cell-array)
% These options may be used for the different methods.
%
% [W, res] = SPARSEAPPROX(X, D, met, 'option',value, ...)
% res is a struct with additional results
%
% The alternative methods are:
% ----------------------------
% (1) Methods that use Matlab standard functions: 'pinv', '\', 'linprog'
% The representation is now exact and usually not sparse unless
% thresholding is done (see below).
% 'MOF', 'MethodOfFrames', or 'pinv'
% 'BackSlash' or '\'
% 'BP' or 'BasisPursuit' or 'linprog'
% (2) Methods implemented in this m-file
% 'FOCUSS' a best basis variant. Use options 'nIt' and 'pFOCUSS'.
% When option 'lambdaFOCUSS' is given Regularized FOCUSS is used.
% For the four methods ('pinv', 'BackSlash', 'linprog' and 'FOCUSS')
% thresholding is done when 'tnz', 'tre' or 'tae' is given as option.
% If 'doOP' the coefficients are set so that D*W(:,i) is an orthogonal
% projection onto the space spanned by used columns of D.
% 'GMP', a Global variant of Matching Pursuit. NOTE that option 'tnz'
% should be given as the total number of non-zeros in W
% 'OMP', Orthogonal Matching Pursuit
% 'ORMP', Order Recursive Matching Pursuit
% (3) Methods implemented in the 'mpv2' java package (by K. Skretting)
% see page: http://www.ux.uis.no/~karlsk/dle/index.html
% They are all variants of Matching Pursuit
% 'javaMP', 'javaMatchingPursuit', 'javaBMP', 'javaBasicMatchingPursuit'
% 'javaOMP' or 'javaOrthogonalMatchingPursuit'
% 'javaORMP' or 'javaOrderRecursiveMatchingPursuit'
% 'javaPS' or 'javaPartialSearch'
% (4) Methods implemented as mex-files in SPAMS (by J. Mairal)
% These are very fast and can be used also for quite large problems.
% see page: http://spams-devel.gforge.inria.fr/
% 'mexLasso', 'LARS', or 'LASSO'
% 'mexOMP' NOTE: this is the same algorithm as ORMP and javaORMP!
%
% The Options may be:
% -------------------
% 'tnz' or 'targetNonZeros' with values as 1x1 or 1xL, gives the target
% number of non-zero coefficients for each column vector in X
% Default is ceil(N/2)
% 'tre' or 'targetRelativeError' with values as 1x1 or 1xL, gives the
% target relative error, i.e. iterations stops when ||r|| < tre*||x||.
% If both tnz and tre is given, the iterations stops when any criterium
% is met. Default is 1e-6.
% 'tae' or 'targetAbsoluteError' with values as 1x1 or 1xL, gives an
% alternative way to set 'tre' on: tre = tae ./ ||x||
% Iterations stops when ||r|| < tae. If used 'tae' overrides 'tre'.
% 'doOP' do Orthogonal Projection when thresholding, default is true
% 'nIt' or 'nofIt' or 'numberOfIterations' is used in FOCUSS, default 20.
% 'p' or 'pFOCUSS' the p-norm to use in FOCUSS. Default 0.5.
% 'l' or 'lambdaFOCUSS' is lambda in Regularized FOCUSS, default is 0
% 'nComb' number of combinations in 'javaPS'. Default is 20.
% 'globalReDist', 'tSSE' or 'targetSSE', 'tSNR' or 'targetSNR' are
% undocumented options which may be used with method 'javaORMP' only.
% 'GMPLoopSize', optional parameter used in GMP (default usually ok)
% 'paramSPAMS', optional parameter to use in mexLasso or mexOMP. If not given
% the following will be used for mexOMP and mexLasso respectively:
% paramSPAMS = struct( 'eps', tae.^2, 'L', int32(tnz) );
% paramSPAMS = struct( 'mode',1, 'lambda', tae.^2, 'L', int32(tnz), 'ols',true );
% 'v' or 'verbose' may be 0/false (default), 1/true or 2 (very verbose).
%
% Examples:
% ---------
% L = 400; N = 16; K = 32; D = randn(N,K); X = randn(N,L); % or
% % Wi = zeros(K,L); for i=1:L; Wi(:,i) = randperm(K); end;
% % Wt = zeros(K,L); Wt(Wi <= 5) = randn(5*L,1); X = D*Wt + 0.01*randn(N,L);
% op = struct('targetNonZeros',5, 'verbose', 2);
% [Wpinv, ra] = sparseapprox(X, D, 'pinv', op);
% [Wlp, rc] = sparseapprox(X, D, 'linprog', op);
% [Wf, rd] = sparseapprox(X, D, 'FOCUSS', op, 'p', 0.8, 'l', 0.4, 'nIt', 100);
% [Womp, rf] = sparseapprox(X, D, 'javaOMP', op);
% [Wormp, rg] = sparseapprox(X, D, 'javaORMP', op);
% [Wps100, rh] = sparseapprox(X, D, 'javaPS', 'nComb',100, op);
% [Wgmp, ri] = sparseapprox(X, D, 'GMP', 'tnz',5*L, 'v',2);
% [Wps, rih] = sparseapprox(X, D, 'javaPS', 'tnz',sum(Wgmp~=0), 'nComb',100);
% fs = ' %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f \n';
% fprintf('\n pinv linprog FOCUSS OMP ORMP PS GMP GMP+PS \n');
% fprintf(['SNR ',fs], ra.snr,rc.snr,rd.snr,rf.snr,rg.snr,rh.snr,ri.snr,rih.snr);
% fprintf(['time ',fs], ra.time,rc.time,rd.time,rf.time,rg.time,rh.time,ri.time,rih.time);
% [Wmomp, rm] = sparseapprox(X, D, 'mexOMP', op); % like ORMP
% [Wlasso, rl] = sparseapprox(X, D, 'mexLasso', op);
% % to show the coefficients for a selected vector
% % i = 11; [Wt(:,i), Womp(:,i), Wormp(:,i), Wmomp(:,i), Wlasso(:,i), Wps100(:,i)]
%----------------------------------------------------------------------
% Copyright (c) 2009. Karl Skretting. All rights reserved.
% University of Stavanger, Signal Processing Group
% Mail: karl.skretting@uis.no Homepage: http://www.ux.uis.no/~karlsk/
%
% HISTORY: dd.mm.yyyy
% Ver. 1.0 10.10.2009 Made function
% Ver. 1.1 04.01.2010 Add GMP method, and globalReDist in javaORMP
% Ver. 1.2 05.01.2010 added targetTotalSumSquaredError and targetSNR
% 06.01.2010 and now it works well.
% Ver. 1.3 26.03.2010 Add LARS and LASSO methods
% Ver. 1.4 07.04.2010 Added smoothed L_0 norm (SL0) method
% Ver. 2.0 06.04.2011 Removed some methods and tried to 'clean' this file
% SLO, LARS/LASSO matlab implementations, and
% globalReDist after javaORMP were removed.
% Ver. 2.1 31.05.2011 globalReDist included again
% Ver. 2.2 08.08.2012 Added affine ORMP method,
% Ver. 2.3 15.04.2013 Removed affine ORMP method (It was not as I hoped it would be)
% Ver. 2.4 09.10.2014 Minor changes to make mexOMP work (SPAMS ver. 2.5)
%----------------------------------------------------------------------
%
% additional documentation:
% - Dictionary Learning Tools: http://www.ux.uis.no/~karlsk/dle/index.html
% alternative functions:
% - GreedLab (sparsify), Thomas Blumensath et al. (Edinburgh)
% - SparseLab, David Donoho et al. (Stanford)
% - SPAMS (Mairal): http://spams-devel.gforge.inria.fr/
% - OMPBox (Ron Rubinstein): http://www.cs.technion.ac.il/~ronrubin/software.html
%%Check if input arguments are given
mfile = 'sparseapprox';
L = 400;
N = 16;
K = 32;
D = randn(N,K);
X = randn(N,L);
Wi = zeros(K,L);
for i=1:L; Wi(:,i) = randperm(K); end
Wt = zeros(K,L);
Wt(Wi <= 5) = randn(5*L,1);
%X = D*Wt + 0.01*randn(N,L);
op = struct('targetNonZeros',5, 'verbose', 2);
[Wpinv, ra] = sparse(X);
%[Wpinv, ra] = sparseapprox(X, D, 'pinv',op);
[Wlp, rc] = sparseapprox(X, D, 'linprog', op);
[Wf, rd] = sparseapprox(X, D, 'FOCUSS', op, 'p', 0.8, 'l', 0.4, 'nIt', 100);
[Womp, rf] = sparseapprox(X, D, 'javaOMP', op);
[Wormp, rg] = sparseapprox(X, D, 'javaORMP', op);
[Wps100, rh] = sparseapprox(X, D, 'javaPS', 'nComb',100, op);
[Wgmp, ri] = sparseapprox(X, D, 'GMP', 'tnz',5*L, 'v',2);
[Wps, rih] = sparseapprox(X, D, 'javaPS', 'tnz',sum(Wgmp~=0), 'nComb',100);
fs = ' %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f \n';
fprintf('\n pinv linprog FOCUSS OMP ORMP PS GMP GMP+PS \n');
fprintf(['SNR ',fs], ra.snr,rc.snr,rd.snr,rf.snr,rg.snr,rh.snr,ri.snr,rih.snr);
fprintf(['time ',fs], ra.time,rc.time,rd.time,rf.time,rg.time,rh.time,ri.time,rih.time);
[Wmomp, rm] = sparseapprox(X, D, 'mexOMP', op); % like ORMP
[Wlasso, rl] = sparseapprox(X, D, 'mexLasso', op);
% % to show the coefficients for a selected vector
i = 11; [Wt(:,i), Womp(:,i), Wormp(:,i), Wmomp(:,i), Wlasso(:,i), Wps100(:,i)]
if (nargin < 3) % just check number of input arguments
t = [mfile,': arguments must be given, see help.'];
disp(t);
if nargout >= 1
varargout{1} = -1;
end
if nargout >= 2
varargout{2} = struct('Error',t);
end
return
end
%%defaults, initial values
tstart = tic;
[N,L] = size(X);
K = size(D,2);
norm2X = sqrt(sum(X.*X)); % ||x(i)||_2 1xL
W = zeros(K,L); % the weights (coefficients)
tnz = ceil(N/2)*ones(1,L); % target number of non-zeros
thrActive = false; % is set to true if tnz, tre or tae is given
% and used for methods: pinv, backslash, linprog and
% FOCUSS
doOP = true; % do Orthogonal Projection when thresholding
relLim = 1e-6;
tre = relLim*ones(1,L); % target relative error: ||r|| <= tre*||x||
nComb = 20; % used only in javaPS
nIt = 20; % used only in FOCUSS
pFOCUSS = 0.5; % used only in FOCUSS
lambdaFOCUSS = 0; % used only in FOCUSS
deltaWlimit = 1e-8; % used only in FOCUSS
GMPLoopSize = 0; % used only in GMP
globalReDist = 0; % may be used with javaORMP
targetSSE = 0; % may be used with javaORMP
verbose = 0;
done = false;
javaClass = 'mpv2.MatchingPursuit'; % the important java class
spams_mex_file = 'mexLasso'; % one of the used SPAMS files
%%get the options
nofOptions = nargin-3;
optionNumber = 1;
fieldNumber = 1;
while (optionNumber <= nofOptions)
if isstruct(varargin{optionNumber})
sOptions = varargin{optionNumber};
sNames = fieldnames(sOptions);
opName = sNames{fieldNumber};
opVal = sOptions.(opName);
% next option is next field or next (pair of) arguments
fieldNumber = fieldNumber + 1; % next field
if (fieldNumber > numel(sNames))
fieldNumber = 1;
optionNumber = optionNumber + 1; % next pair of options
end
elseif iscell(varargin{optionNumber})
sOptions = varargin{optionNumber};
opName = sOptions{fieldNumber};
opVal = sOptions{fieldNumber+1};
% next option is next pair in cell or next (pair of) arguments
fieldNumber = fieldNumber + 2; % next pair in cell
if (fieldNumber > numel(sOptions))
fieldNumber = 1;
optionNumber = optionNumber + 1; % next pair of options
end
else
opName = varargin{optionNumber};
opVal = varargin{optionNumber+1};
optionNumber = optionNumber + 2; % next pair of options
end
% interpret opName and opVal
if strcmpi(opName,'targetNonZeros') || strcmpi(opName,'tnz')
if strcmpi(met,'GMP')
tnz = opVal; % GMP will distribute the non-zeros
else
if numel(opVal)==1
tnz = opVal*ones(1,L);
elseif numel(opVal)==L
tnz = reshape(opVal,1,L);
else
error([mfile,': illegal size of value for option ',opName]);
end
end
thrActive = true;
end
if strcmpi(opName,'targetRelativeError') || strcmpi(opName,'tre')
if numel(opVal)==1
tre = opVal*ones(1,L);
elseif numel(opVal)==L
tre = reshape(opVal,1,L);
else
error([mfile,': illegal size of value for option ',opName]);
end
thrActive = true;
end
if strcmpi(opName,'targetAbsoluteError') || strcmpi(opName,'tae')
if numel(opVal)==1
tae = opVal*ones(1,L);
elseif numel(opVal)==L
tae = reshape(opVal,1,L);
else
error([mfile,': illegal size of value for option ',opName]);
end
thrActive = true;
end
if ( strcmpi(opName,'nIt') || strcmpi(opName,'nofIt') || ...
strcmpi(opName,'numberOfIterations') )
if (isnumeric(opVal) && numel(opVal)==1)
nIt = max(floor(opVal), 1);
else
error([mfile,': illegal size of value for option ',opName]);
end
end
if strcmpi(opName,'p') || strcmpi(opName,'pFOCUSS')
if (isnumeric(opVal) && numel(opVal)==1)
pFOCUSS = min(opVal, 1);
else
error([mfile,': illegal size of value for option ',opName]);
end
end
if strcmpi(opName,'l') || strcmpi(opName,'lambdaFOCUSS')
if (isnumeric(opVal) && numel(opVal)==1)
lambdaFOCUSS = abs(opVal);
else
error([mfile,': illegal size of value for option ',opName]);
end
end
if strcmpi(opName,'nComb')
if (isnumeric(opVal) && numel(opVal)==1)
nComb = max(floor(opVal), 2);
else
error([mfile,': illegal size of value for option ',opName]);
end
end
if strcmpi(opName,'paramSPAMS')
if (isstruct(opVal))
paramSPAMS = opVal;
else
error([mfile,': option paramSPAMS is not a struct as it should be, see SPAMS help.']);
end
end
if strcmpi(opName,'globalReDist')
if (isnumeric(opVal) && numel(opVal)==1)
globalReDist = min(max(floor(opVal), 0), 2); % 0, 1 or 2
else
error([mfile,': illegal size of value for option ',opName]);
end
end
if strcmpi(opName,'doOP')
if (islogical(opVal)); doOP = opVal; end;
if isnumeric(opVal); doOP = (opVal ~= 0); end;
end
if strcmpi(opName,'GMPLoopSize')
if (isnumeric(opVal) && numel(opVal)==1)
GMPLoopSize = max(floor(opVal), 2);
else
error([mfile,': illegal size of value for option ',opName]);
end
end
if strcmpi(opName,'tSSE') || strcmpi(opName,'targetSSE')
if (isnumeric(opVal) && numel(opVal)==1)
targetSSE = min(max(opVal, 0), sum(sum(X.*X)));
else
error([mfile,': illegal size of value for option ',opName]);
end
end
if strcmpi(opName,'tSNR') || strcmpi(opName,'targetSNR')
if (isnumeric(opVal) && numel(opVal)==1)
targetSSE = 10^(-abs(opVal)/10) * sum(sum(X.*X));
else
error([mfile,': illegal size of value for option ',opName]);
end
end
if strcmpi(opName,'verbose') || strcmpi(opName,'v')
if (islogical(opVal) && opVal); verbose = 1; end;
if isnumeric(opVal); verbose = opVal(1); end;
end
end
if exist('tae','var') % if both exist 'tae' overrules 'tre'
tre = tae./norm2X;
elseif exist('tre','var') % 'tre' was given a default value
tae = tre.*norm2X;
else % so this case is redundant
disp(' ??? This is never printed.');
end
%%Display info
if (verbose > 1) % very verbose
disp(' ');
disp([mfile,' with method ',met,' started ',datestr(now)]);
disp(['Size of X is ',int2str(size(X,1)),'x',int2str(size(X,2)),...
', D is ',int2str(size(D,1)),'x',int2str(size(D,2)),...
', and W is ',int2str(size(W,1)),'x',int2str(size(W,2))]);
end
Method of Frames
if strcmpi(met,'MOF') || strcmpi(met,'MethodOfFrames') || strcmpi(met,'pinv')
textMethod = 'Method of Frames with pseudoinverse (pinv).';
if (verbose >= 1); disp([mfile,': ',textMethod]); end;
W = pinv(D)*X;
if thrActive % then adjust w by setting more to zero
W = setSmallWeightsToZero(X,D,W,tnz,tae,doOP);
end
done = true;
end
Backslash method just find an exact solution with N non-zeros
if strcmpi(met,'BackSlash') || strcmpi(met,'\')
textMethod = 'Matlab backslash operator.';
if (verbose >= 1); disp([mfile,': ',textMethod]); end;
W = D\X;
if thrActive % then adjust w by setting more to zero
W = setSmallWeightsToZero(X,D,W,tnz,tae,doOP);
end
done = true;
end
Basis Pursuit
if strcmpi(met,'BP') || strcmpi(met,'BasisPursuit') || strcmpi(met,'linprog')
textMethod = 'Basis Pursuit with Matlab function linprog.';
if (verbose >= 1); disp([mfile,': ',textMethod]); end;
f = ones(2*K,1);
A = [D,-D];
LB = zeros(2*K,1); % lower bound for w
UB = ones(2*K,1)*inf; % upper bound for w
for columnNumber = 1:L
x = X(:,columnNumber);
w2 = linprog(f,A,x,A,x,LB,UB); % minimize 1-norm of w
W(:,columnNumber) = w2(1:K)-w2((1+K):(2*K));
end
if thrActive % then adjust w by setting more to zero
W = setSmallWeightsToZero(X,D,W,tnz,tae,doOP);
end
done = true;
end
the original (and possible regularized) FOCUSS algorithm
if strcmpi(met,'FOCUSS')
W = pinv(D)*X; % initial values
if (lambdaFOCUSS > 0) % Regularized FOCUSS
textMethod = ['FOCUSS with p=',num2str(pFOCUSS),...
', and regularization (lambda = ',....
num2str(lambdaFOCUSS),')',...
' and ',int2str(nIt),' iterations.'];
else
textMethod = ['FOCUSS with p=',num2str(pFOCUSS),...
' and ',int2str(nIt),' iterations.'];
end
if thrActive
textMethod = char(textMethod, ...
' Thresholding is done in the end.');
end
if (nargout >= 2) % keep track of sparseness
sparseInW = zeros(5,nIt);
edges = [0, 0.0001, 0.001, 0.01, 0.1, inf];
changeInW = zeros(nIt,L);
numberOfIterations = nIt*ones(1,L);
end
if (verbose >= 1); disp(char([mfile,': '],textMethod)); end;
for columnNumber = 1:L
w = W(:,columnNumber);
w0 = w;
x = X(:,columnNumber);
for i=1:nIt
Qdiagonal = abs(w).^(1-pFOCUSS/2);
F = D.*(ones(N,1)*(Qdiagonal')); % F^{k+1} in Engan's PhD
if (lambdaFOCUSS > 0) % Regularized FOCUSS
q = F'*( ((F*F'+lambdaFOCUSS*eye(N)) \ x) );
else % original FOCUSS
q = pinv(F) * x;
end
w = Qdiagonal.*q;
change = norm(w-w0);
if (nargout >= 2) % keep track of sparseness and more
m = max(abs(w));
I = histc(abs(w)/m, edges);
sparseInW(:,i) = sparseInW(:,i)+I(1:5);
changeInW(i,columnNumber) = change;
if (change < deltaWlimit)
sparseInW(:,(i+1):nIt) = sparseInW(:,(i+1):nIt) + I(1:5)*ones(1,nIt-i);
changeInW((i+1):nIt,columnNumber) = change;
numberOfIterations(columnNumber) = i;
end
end
if (change < deltaWlimit); break; end;
w0 = w;
end
W(:,columnNumber) = w;
end
if thrActive % then adjust w by setting more to zero
W = setSmallWeightsToZero(X,D,W,tnz,tae,doOP);
end
done = true;
end
mexOMP or mexLasso algorithm from SPAMS
if (strcmpi(met,'mexOMP') || strcmpi(met,'mexORMP') || ...
strcmpi(met,'mexLasso') || strcmpi(met,'LARS') || strcmpi(met,'LASSO'))
% check matlab version
t = version('-release');
if (eval(t(1:4)) < 2009) || strcmpi(t,'2009a')
t = [mfile,': mexLasso and mexOMP need Matlab version >= 2009b. (?)'];
disp(t);
if nargout >= 1
varargout{1} = -1;
end
if nargout >= 2
varargout{2} = struct('Error',t);
end
return
end
%
% The way access to SPAMS is checked on the computeres I may use
if ~(exist(spams_mex_file,'file') == 3) % mex-file not available (yet)
start_spams; % a m-file that comes with SPAMS and is located on a common path
end
if ~(exist(spams_mex_file,'file') == 3) % mex-file still not available
t = [mfile,': can not access mexLasso and mexOMP on this computer.'];
disp(t);
if nargout >= 1
varargout{1} = -1;
end
if nargout >= 2
varargout{2} = struct('Error',t);
end
return
end
%
if (strcmpi(met,'mexOMP') || strcmpi(met,'mexORMP'))
textMethod = 'mexOMP in SPAMS package (by J. Mairal).';
else
textMethod = 'mexLasso (mode=1) in SPAMS package (by J. Mairal).';
end
if (verbose >= 1); disp([mfile,': ',textMethod]); end;
if ~exist('paramSPAMS','var')
if (strcmpi(met,'mexOMP') || strcmpi(met,'mexORMP'))
% Example from SPAMS documentation, mexOMP
% parameter of the optimization procedure are chosen
% param.L=10; % not more than 10 non-zeros coefficients
% param.eps=0.1; % squared norm of the residual should be less than 0.1
% param.numThreads=-1; % number of processors/cores to use; the default choice is -1
% and uses all the cores of the machine
paramSPAMS = struct(...
'eps', tae.^2, ...
'L', int32(tnz), ...
'numThreads', -1 );
else
paramSPAMS = struct(...
'mode', 1, ...
'lambda', tae.^2, ...
'L', int32(tnz), ...
'ols', true );
end
if (verbose >= 1); disp(' and uses ''default'' param when calling SPAMS function.'); end;
else
if (verbose >= 1); disp(' and uses user supplied param when calling SPAMS function.'); end;
end
%
d_norm = sqrt(sum(D.^2));
if sum(abs(d_norm - ones(1,K))) > 1e-6
D = D./repmat(d_norm,[N 1]);
end
%
if (strcmpi(met,'mexOMP') || strcmpi(met,'mexORMP'))
W = mexOMP(X, D, paramSPAMS);
else
% x_norm = sqrt(sum(X.^2));
% X = X./repmat(x_norm,[N 1]);
W = mexLasso(X, D, paramSPAMS);
% W = W.*repmat(x_norm,[K 1]);
% X = X.*repmat(x_norm,[N 1]);
end
W = full(W);
if sum(abs(d_norm - ones(1,K))) > 1e-6
W = W./repmat(d_norm(:),[1 L]);
D = D.*repmat(d_norm,[N 1]);
end
done = true;
end
%%Algorithms implemented in Java
if strcmpi(met(1:min(numel(met),3)),'aff')
disp('Affine approximation is essentiallaly that sum of coefficients should be 1.');
disp('This is achieved by adding a row of constants to both D and X and ');
disp('then doing ordinary aparse approximation. Do this outside this function ');
disp('to have better control. Set a to an approriate value, 5 perhaps.');
disp('Ex: D = [a*ones(1,K); D]; X = [a*ones(1,L); X];');
% met = ['java',met(4:end)]; % use java method
% D = [5*ones(1,K); D];
% X = [5*ones(1,L); X];
done = false;
end
if strcmpi(met(1:min(numel(met),4)),'java')
if (not(exist(javaClass,'class')) && exist('java_access.m','file'))
java_access;
end
if (not(exist(javaClass,'class')) && exist('javaAccess.m','file'))
javaAccess; % an older version of java_access (may work if it exist)
end
if not(exist(javaClass,'class'))
javaErrorMessage = ['No access to ',javaClass,' in static or dynamic Java path.'];
disp(javaErrorMessage);
met = met(5:end);
disp(['Use method ',met,' instead.']);
else
jD = mpv2.SimpleMatrix( D );
if (L == 1)
jMP = mpv2.MatchingPursuit(jD);
else
jDD = mpv2.SymmetricMatrix(K,K);
jDD.eqInnerProductMatrix(jD);
jMP = mpv2.MatchingPursuit(jD,jDD);
end
end
end
if ( strcmpi(met,'javaMP') || strcmpi(met,'javaMatchingPursuit') || ...
strcmpi(met,'javaBMP') || strcmpi(met,'javaBasicMatchingPursuit') )
textMethod = 'Basic Matching Pursuit, Java implementation.';
if (verbose >= 1); disp([mfile,': ',textMethod]); end;
% note the 'tre' is not used for BMP
for j=1:L
if (tnz(j) > 0)
W(:,j) = jMP.vsBMP(X(:,j), int32(tnz(j)));
end
end
done = true;
end
if strcmpi(met,'javaOMP') || strcmpi(met,'javaOrthogonalMatchingPursuit')
textMethod = 'Orthogonal Matching Pursuit, Java implementation.';
if (verbose >= 1); disp([mfile,': ',textMethod]); end;
for j=1:L
if (tnz(j) > 0) && (tre(j) < 1)
W(:,j) = jMP.vsOMP(X(:,j), int32(tnz(j)), tre(j));
end
end
done = true;
end
if strcmpi(met,'javaORMP') || strcmpi(met,'javaOrderRecursiveMatchingPursuit')
%
% This could be as simple as javaOMP, but since globalReDist was
% reintroduced it is now quite complicated here.
if (targetSSE > 0)
% This is initialization of tre (and tnz ?) for the special case of
% global distribution of non-zeros where a target sum og squared
% errors is given as an input argument.
% Perhaps tnz also should be set to an appropriate value
% tnz = 2*ones(1,L);
tre = sqrt(targetSSE/L)./norm2X;
globalReDist = 2;
textMethod = ['javaORMP with global distribution of non-zeros ',...
'given target SSE (or SNR).'];
elseif (globalReDist == 1)
textMethod = ['javaORMP with global distribution of non-zeros ',...
'keeping the total number of non-zeros fixed.'];
elseif (globalReDist == 2)
textMethod = ['javaORMP with global distribution of non-zeros ',...
'keeping the total SSE fixed.'];
else
textMethod = 'Order Recursive Matching Pursuit, Java implementation.';
end
%
% below is the javaORMP lines
if (verbose >= 1); disp([mfile,': ',textMethod]); end;
for j=1:L
if (tnz(j) > 0) && (tre(j) < 1)
W(:,j) = jMP.vsORMP(X(:,j), int32(tnz(j)), tre(j));
end
end
%
% below is the globalReDist lines
% ******* START Global distribution of non-zeros.*****
% The structure is:
% 1. Initializing: Sm1 <= S <= Sp1 and SEm1 >= SE >= SEp1
% 2. Add atoms until SSE is small enough
% 3. or remove atoms until SSE is large enough
% 4. Add one atom as long as one (or more) may be removed and the
% SSE is reduced
if (globalReDist > 0)
% part 1
R = X - D * W; % representation error
S = sum(W ~= 0); % selected number of non-zeros for each column
SE = sum(R.*R); % squared error for each (when S is selected)
sumSinit = sum(S);
SSE = sum(SE);
SSEinit = SSE; % store initial SSE
Sp1 = S + 1; % selected number of non-zeros plus one
Sp1(Sp1 > N) = N;
Sm1 = S - 1; % selected number of non-zeros minus one
Sm1(Sm1 < 0) = 0;
SEp1 = zeros(1,L); % initializing corresponding squared error
SEm1 = zeros(1,L);
for j=1:L
x = X(:,j);
if Sp1(j) == S(j) % == N
w = W(:,j);
else
w = jMP.vsORMP(x, Sp1(j), relLim);
end
r = x-D*w;
SEp1(j) = r'*r;
if Sm1(j) == 0
w = zeros(K,1);
else
w = jMP.vsORMP(x, Sm1(j), relLim);
end
r = x-D*w;
SEm1(j) = r'*r;
end
SEdec = SE-SEp1; % the decrease in error by selectiong one more
SEinc = SEm1-SE; % the increase in error by selectiong one less
SEinc(S == 0) = inf; % not possible to select fewer than zero
addedS = 0;
removedS = 0;
addedSE = 0;
removedSE = 0;
[valinc, jinc] = min(SEinc); % min increase in SE by removing one atom
[valdec, jdec] = max(SEdec); % max reduction in SE by adding one atom
if (targetSSE > 0)
if (SSEinit > targetSSE) % part 2
if (verbose > 2)
disp(['(part 2 add atoms, target SSE = ',num2str(targetSSE),...
' and initial SSE = ',num2str(SSEinit),')']);
end
while (SSE > targetSSE)
j = jdec; % an atom is added to vector j
addedS = addedS+1;
removedSE = removedSE + valdec;
SSE = SSE - valdec;
% shift in Sm1,S,Sp1 and SEm1,SE,SEp1
[Sm1(j), S(j), Sp1(j)] = assign3(S(j), Sp1(j), min(Sp1(j)+1, N));
[SEm1(j), SE(j)] = assign2(SE(j), SEp1(j)); % and SEp1(j)=SEp1(j)
if (Sp1(j) > S(j)) % the normal case, find new SEp1(j)
w = jMP.vsORMP(X(:,j), Sp1(j), relLim);
r = X(:,j) - D*w;
SEp1(j) = r'*r;
end
SEinc(j) = SEdec(j); % SE cost by removing this again
SEdec(j) = SE(j)-SEp1(j); % SE gain by adding one more atom
%
W(:,j) = jMP.vsORMP(X(:,j), S(j), relLim);
[valdec, jdec] = max(SEdec);
end
[valinc, jinc] = min(SEinc);
elseif ((SSEinit+valinc) < targetSSE) % part 3
if (verbose > 2)
disp(['(part 3 remove atoms, target SSE = ',num2str(targetSSE),...
' and initial SSE = ',num2str(SSEinit),')']);
end
while ((SSE+valinc) < targetSSE)
j = jinc; % an atom is removed from vector j
removedS = removedS+1;
addedSE = addedSE + valinc;
SSE = SSE + valinc;
% shift in Sm1,S,Sp1 and SEm1,SE,SEp1
[Sm1(j), S(j), Sp1(j)] = assign3(max(Sm1(j)-1, 0), Sm1(j), S(j));
[SE(j), SEp1(j)] = assign2(SEm1(j), SE(j)); % and SEm1(j)=SEm1(j)
if (Sm1(j) > 0)
w = jMP.vsORMP(X(:,j), Sm1(j), relLim);
r = X(:,j) - D*w;
else
r = X(:,j);
end
SEm1(j) = r'*r;
%
SEdec(j) = SEinc(j); % SE gain by adding this atom again
if (S(j) > 0) % SE cost by removing another atom
W(:,j) = jMP.vsORMP(X(:,j), S(j), relLim);
SEinc(j) = SEm1(j)-SE(j);
else
W(:,j) = zeros(K,1);
SEinc(j) = inf; % can not select fewer and increase error
end
[valinc, jinc] = min(SEinc);
end
[valdec, jdec] = max(SEdec);
else %
if (verbose > 2)
disp(['(target SSE = ',num2str(targetSSE),...
' is close to initial SSE = ',num2str(SSEinit),')']);
end
end
else
targetSSE = SSEinit;
end
%
% part 4
while ((valinc < valdec) && (jinc ~= jdec))
j = jdec;
addedS = addedS+1;
removedSE = removedSE + valdec;
SSE = SSE - valdec;
% shift in Sm1,S,Sp1 and SEm1,SE,SEp1
[Sm1(j), S(j), Sp1(j)] = assign3(S(j), Sp1(j), min(Sp1(j)+1, N));
[SEm1(j), SE(j)] = assign2(SE(j), SEp1(j)); % and SEp1(j)=SEp1(j)
if (Sp1(j) > S(j)) % the normal case, find new SEp1(j)
w = jMP.vsORMP(X(:,j), Sp1(j), relLim);
r = X(:,j) - D*w;
SEp1(j) = r'*r;
end
SEinc(j) = SEdec(j); % SE cost by removing this again
SEdec(j) = SE(j)-SEp1(j); % SE gain by adding one more atom
W(:,j) = jMP.vsORMP(X(:,j), S(j), relLim);
[valinc, jinc] = min(SEinc);
%
while ((SSE+valinc) < targetSSE)
j = jinc;
removedS = removedS+1;
addedSE = addedSE + valinc;
SSE = SSE + valinc;
% shift in Sm1,S,Sp1 and SEm1,SE,SEp1
[Sm1(j), S(j), Sp1(j)] = assign3(max(Sm1(j)-1, 0), Sm1(j), S(j));
[SE(j), SEp1(j)] = assign2(SEm1(j), SE(j)); % and SEm1(j)=SEm1(j)
if (Sm1(j) > 0)
w = jMP.vsORMP(X(:,j), Sm1(j), relLim);
r = X(:,j) - D*w;
else
r = X(:,j);
end
SEm1(j) = r'*r;
%
SEdec(j) = SEinc(j); % SE gain by adding this atom again
if (S(j) > 0) % SE cost by removing another atom
W(:,j) = jMP.vsORMP(X(:,j), S(j), relLim);
SEinc(j) = SEm1(j)-SE(j);
else
W(:,j) = zeros(K,1);
SEinc(j) = inf; % can not select fewer and increase error
end
[valinc, jinc] = min(SEinc);
if (globalReDist == 1); break; end;
end
[valdec, jdec] = max(SEdec); % next now
end
%
if (verbose > 2)
disp(['Using globalReDist=',int2str(globalReDist),': ',...
'non-zeros in W changed as ',int2str(sumSinit),...
' + ',int2str(addedS),' - ',int2str(removedS),...
' = ',int2str(sum(sum(W ~= 0)))]);
disp(['SSE changed as ',num2str(SSEinit),...
' + ',num2str(addedSE),' - ',num2str(removedSE),...
' = ',num2str(SSE),' = ',num2str(SSEinit+addedSE-removedSE)]);
disp(['(target SSE = ',num2str(targetSSE),...
' and actual SSE = ',num2str(sum(sum((X - D * W).^2))),')']);
end
end
%
% ******* END Global distribution of non-zeros.*****
%
%
done = true;
end
if strcmpi(met,'javaPS') || strcmpi(met,'javaPartialSearch')
textMethod = ['Partial Search with ',...
int2str(nComb),' number of combinations.'];
if (verbose >= 1); disp([mfile,': ',textMethod]); end;
for j=1:L
if (tnz(j) > 0) && (tre(j) < 1)
if (tnz(j) < 2)
W(:,j) = jMP.vsORMP(X(:,j), int32(tnz(j)), tre(j));
else
W(:,j) = jMP.vsPS(X(:,j), int32(tnz(j)), tre(j), int32(nComb));
end
end
end
done = true;
end
%%GMP is Global Matching Pursuit
if strcmpi(met,'GMP')
if (GMPLoopSize <= 0)
GMPLoopSize = floor(tnz/N);
end
if (GMPLoopSize > 0.9*L)
GMPLoopSize = floor(0.9*L);
end
textMethod = ['Global Matching Pursuit with ',...
int2str(tnz),' non-zeros, (N=',int2str(N),...
', L=',int2str(L),' GMPLoopSize=',int2str(GMPLoopSize),').'];
if verbose
disp(textMethod);
end
Gd = (1./sqrt(sum(D.*D))); % G = diag(Gd)
F = D.*(ones(size(D,1),1)*Gd); % normalize D, F = D*G
%
nzW = 0;
R = X;
while (nzW < tnz)
if (verbose > 2)
disp(['GMP: nzW = ',int2str(nzW),' and ||R||^2 = ',...
num2str(sum(sum((R.*R))))]);
end
U = F'*R; % inner products = G*D'*R
[um,iK] = max(abs(U));
[temp,iL] = sort(um,2,'descend'); %#ok<ASGLU>
for i=1:min(tnz-nzW, GMPLoopSize)
il = iL(i);
ik = iK(il);
W(ik,il) = W(ik,il) + U(ik,il)*Gd(ik);
end
nzW = sum(W(:) ~= 0);
R = X - D*W;
end
if (verbose > 1)
disp(['GMP: nzW = ',int2str(nzW),' and ||R||^2 = ',...
num2str(sum(sum((R.*R))))]);
end
done = true;
end
%%OMP and ORMP are almost equal.
% function is from FrameTools\VSormp.m (ver 12.06.2003 with some few
% modifications). Variables have short names so there is a risk for
% variable mixup when using the algoritm within a larger context.
if strcmpi(met,'OMP') || strcmpi(met,'ORMP')
if strcmpi(met,'OMP')
textMethod = 'Orthogonal Matching Pursuit.';
end
if strcmpi(met,'ORMP')
textMethod = 'Order Recursive Matching Pursuit.';
end
if (verbose >= 1); disp([mfile,': ',textMethod]); end;
F = D.*(ones(size(D,1),1)*(1./sqrt(sum(D.*D)))); % normalize D
FF = F'*F;
for columnNumber = 1:L % i.e. for each data vector
% ********************** INITIALIZE **********************
x = X(:,columnNumber);
S = tnz(columnNumber);
r = zeros(S,K);
w = zeros(K,1);
T = 1:K;
e = ones(K,1);
u = ones(K,1);
c = x'*F;
n2x = x'*x;
n2xLim = n2x*tre(columnNumber);
% select the first frame vector
[cm,km] = max(abs(c));
s = 1;
J = km;
T(km) = -1;
r(1,km) = u(km);
n2x = n2x-cm*cm;
% ********************** THE MAIN LOOP **********************
while ((s<S) && (n2x>n2xLim))
for k=1:K
if (T(k)>=0)
r(s,k) = FF(km,k);
for n=1:(s-1)
r(s,k) = r(s,k)-r(n,km)*r(n,k);
end
if (u(km)~=0); r(s,k) = r(s,k)/u(km); end;
c(k) = c(k)*u(k)-c(km)*r(s,k);
if strcmpi(met,'OMP') % use next line for OMP
w(k) = abs(c(k)); % use w here (instead of a new variable d)
end
e(k) = e(k)-r(s,k)*r(s,k);
u(k) = sqrt(abs(e(k))); % abs kun i matlab!
if (u(k)~=0); c(k) = c(k)/u(k); end;
if strcmpi(met,'ORMP') % use next line for ORMP
w(k) = abs(c(k)); % use w here (instead of a new variable d)
end
end
end
w(km) = 0; % w(J) = 0;
% select the next frame vector
[temp,km] = max(abs(w)); %#ok<ASGLU>
s = s+1;
J(s) = km;
T(km) = -1;
r(s,km) = u(km);
cm = c(km);
n2x = n2x-cm*cm;
end % ******** END OF MAIN LOOP **********************************
% ************ BACK-SUBSTITUTION *************
w = zeros(K,1);
for k=s:(-1):1
Jk=J(k);
for n=s:(-1):(k+1)
c(Jk) = c(Jk)-c(J(n))*r(k,J(n));
end
if (r(k,Jk)~=0); c(Jk) = c(Jk)/r(k,Jk); end;
w(Jk) = c(Jk);
end
%
W(:,columnNumber) = w;
end
W = W .* ((1./sqrt(sum(D.*D)))'*ones(1,L)); % rescale W
done = true;
end
%%Now we are finished, 'done' should be true
% but test this before finding (and/or) displaying results
W = full(W);
% W = sparse(W); % is a good alternative
%%may display info before returning
if done && ((verbose > 1) || (nargout >= 2)) % need some results
R = X - D*W;
varX = var(X(:));
varR = var(R(:));
if (varR > 0)
snr = 10*log10(varX/varR);
else
snr = inf;
end
norm0X = sum(X ~= 0);
norm1X = sum(abs(X));
normiX = max(abs(X));
norm0R = sum(R ~= 0);
norm1R = sum(abs(R));
norm2R = sqrt(sum(R.*R));
normiR = max(abs(R));
norm0W = sum(W ~= 0);
norm1W = sum(abs(W));
norm2W = sqrt(sum(W.*W));
normiW = max(abs(W));
end
if done && (verbose >= 2) % very verbose
if (snr < 100)
disp([mfile,': SNR for the reconstruction is ',...
num2str(snr,'%7.4f')]);
elseif (snr < 500)
disp([mfile,': Almost perfect reconstruction, SNR > 100.']);
else
disp([mfile,': Perfect reconstruction, X = D*W.']);
end
disp(['Number of non-zeros in W is ',int2str(sum(norm0W)),...
' (sparseness factor is ',num2str(sum(norm0W)/(N*L)),')']);
if exist('numberOfIterations','var');
disp(['Average number of iterations for each column : ',...
num2str(mean(numberOfIterations),'%5.1f')]);
end
%
disp(['X: ',num2str(min(norm0X)),' <= ||x||_0 <= ',...
num2str(max(norm0X)),' and mean is ',num2str(mean(norm0X))]);
disp([' ',num2str(min(norm1X)),' <= ||x||_1 <= ',...
num2str(max(norm1X)),' and mean is ',num2str(mean(norm1X))]);
disp([' ',num2str(min(norm2X)),' <= ||x||_2 <= ',...
num2str(max(norm2X)),' and mean is ',num2str(mean(norm2X))]);
disp([' ',num2str(min(normiX)),' <= ||x||_inf <= ',...
num2str(max(normiX)),' and mean is ',num2str(mean(normiX))]);
disp(['R: ',num2str(min(norm0R)),' <= ||r||_0 <= ',...
num2str(max(norm0R)),' and mean is ',num2str(mean(norm0R))]);
disp([' ',num2str(min(norm1R)),' <= ||r||_1 <= ',...
num2str(max(norm1R)),' and mean is ',num2str(mean(norm1R))]);
disp([' ',num2str(min(norm2R)),' <= ||r||_2 <= ',...
num2str(max(norm2R)),' and mean is ',num2str(mean(norm2R))]);
disp([' ',num2str(min(normiR)),' <= ||r||_inf <= ',...
num2str(max(normiR)),' and mean is ',num2str(mean(normiR))]);
disp(['W: ',num2str(min(norm0W)),' <= ||w||_0 <= ',...
num2str(max(norm0W)),' and mean is ',num2str(mean(norm0W))]);
disp([' ',num2str(min(norm1W)),' <= ||w||_1 <= ',...
num2str(max(norm1W)),' and mean is ',num2str(mean(norm1W))]);
disp([' ',num2str(min(norm2W)),' <= ||w||_2 <= ',...
num2str(max(norm2W)),' and mean is ',num2str(mean(norm2W))]);
disp([' ',num2str(min(normiW)),' <= ||w||_inf <= ',...
num2str(max(normiW)),' and mean is ',num2str(mean(normiW))]);
disp(' ');
disp([mfile,' with ',met,' done. Used time is ',num2str(toc(tstart))]);
disp(' ');
end
%%Return Outputs
if done
if (nargout >= 1)
varargout{1} = W;
end
if (nargout >= 2)
varargout{2} = struct( 'time', toc(tstart), 'snr', snr, ...
'textMethod', textMethod, ...
'norm0X', norm0X, 'norm1X', norm1X, ...
'norm2X', norm2X, 'normiX', normiX, ...
'norm0R', norm0R, 'norm1R', norm1R, ...
'norm2R', norm2R, 'normiR', normiR, ...
'norm0W', norm0W, 'norm1W', norm1W, ...
'norm2W', norm2W, 'normiW', normiW );
% extra output arguments for standard and regularized FOUCUSS
if exist('sparseInW','var');
varargout{2}.sparseInW = sparseInW;
varargout{2}.edges = edges;
end
if exist('changeInW','var');
varargout{2}.changeInW = changeInW;
end
if exist('numberOfIterations','var');
varargout{2}.numberOfIterations = numberOfIterations;
end
if exist('sol','var');
varargout{2}.sol = sol;
end
if exist('paramSPAMS','var');
varargout{2}.paramSPAMS = paramSPAMS;
end
end
else % ~done did not find W
textMethod = ['Unknown method ',met,'. Coefficients W not found.'];
disp([mfile,': ',textMethod]);
if (nargout >= 1)
varargout{1} = [];
end
if (nargout >= 2)
varargout{2} = struct( 'time', toc(tstart), ...
'textMethod', textMethod );
end
end
if exist('javaErrorMessage','var');
varargout{2}.Error = javaErrorMessage;
end
return
function W = setSmallWeightsToZero(X,D,W,tnz,tae,doOP)
[K,L] = size(W); % K=size(D,2), L=size(X,2)
for i = 1:L
[absw,I] = sort(abs(W(:,i)),'descend');
w = zeros(K,1);
for s = 1:tnz(i)
if (absw(s) < 1e-10); break; end;
if doOP % do orthogonal projection onto the selected columns of D
Di = D(:,I(1:s));
w(I(1:s)) = (Di'*Di)\(Di'*X(:,i));
else % simple thresholding keeps the largest coefficients unaltered
w(I(s)) = W(I(s),i);
end
r = X(:,i) - D*w;
if (sqrt(r'*r) < tae(i)); break; end;
end
W(:,i) = w;
end
return
function [a, b] = assign2(c, d)
a = c;
b = d;
return
function [a, b, c] = assign3(d, e, f)
a = d;
b = e;
c = f;
return

Sign in to answer this question.

Answers (1)

OCDER
OCDER on 22 Aug 2018
Edited: OCDER on 22 Aug 2018
The code is very weird. It immediately summons itself with no switch to break the recursion it seems. Also, inputs D and X are just overridden immediately with randn(N, k) and randn(N, L).
Essentially, it does this as shown below. Not sure how to fix it to give the "correct" results.
function varargout = sparseapprox(X, D, met, varargin)
L = 400;
N = 16;
K = 32;
D = randn(N,K);
X = randn(N,L);
Wi = zeros(K,L);
for i=1:L; Wi(:,i) = randperm(K); end
op = struct('targetNonZeros',5, 'verbose', 2);
[Wlp, rc] = sparseapprox(X, D, 'linprog', op);
if (nargout >= 1)
varargout{1} = Wlb;
end
if (nargout >= 2)
varargout{2} = struct('error', 'failed');
end
Try contacting the author of this work for help, as this is too complex to fix "correctly", as in we could stop the infinite recursion, but don't know when it was suppose to do stop the recursion.
Copyright (c) 2009. Karl Skretting. All rights reserved.
University of Stavanger, Signal Processing Group
Mail: karl.skretting@uis.no Homepage: http://www.ux.uis.no/~karlsk/
-----------------------------------------------------
It is unlikely that we'll read every code you've uploaded in your post. The best answer is probably just the link below:
https://www.mathworks.com/help/matlab/matlab_prog/resolving-out-of-memory-errors.html
See this link to begin your debugging journey.
https://www.mathworks.com/help/matlab/matlab_prog/debugging-process-and-features.html
If this is not the answer you seek, read this before editing your question:
https://www.mathworks.com/matlabcentral/answers/6200-tutorial-how-to-ask-a-question-on-answers-and-get-a-fast-answer

Categories

Find more on Units of Measurement in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!