diff --git a/gurls/defopt.m b/gurls/defopt.m
index 816eb6d..40b5828 100644
--- a/gurls/defopt.m
+++ b/gurls/defopt.m
@@ -67,6 +67,10 @@
     opt.epochs = 4;
     opt.subsize = 50;
     opt.calibfile = 'foo';
+    
+    %% Subgradient
+    opt.Niter= 100;
+    opt.gammafunc = @(x) power(x,-1);
 
     %% Random features options
     opt.randfeats = struct();
diff --git a/gurls/demo/demo_svmsubgradient.m b/gurls/demo/demo_svmsubgradient.m
new file mode 100644
index 0000000..bcae52a
--- /dev/null
+++ b/gurls/demo/demo_svmsubgradient.m
@@ -0,0 +1,99 @@
+% DEMO OF THE SVM SUBGRADIENT METHOD COMPARED TO DUAL FORMULATION
+% adapted from gurlsdemo.m
+
+% Use either the 'breastcancer' or 'yeast' datasets
+dataset = 'yeast';
+if strcmp(dataset,'breastcancer')
+    load(fullfile(gurls_root, 'demo/data/breastcancer_data.mat')); 
+    ytr = ytr*2-1; yte = yte*2-1; % Convert output to {-1,1}  
+elseif strcmp(dataset,'yeast')
+    load(fullfile(gurls_root, 'demo/data/yeast_data.mat'));
+end
+
+% List the methods to use
+% pipelines = {'linsvmsub','rbfsvmsub','rbfho','rbfloo'};
+pipelines = {'linsvmsub','rbfsvmsub','rbfho'};
+
+% Number of trials to run each method
+n = 5;
+
+res_root = fullfile(gurls_root, 'demo'); % location where res files are stored
+
+for r = 1:n
+    strind = strmatch('linsvmsub',pipelines);
+    if ~isempty(strind)
+        % Gaussian kernel, SVM Subgradient, Hold Out cross validation to select lambda
+        name = [pipelines{strind} '_' num2str(r)];
+        opt = defopt(name);
+        opt.seq = {'split:ho', 'kernel:linear', 'paramsel:hodual', 'rls:svmsubgradient', 'pred:dual', 'perf:macroavg', 'perf:precrec'};
+        opt.process{1} = [2,2,2,2,0,0,0];
+        opt.process{2} = [3,3,3,3,2,2,2];
+        opt.Niter = 200;
+%         opt.gammafunc = @(x) power(x,-.5);
+        gurls(Xtr, ytr, opt, 1);
+        gurls(Xte, yte, opt, 2);
+    end
+
+    strind = strmatch('rbfsvmsub',pipelines);
+    if ~isempty(strind)
+        % Gaussian kernel, SVM Subgradient, Hold Out cross validation to select lambda and the Kernel width sigma
+        name = [pipelines{strind} '_' num2str(r)];
+        opt = defopt(name);
+        opt.seq = {'split:ho', 'paramsel:siglamho', 'kernel:rbf', 'rls:svmsubgradient', 'predkernel:traintest', 'pred:dual', 'perf:macroavg', 'perf:precrec'};
+        opt.process{1} = [2,2,2,2,0,0,0,0];
+        opt.process{2} = [3,3,3,3,2,2,2,2];
+        opt.Niter = 200;
+%         opt.gammafunc = @(x) power(x,-.5);
+        gurls(Xtr, ytr, opt, 1);
+        gurls(Xte, yte, opt, 2);
+    end
+
+    strind = strmatch('rbfho',pipelines);
+    if ~isempty(strind)
+        % Gaussian kernel, (dual formulation), Hold Out cross validation to select lambda and the Kernel width sigma
+        name = [pipelines{strind} '_' num2str(r)];
+        opt = defopt(name);
+        opt.seq = {'split:ho', 'paramsel:siglamho', 'kernel:rbf', 'rls:dual', 'predkernel:traintest', 'pred:dual', 'perf:macroavg', 'perf:precrec'};
+        opt.process{1} = [2,2,2,2,0,0,0,0];
+        opt.process{2} = [3,3,3,3,2,2,2,2];
+        gurls(Xtr, ytr, opt, 1);
+        gurls(Xte, yte, opt, 2); 
+    end
+    
+    strind = strmatch('rbfloo',pipelines);
+    if ~isempty(strind)
+        % Gaussian kernel, (dual formulation), Leave One Out cross validation to select lambda and the Kernel width sigma
+        name = [pipelines{strind} '_' num2str(r)];
+        opt = defopt(name);
+        opt.seq = {'paramsel:siglam', 'kernel:rbf', 'rls:dual', 'predkernel:traintest', 'pred:dual', 'perf:macroavg', 'perf:precrec'};
+        opt.process{1} = [2,2,2,0,0,0,0];
+        opt.process{2} = [3,3,3,2,2,2,2];
+        gurls(Xtr, ytr, opt, 1);
+        gurls(Xte, yte, opt, 2);
+    end
+end
+
+%% Visualization:
+%	- filestr: cell with names of the experiments
+%	- nRuns: number of runs for each experiment
+%	- fields: which fields of opt to display (as many plots as the elements of fields will be generated).
+%	- plotopt: a structure containing various text labels for the plots.
+
+nRuns = n*ones(1,size(pipelines,2));
+fields = {'perf.ap', 'perf.acc'};
+plotopt.titles = {'Model Comparison - Accuracy', 'Model Comparison - Precision'};
+
+% Label based on dataset used
+if strcmp(dataset,'breastcancer')
+    plotopt.class_names = {'Positive'};
+elseif strcmp(dataset,'yeast')
+    plotopt.class_names = {'CYT', 'NUC', 'MIT', 'ME3'};
+end
+
+% Generates "per-class" plots
+summary_plot(pipelines, fields, nRuns, plotopt, res_root)
+% Generates "global" plots
+summary_overall_plot(pipelines, fields, nRuns, plotopt, res_root)
+% Plots times taken by each step of the pipeline for performance reference.
+plot_times(pipelines, nRuns, res_root)
+
diff --git a/gurls/gurls_defopt.m b/gurls/gurls_defopt.m
index 97c80d5..c0eab29 100644
--- a/gurls/gurls_defopt.m
+++ b/gurls/gurls_defopt.m
@@ -62,6 +62,10 @@
     opt.newprop( 'epochs', 4);
     opt.newprop( 'subsize', 50);
     opt.newprop( 'calibfile', 'foo');
