Sparsity

gurls/gurls_defopt.m

@@ -76,4 +76,10 @@
     opt.newprop( 'jobid', 1);
     opt.newprop( 'seq', {});
     opt.newprop( 'process', {});
+    %% ISTA options
+    opt.newprop( 'ISTAalpha',1);
+    opt.newprop( 'ISTAniter',inf);
+    opt.newprop( 'ISTArelthre',1e-4);
+    opt.newprop( 'ISTASparsitylvl',inf);
+    opt.newprop( 'ISTASAccuReq',0);
 end

gurls/optimizers/rls_fista.m

@@ -0,0 +1,80 @@
+function [cfr] = rls_fista (X, y, opt)
+% rls_fista(X,y,opt)
+% computes a classifier for elastic nerwork using FISTA.
+% The regularization parameter is set to the one found in opt.paramsel.
+% In case of multiclass problems, the regularizers need to be combined with the opt.singlelambda function.
+%
+% INPUTS:
+% -OPT: struct of options with the following fields:
+%   fields that need to be set through previous gurls tasks:
+%		- paramsel.lambdas (set by the paramsel_* routines)
+%   fields with default values set through the defopt function:
+%		- singlelambda
+%   fields with default values set through the defopt function:
+%       - ISTAalpha (paramters for balancing l1-norm and l-2 norm. 1 for LASSO and 0 for ridge)
+%       - niter (maximun number for iteration. Set to either negative number or inf for using threshold rule only)
+%       - relthre (relative convergence threshold for iteration to stop)
+%
+%   For more information on standard OPT fields
+%   see also defopt
+% 
+% OUTPUT: struct with the following fields:
+% -W: matrix of coefficient vectors of rls estimator for each class
+% -C: matrix of coefficient vectors of rls estimator for each class
+% -X: empty matrix
+
+lambda = opt.singlelambda(opt.paramsel.lambdas);
+
+n = size(y,1);
+
+% load in parameters alpha
+if isprop(opt,'ISTAalpha')
+    ISTAalpha=opt.ISTAalpha;
+    if ISTAalpha <= 0 || ISTAalpha > 1
+        error('Invalid alpha');
+    end
+else
+    if opt.verbose
+            fprintf('\t...alpha not defined. Use default value alpha=1 for LASSO\n');
+            ISTAalpha=1;
+    end
+end
+
+% load in number of iterations or relative tolerance
+Niter=inf;
+relthre=1e-4;
+if isprop(opt, 'ISTAniter')
+    Niter=opt.ISTAniter;
+end
+if isprop(opt, 'ISTArelthre')
+    relthre=opt.ISTArelthre;
+end
+
+% check if matrices XtX and Xty have been previously computed during
+% parameter selection
+if isfield(opt.paramsel,'XtX');
+    XtX = opt.paramsel.XtX;
+else
+    XtX = X'*X; % d x d matrix.
+end
+if isfield(opt.paramsel,'XtX');
+    Xty = opt.paramsel.Xty;
+else
+    Xty = X'*y; % d x T matrix.
+end
+
+
+%redo OLR based on non-sparsy components
+
+
+w = rls_fista_driver( XtX, Xty, n, lambda,ISTAalpha,Niter,relthre,opt.verbose);
+
+cfr.IndexFlag = ~~(w);
+% Xs=X(:,~~w);
+% cfr.Wr=zeros(size(w));
+% cfr.Wr(~~w) = rls_primal_driver(Xs'*Xs,Xs'*y,n,0);
+cfr.W = w;
+cfr.C = [];
+cfr.X = [];
+%plot(w)
+

gurls/optimizers/rls_ista.m

