% DSVRG implements distributed SVRG algorithm for Example 26.3 where
% a logistic loss is minimized for binary classification. 
% The key ingredient of the algorithm is the stochastic variance reduced
% gradient (SVRG) iteration by Johnson and Zhang (2013) that is applied
% in a distributed setting.
% Input:
% X: Training samples.
% y: labels of training samples.
% mu: denotes parameter \mu.
% a: stepsize for SVRG iteration.
% T: length of inner loop.
% N: number of nodes.
% st: initial randon states.
% S: length of outer loop. 
% output:
% ws: optimized weights.
% f: profile of softmax loss.
% rtm: profile of rate of misclassification.
% ==============================================
% To run the code, prepare data sets as follows:
% Example 1:
% load X2; load X7; load T2; load T7; X = [X2 X7]; 
% y = [ones(2000,1); -ones(2000,1)];
% [ws,f,rtm] = DSVRG_rlog(X,y,T2,T7,0,0.3,20,200,17,66);
% Example 2:
% load X3; load X5; load T3; load T5; X = [X3 X5]; 
% y = [ones(2000,1); -ones(2000,1)];
% [ws,f,rtm] = DSVRG_rlog(X,y,T3,T5,0,0.22,20,200,17,101);
function [ws,f,rtm] = DSVRG_rlog(X,y,Xte1,Xte2,mu,a,T,N,st,S)
format long
format compact
[N1,P] = size(X);
Xh = [X; ones(1,P)];
t1 = size(Xte1,2);
t2 = size(Xte2,2);
A1 = [Xte1; ones(1,t1)];
A2 = [Xte2; ones(1,t2)];
w0 = zeros(N1+1,1);
wt = w0;
f0 = f_rlog(wt,Xh,y,0);
f = f0;
y1 = wt'*A1;
y2 = wt'*A2;
C = zeros(2,2);
C(1,1) = length(find(y1 > 0));
C(2,1) = t1 - C(1,1);
C(2,2) = length(find(y2 < 0));
C(1,2) = t2 - C(2,2);
rtm = (1-((C(1,1)+C(2,2))/sum(sum(C))))*100;
current = [0 f rtm]
z = floor(P/N);
for s = 1:S
  gt = g_rlog(wt,Xh,y,mu);
  we = 0;
  for k = 1:N
      wk = wt;
      Pk = Xh(:,(k-1)*z+1:k*z);
      yk = y((k-1)*z+1:k*z);
      for t = 1:T
        rand('state',st+(k-1)*T+t)
        ind2 = randperm(z);
        i = ind2(1);
        xi = Pk(:,i);
        yi = y(i);
        gik = g_rlog(wk,xi,yi,mu);
        git = g_rlog(wt,xi,yi,mu);
        gk = gik - git + gt;
        wk = wk - a*gk;
      end
      we = we + wk - wt;
  end
  wt = wt + we/N;
  y1 = wt'*A1;
  y2 = wt'*A2;
  C(1,1) = length(find(y1 > 0));
  C(2,1) = t1 - C(1,1);
  C(2,2) = length(find(y2 < 0));
  C(1,2) = t2 - C(2,2);
  rk = (1-((C(1,1)+C(2,2))/sum(sum(C))))*100;
  rtm = [rtm rk];
  fs = f_rlog(wt,Xh,y,0);
  f = [f fs];
  current = [s fs rk]
end
ws = wt;
disp('confusion matrix')
C
disp('rate of misclassification in percentage:')
rk
figure(1)
set(gca,'fontsize',14,'fontname','times')
semilogy(0:1:S,f,'b-','linewidth',1.5)
xlabel('Iterations \its','fontsize',15,'fontname','times')
ylabel('Objective  {\itE_L}{\rm(}{\itw_k}{\rm)}','fontsize',15,'fontname','times')
axis([0 S 3e-2 1])
grid