How to solve problem of out memory

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:  [e-mail]   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

Answer 1

[OCDER]

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:  [e-mail]   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