+    
+    %% Subgradient
+    opt.newprop( 'Niter', 100);
+    opt.newprop( 'gammafunc', @(x) power(x,-1));
 
     %% Random features options
     opt.newprop( 'randfeats', struct());
diff --git a/gurls/optimizers/rls_svmsubgradient.m b/gurls/optimizers/rls_svmsubgradient.m
new file mode 100644
index 0000000..c7f777c
--- /dev/null
+++ b/gurls/optimizers/rls_svmsubgradient.m
@@ -0,0 +1,67 @@
+function [rls] = rls_svmsubgradient(X,y, opt)
+
+% rls_svmsubgradient(X, y, opt)
+% computes the regression function for svm subgradient method.
+% The regularization parameter is set to the one found in opt.paramsel (set by the paramsel_* routines)
+
+% 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
+%       - Niter
+%       - gammafunc - the function that maps niter to step factor, gamma
+%                     the default is 1/niter --> [1, 1/2, 1/3, ... 1/Niter]
+% 
+%   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: empty matrix
+% -X: empty matrix
+
+lambda = opt.singlelambda(opt.paramsel.lambdas);
+Niter = opt.Niter;
+
+[n,T] = size(y);
+
+if isfield(opt.paramsel,'niter');
+    niter = max(1,ceil(opt.paramsel.niter));
+else
+    niter = 1;
+end
+
+if isfield(opt.paramsel,'f0') && niter <= Niter;
+    alpha = opt.paramsel.f0;
+else
+    niter = 1;
+    alpha=zeros(n,T);
+end
+
+iter = niter:(Niter);
+gamma = opt.gammafunc(iter);
+
+for i = iter;
+    yhat = opt.kernel.K*alpha;
+    yxyhat = y.*yhat;
+    indic = (yxyhat <= 1);
+    alpha = alpha*(1-2*gamma(i)*lambda) + (gamma(i)/n)*(y.*indic);
+end
+
+opt.paramsel.f0 = alpha;
+opt.paramsel.niter = Niter;
+
+if strcmp(opt.kernel.type, 'linear')
+	rls.W = X'*alpha;
+	rls.C = [];
+	rls.X = [];
+else
+	rls.W = [];
+    rls.C = alpha;
+	rls.X = X;
+end
+
+
+