@@ -0,0 +1,79 @@
+function [cfr] = rls_ista (X, y, opt)
+
+% rls_ista(X,y,opt)
+% computes a classifier for elastic nerwork using ISTA.
+% The regularization parameter is set to the one found in opt.paramsel.
+% In case of multiclass problems, the regularizers need to be combined with the opt.singlelambda function.
+%
+% INPUTS:
+% -OPT: struct of options with the following fields:
+%   fields that need to be set through previous gurls tasks:
+%		- paramsel.lambdas (set by the paramsel_* routines)
+%   fields with default values set through the defopt function:
+%		- singlelambda
+%   fields with default values set through the defopt function:
+%       - ISTAalpha (paramters for balancing l1-norm and l-2 norm. 1 for LASSO and 0 for ridge)
+%       - niter (maximun number for iteration. Set to either negative number or inf for using threshold rule only)
+%       - relthre (relative convergence threshold for iteration to stop)
+%
+%   For more information on standard OPT fields
+%   see also defopt
+% 
+% OUTPUT: struct with the following fields:
+% -W: matrix of coefficient vectors of rls estimator for each class
+% -C: matrix of coefficient vectors of rls estimator for each class
+% -X: empty matrix
+
+lambda = opt.singlelambda(opt.paramsel.lambdas);
+
+n = size(y,1);
+
+% load in parameters alpha
+if isprop(opt,'ISTAalpha')
+    ISTAalpha=opt.ISTAalpha;
+    if ISTAalpha <= 0 || ISTAalpha > 1
+        error('Invalid alpha');
+    end
+else
+    if opt.verbose
+            fprintf('\t...alpha not defined. Use default value alpha=1 for LASSO\n');
+            ISTAalpha=1;
+    end
+end
+
+% load in number of iterations or relative tolerance
+Niter=inf;
+relthre=1e-4;
+if isprop(opt, 'ISTAniter')
+    Niter=opt.ISTAniter;
+end
+if isprop(opt, 'ISTArelthre')
+    relthre=opt.ISTArelthre;
+end
+
+% check if matrices XtX and Xty have been previously computed during
+% parameter selection
+if isfield(opt.paramsel,'XtX');
+    XtX = opt.paramsel.XtX;
+else
+    XtX = X'*X; % d x d matrix.
+end
+if isfield(opt.paramsel,'XtX');
+    Xty = opt.paramsel.Xty;
+else
+    Xty = X'*y; % d x T matrix.
+end
+
+%redo OLR based on non-sparsy components
+
+w = rls_ista_driver( XtX, Xty, n, lambda,ISTAalpha,Niter,relthre,opt.verbose);
+
+cfr.IndexFlag = ~~(w);
+% Xs=X(:,~~w);
+% cfr.Wr=zeros(size(w));
+% cfr.Wr(~~w) = rls_primal_driver(Xs'*Xs,Xs'*y,n,0);
+cfr.W = w;
+cfr.C = [];
+cfr.X = [];
+%plot(w)
+

gurls/paramsel/paramsel_hoista.m

