% Logistic regression for binary classification by Stochastic variance reduced 
% gradient(SVRG). 
% Main reference: R. Johnson and T. Zhang, "Accelerating stochastic
% gradient descent using predictive variance reduction" in Advnces in
% Neural Information Processing Systems, pp. 315-323, 2013. 
% Example 1: 
% load X2; load X7; load T2; load T7; X = [X2 X7]; 
% y = [ones(2000,1); -ones(2000,1)];
% [ws,f,rtm] = SVRG_rlog(X,y,T2,T7,0,0.12,340,9,34);
% Example 2 
% load X3; load X5; load T3; load T5; X = [X3 X5]; 
% y = [ones(2000,1); -ones(2000,1)];
% [ws,f,rtm] = SVRG_rlog(X,y,T3,T5,0.001,0.11,340,9,55);
function [ws,f,rtm] = SVRG_rlog(X,y,Xte1,Xte2,mu,a,T,st,K)
format long
format compact
[N,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(N+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;
fprintf('iter %1i: obj =%9.2e\n',0,f);
tt = 1;
for k = 1:K
    gt = g_rlog(wt,Xh,y,mu);
    wk = wt;
    for t = 1:T
        rand('state',st+(k-1)*T+t)
        ind = randperm(P);
        it = ind(1:tt);
        xi = Xh(:,it);
        yi = y(it);
        gik = g_rlog(wk,xi,yi,mu);
        git = g_rlog(wt,xi,yi,mu);
        gk = gik - git + gt;
        wk = wk - a*gk;
    end
    wt = wk;
    fk = f_rlog(wt,Xh,y,0);
    f = [f; fk];
    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];
    fprintf('iter %1i: obj =%9.2e\n',k,fk);
end
ws = wt;
disp('Objective function at the solution point:')
fs = f(end)
disp('confusion matrix')
C
disp('rate of misclassification in percentage:')
rk
figure(1)
set(gca,'fontsize',14,'fontname','times')
semilogy(0:1:K,f,'b-','linew',1.5)
xlabel('Iterations \itk','fontsize',15,'fontname','times')
ylabel('Objective  {\itE_L}{\rm(}{\itw_k}{\rm)}','fontsize',15,'fontname','times')
axis([0 K 1.5e-2 1])
grid