@@ -0,0 +1,134 @@
+function vout = paramsel_hoista(X, y, opt)
+% paramsel_hoista(X,Y,OPT)
+% Performs parameter selection ISTA for elastic net is used.
+% The hold-out approach is used. 
+% The performance measure specified by opt.hoperf is maximized.
+%
+% INPUTS:
+% -OPT: struct of options with the following fields:
+%   fields that need to be set through previous gurls tasks:
+%		- split (set by the split_* routine)
+%   fields with default values set through the defopt function:
+%		- nlambda
+%		- hoperf
+%       - nholdouts
+%
+%   For more information on standard OPT fields
+%   see also defopt
+% 
+% OUTPUTS: structure with the following fields:
+% -lambdas_round: cell array (opt.nholdoutsX1). For each split a cell contains the 
+%       values of the regularization parameter lambda minimizing the 
+%       validation error for each class.
+% -perf: cell array (opt.nholdouts). For each split a cell contains a matrix 
+%       with the validation error for each lambda guess and for each class
+% -guesses: cell array (opt.nholdoutsX1). For each split a cell contains an 
+%       array of guesses for the regularization parameter lambda
+% -lambdas: mean of the optimal lambdas across splits
+
+if isprop(opt,'paramsel')
+    vout = opt.paramsel; % lets not overwrite existing parameters.
+    % unless they have the same name
+
+    if isfield(opt.paramsel,'perf')
+        vout = rmfield(vout,'perf');
+    end
+    if isfield(opt.paramsel,'guesses')
+        vout = rmfield(vout,'guesses');
+    end
+else
+    opt.newprop('paramsel', struct());
+end
+
+% load in number of iterations or relative tolerance
+Niter=inf;
+relthre=1e-4;
+if isprop(opt, 'ISTAniter')
+    Niter=opt.ISTAniter;
+end
+if isprop(opt, 'ISTArelthre')
+    relthre=opt.ISTArelthre;
+end
+% load in alpha for elastic net
+if isprop(opt,'ISTAalpha')
+    ISTAalpha=opt.ISTAalpha;
+    if ISTAalpha <= 0 || ISTAalpha > 1
+        error('Invalid alpha');
+    end
+else
+    if opt.verbose
+            fprintf('\t...alpha not defined. Use default value alpha=1 for LASSO\n');
+            ISTAalpha=1;
+    end
+end
+
+
+Xtytot = X'*y;
+
+d = size(X,2);
+T = size(y,2);
+
+% Verify the parameter selection constraint is defined proper
+opt.ISTASparsitylvl=ceil(opt.ISTASparsitylvl);
+if (opt.ISTASAccuReq>0&&opt.ISTASparsitylvl<d)
+    error('Please define parameter selection constraint correctly')
+end
+
+for nh = 1:opt.nholdouts
+	if iscell(opt.split)
+		tr = opt.split{nh}.tr;
+		va = opt.split{nh}.va;
+	else	
+		tr = opt.split.tr;
+		va = opt.split.va;
+	end	
+
+	n = length(tr);
+
+	XtXtr = (X(tr,:))'*X(tr,:);
+	Xtytr =(X(tr,:))'*y(tr,:);
+	tot = opt.nlambda;
+
+
+    L = 2*max(abs(Xtytot),[],2)/(n*ISTAalpha);
+	guesses = n*paramsel_lambdaguesses(L, min(n,d), n, opt);
+
+    if ~isprop(opt, 'rls')
+        opt.newprop('rls', struct());
+    end
+
+	ap = zeros(tot,T);
+    apt=ap;
+    Sparsity = zeros(tot,1);
+    for i = 1:tot
+        W = rls_ista_driver( XtXtr, Xtytr, n, guesses(i),ISTAalpha,Niter,relthre,0);
+        Sparsity(i) = sum(~~W);
+        opt.rls.W=W;
+        opt.newprop('pred', pred_primal(X(va,:),y(va,:),opt));
+        opt.newprop('perf', opt.hoperf(X(va,:),y(va,:),opt));
+        for t = 1:T
+            ap(i,t) = opt.perf.forho(t);
+        end	
+    end	
+    if opt.ISTASparsitylvl<d
+        apt((Sparsity>opt.ISTASparsitylvl),:)=-1;
+        [~, idx] = max(apt,[],1);	
+    else
+        Sparsity(min(ap,[],2)<opt.ISTASAccuReq)=inf;
+        [~, idx] = min(Sparsity);
+    end
+
+    if (max(apt,[],1)<0)||isinf(min(Sparsity))
+        warning('...Cannot find sych lambda that fit the requirement. Use default requirment instead')
+        [~, idx] = max(ap,[],1);
+    end
+    vout.lambdas_round{nh} = guesses(idx);
+	vout.perf{nh} = ap;
+	vout.guesses{nh} = guesses;
+end	
+if numel(vout.lambdas_round) > 1
+	lambdas = cell2mat(vout.lambdas_round');
+	vout.lambdas = mean(lambdas);
+else
+	vout.lambdas = vout.lambdas_round{1};
+end

gurls/utils/rls_fista_driver.m

@@ -0,0 +1,62 @@
+function w=rls_fista_driver( XtX, Xty, n, lambda,inst_alpha,Niter,relthre,verbose)
+% rls_ista_driver( XtX, Xty, n, lambda,alpha,Niter,reltol,verbose)
+% Utility function used by rls_ista
+% 
+% INPUTS:
+% -XtX: symmetric dxd square matrix
+% -Xty: dxT matrix
+% -n: number of training samples
+% -lambda: regularization parameter
+% -alpha: l1 norm l2 norm weight parameter
+% -Niter: maximum number of iteration
+% -retol: relative tolernce for convergence
+% -verbose: output option
+%
+% OUTPUTS:
+% -W: matrix of coefficient vector for linear RLS classifier
+
+    gamma = n/(2*eigs(XtX,1)+lambda*(1-inst_alpha)); %Step size
+
+    % make sure Niter, retol are proper
+    if (Niter-round(Niter))~=0 && ~isinf(Niter)
+        if verbose
+        fprintf('\t...Interation number must be integer. Rounding opt.paramsel.niter\n');
+        end
+        Niter = ceil(Niter);
+    end
+    if Niter<=0 && relthre<0
+        if verbose
+            fprintf(['\t...Unvalid stopping rule for ISTA.',...
+            'Using default relative tolerance = 1e-4\n']);
+        end
+        Niter = -1;
+        relthre=1e-4;      
+    end
+
+    % start ista
+    t0=1;t1=1;
+    w2=rls_primal_driver(XtX, Xty, n, 0);
+    w2(abs(w2)<lambda*inst_alpha*gamma)=0;
+    w2=w2-sign(w2)*lambda*inst_alpha*gamma;
+    w1=0*w2;
+    %dw = mean(abs(w1-w2));
+    k=0;
+
+    while mean(abs(w1-w2))>relthre*max([mean(abs(w1)),mean(abs(w2))])
+        yk=w2+(t0-1)*(w2-w1)/t1;
+        w1=w2;
+        w2=(1-lambda*gamma*(1-inst_alpha))*yk+2*gamma*(Xty-XtX*yk)/n;
+        w2(abs(w2)<lambda*inst_alpha*gamma)=0;
+ %       plot(w2)
+        w2=w2-sign(w2)*lambda*inst_alpha*gamma;
+        t0=t1;
+        t1=(1+(1+4*t1^2)^0.5)/2;
+        k=k+1;
+        if k==Niter
+            break;
+        end
+
+    end
+
+    w=w2;
+end