run_rsg_svm_bank.m

close all; clear all; clc
%% Input Parameters
    
ra = rng(105,'twister'); % Setting some random seed for reproducibility
step = 0; % Setting type of stepsize; value of 0 mean min{1/L, 1/(d^2N)^2/3}, 1 means min{1/L, c/sqrt(N)}, 2 means 1\L, 3 means 1\LK with eta = 1

initial = 5; % Initial point

lambda = 0.01; % Constant Lagrange multiplier
NList = [10, 100, 250, 500, 1000, 1850, 3000, 5000]; % Iterations list
T = 50; % Replications

%% Loading Dataset

load('bank_dataset.mat');

u = dataset(:,1:4);
v = dataset(:,5);
v(v == 0) = -1;


d = size(u,2); % Dimension of the problem
dataset_len = size(u,1); % Dataset length
reduce_dim = 0; % Reduce dimension by SVD?

x_1 = initial.*unifrnd(-1,1,d,1); % Initial point of d dimension


% Performing min max normalization
u_min = min(u);
u_max = max(u);
diff = u_max - u_min;

for i = 1:dataset_len
    u(i,:) = (u(i,:) - u_min)./diff;
end

% Randomly shuffle dataset
dataset = [u, v];
dataset = dataset(randperm(size(dataset,1)),:);

u = dataset(:,1:d);
v = dataset(:,d+1);

%% Performing SVD to reduce dimension

if reduce_dim == 1
    original_norm = norm(u,'fro');
    [U,S,V,flag] = svds(u,min(size(u)));

    % Extract singular values
    singvals = diag(S);
    originalss = sum(singvals.^2);

    % Find out where to truncate the U, S, V matrices to get 90% accuracy
    tmp = 0;
    for i = 1:d
        tmp = tmp + singvals(i).^2;
        if tmp > 0.9*originalss
            newd = i;
            break
        end
    end

    reconstructed = U(:,1:newd)*S(1:newd,1:newd);
    reconstructed_norm = norm(reconstructed,'fro');
    err = abs(reconstructed_norm - original_norm)*100/original_norm;

    str = sprintf('Reduced dimension is %d and frobenious norm error is %f',newd,err); 
    disp(str);

    u = reconstructed;
    d = newd;
end

%% Splitting training and testing dataset  
training_len = round(0.6*dataset_len);
test_len = round(0.4*dataset_len);

training_u = u(1:training_len,:);
training_v = v(1:training_len);

test_u = u(training_len+1:end,:);
test_v = v(training_len+1:end);

%% Paramater Estimation L and sigma

N_0 = 200; % Using N_0 i.i.d. samples 
approx_samples = 200;

gradient_var = zeros(N_0,1);
B_list = zeros(N_0,1);
hessian_norm = zeros(N_0,1);
third_der_mag = zeros(N_0,1);
f_sigma = zeros(N_0,1);


x = randi(10).*unifrnd(-1,1,d,N_0); % Generating random points

f_x_1 = 0;
iid_samples_index = randi([1 training_len],N_0,1);
u_iid = training_u(iid_samples_index,:);
v_iid = training_v(iid_samples_index,:);

for i = 1:N_0 % looping over N_0 random points x
    
    hessian = zeros(d,d);
    gradient_norm = zeros(approx_samples,1);
    grad_l2norm = 0;
    third_der = zeros(approx_samples,1);
    f_x = zeros(approx_samples,1);
   
    
    for j = 1:approx_samples % looping over N_0 iid samples (u,v)
       gradient = SVM( x(:,i), u_iid(j,:)', v_iid(j), lambda, 0, 'ub', 1 ); % Calling SFO
       gradient_norm(j) = norm(gradient,2)^2;
       grad_l2norm = grad_l2norm + norm(gradient,2);

       
       hessian = hessian + SVM( x(:,i), u_iid(j,:)', v_iid(j), lambda, 0, 0, 2 ); % Getting hessian
       
       third_der(j) = SVM( x(:,i), u_iid(j,:)', v_iid(j), lambda, 0, 0, 3 ); % Getting third derivative at x
       
       f_x(j) = SVM( x(:,i), u_iid(i,:)', v_iid(i), lambda, 0, 0, 0); % Calling SZO
    end
    
    gradient_var(i) = var(gradient_norm);
    B_list(i) = grad_l2norm/N_0;
    
    hessian = hessian./approx_samples;
    hessian_norm(i) = norm(hessian,2);
    
    third_der_mag(i) = max(third_der);
    
    f_sigma(i) = var(f_x);
    
    
    f_x_1 = f_x_1 + SVM( x_1, u_iid(i,:)', v_iid(i), lambda, 0, 0, 0);
   
end

% Estimate of f(x1)
f_x_1 = f_x_1/N_0;

L = max(hessian_norm); % Estimate of L
sigma = max(gradient_var); % Estimate of sigma
f_sigma2 = max(f_sigma); % Estimate of noise in f
B = max(B_list); % Estimation of B

alpha_0 = max(third_der_mag); % Estimate of third derivative upper bound
alpha_1 = (2 + lambda*d*3^2)^2; % Estimate of E[f(x + eta*delta)^2] upper bound

%% Declaring variables to store results

approx_samples = test_len; % Samples to estimate f or gradient, % To find quality of solution at x_R

ub = [NList', zeros(length(NList),4)]; % Itr, mean_sng, std. err sng, classification acc. mean, classification acc. std. err
gs = [NList', zeros(length(NList),4)];
spsa = [NList', zeros(length(NList),4)];

rdsa_u = [NList', zeros(length(NList),4)];
rdsa_ab = [NList', zeros(length(NList),4)];


rdsa_lex = [NList', zeros(length(NList),4)];
rdsa_perm = [NList', zeros(length(NList),4)];
rdsa_kw = [NList', zeros(length(NList),4)];

%% Saving dataset and parameter
% save(['bank_dataset']);

%% Loading dataset and parameter
% load(['bank_dataset']);

%% RSG_Unbiased algorithm
rng(ra);
algo = 'ub';
i = 1;
 disp('--------------------------------Starting Unbiased----------------------------------');
for N = NList
    str=sprintf('N = %d \n', N); disp(str);

    % Setting step size and Probability mass function
    eta = 0*ones(N,1);
    c = sqrt(2*(f_x_1)/(L*sigma));

    if step == 0
       gamma = ones(N,1)*min(1/L, c/(sqrt(N)));
    elseif step == 1
       gamma = ones(N,1).*min(1/L, c/(sqrt(N)));
    elseif step == 2
       gamma = ones(N,1).*1/L;
    elseif step == 3
       gamma = 1/L.*1./(1:N); % Decreasing step size 
    end

    P_R = gamma./sum(gamma);
    
    ub(i,2:5) = rsg_svm(algo, N, T, approx_samples, x_1, gamma, P_R, eta, lambda, training_u, training_v, test_u, test_v);
    i = i+1;
end
disp(ub(:,[2,4]));
disp('--------------------------------End Unbiased----------------------------------');

plot(NList,ub(:,2));
save(['svm_bank_ub_step_',num2str(step),'_x1_',num2str(initial),'.txt'],'ub','-ascii')

%% GS
rng(ra);
algo = 'gs';
i = 1;
sigma = f_sigma2;

% Setting constants c1 and c2
c1 = L*(d+3)^(3/2)/2;
c2 = L^2*(d+3)^(3)/2;
c3 = 2*(d+5)*(B^2 + sigma^2);


 disp('--------------------------------Starting GS----------------------------------');
for N = NList
    str=sprintf('N = %d', N); disp(str);
    
    % Setting step size, Probability mass function and eta
    eta_0 = (d)^(-1);   
    eta = eta_0/(N^(1/2)).*ones(N,1);
    c = sqrt((2*(f_x_1)*N)/(L*(c1^2*eta_0^2 + eta_0^2*c2 + N*c3)));
 
   
    if step == 0
       gamma = ones(N,1)*min(1/L, c/(sqrt(N)));
    elseif step == 1
       gamma = ones(N,1)*min(1/L, c/(N^(1/2)));
    elseif step == 2
       gamma = ones(N,1).*1/L;
    elseif step == 3
       gamma = 1/L.*1./(1:N); % Decreasing step size 
       eta = ones(N,1);
    end
    P_R = gamma./sum(gamma);
    
   gs(i,2:5) = rsg_svm(algo, N, T, approx_samples, x_1, gamma, P_R, eta, lambda, training_u, training_v, test_u, test_v);
    i = i+1;
end
disp(gs(:,[2,4]));
 disp('--------------------------------End GS----------------------------------');

plot(NList,gs(:,2));
save(['svm_bank_gs_step_',num2str(step),'_x1_',num2str(initial),'.txt'],'gs','-ascii')

%% 1SPSA
rng(ra);
algo = 'spsa';
i = 1;
sigma = f_sigma2;

% Setting constants c1 and c2
c1 = alpha_0*d^3/6;
c2 = 2*(alpha_1 + sigma^2)*d;

disp('--------------------------------Starting 1SPSA----------------------------------');
for N = NList
    str=sprintf('N = %d', N); disp(str);
    
    % Setting step size, Probability mass function and eta
    eta_0 = d^(-5/6);
    eta = eta_0/(N^(1/6)).*ones(N,1);
    c = sqrt((2*(f_x_1)*N*eta_0^2)/(L*c1^2*eta_0^6 + N*L*c2));
    
    if step == 0
       gamma = ones(N,1)*min(1/L, 1/(d^2*N)^(2/3));
    elseif step == 1
       gamma = ones(N,1)*min(1/L, c/(N^(2/3)));
    elseif step == 2
       gamma = ones(N,1).*1/L;
    elseif step == 3
       gamma = 1/L.*1./(1:N); % Decreasing step size 
       eta = ones(N,1);
    end
    
    P_R = gamma./sum(gamma);
    
    spsa(i,2:5) = rsg_svm(algo, N, T, approx_samples, x_1, gamma, P_R, eta, lambda, training_u, training_v, test_u, test_v);
    i = i+1;
end
disp(spsa(:,[2,4]));
disp('--------------------------------End 1SPSA----------------------------------');

plot(NList,spsa(:,2));
save(['svm_bank_spsa_step_',num2str(step),'_x1_',num2str(initial),'.txt'],'spsa','-ascii')

%% 1RDSA_Uniform
rng(ra);
algo = 'rdsa_u';
i = 1;
u = 1; % Setting 1RDSA_Uniform distribution parameter

% Setting constants c1 and c2
c1 = alpha_0*d^3*u^4/(6*1);
c2 = (alpha_1 + sigma^2)*(u)^2*d*2/1; 

 disp('--------------------------------Starting 1RDSA_Uniform----------------------------------');
for N = NList
    str=sprintf('N = %d \n', N); disp(str);
    
    % Setting step size, Probability mass function and eta
    eta_0 = d^(-5/6);
    eta = eta_0/(N^(1/6)).*ones(N,1);
    c = sqrt((2*(f_x_1)*N*eta_0^2)/(L*c1^2*eta_0^6 + N*L*c2));
    
    if step == 0
       gamma = ones(N,1)*min(1/L, 1/(d^2*N)^(2/3));
    elseif step == 1
       gamma = ones(N,1)*min(1/L, c/(N^(2/3)));
    elseif step == 2
       gamma = ones(N,1).*1/L;
    elseif step == 3
       gamma = 1/L.*1./(1:N); % Decreasing step size 
       eta = ones(N,1);
    end
    
    P_R = gamma./sum(gamma);
    
    rdsa_u(i,2:5) = rsg_svm(algo, N, T, approx_samples, x_1, gamma, P_R, eta, lambda, training_u, training_v, test_u, test_v);
    i = i+1;
end
disp(rdsa_u(:,[2,4]));
disp('--------------------------------End 1RDSA_Uniform----------------------------------');

plot(NList,rdsa_u(:,2));
save(['svm_bank_rdsa_u_step_',num2str(step),'_x1_',num2str(initial),'.txt'],'rdsa_u','-ascii')

%% 1RDSA_Asymber
rng(ra);
algo = 'rdsa_ab';
i = 1;

ab = 0.0001; % Setting 1RDSA_Asymber distribution parameter

% Setting constants c1 and c2
c1 = alpha_0*d^3*(1 + ab)^4/6;
c2 = (alpha_1 + sigma^2)*(1 + ab)*d*2;

disp('--------------------------------Starting 1RDSA_Asymber----------------------------------');
for N = NList+1
    str=sprintf('N = %d \n', N); disp(str);
    
    % Setting step size, Probability mass function and eta
    eta_0 = d^(-5/6);
    eta = eta_0/(N^(1/6)).*ones(N,1);
    c = sqrt((2*(f_x_1)*N*eta_0^2)/(L*c1^2*eta_0^6 + N*L*c2));
    
    if step == 0
       gamma = ones(N,1)*min(1/L, 1/(d^2*N)^(2/3));
    elseif step == 1
       gamma = ones(N,1)*min(1/L, c/(N^(2/3)));
    elseif step == 2
       gamma = ones(N,1).*1/L;
    elseif step == 3
       gamma = 1/L.*1./(1:N); % Decreasing step size 
       eta = ones(N,1);
    end
    
    P_R = gamma./sum(gamma);
    
    rdsa_ab(i,2:5) = rsg_svm(algo, N, T, approx_samples, x_1, gamma, P_R, eta, lambda, training_u, training_v, test_u, test_v);
    i = i+1;
end
 disp(rdsa_ab(:,[2,4]));
 disp(' -------------------------------End 1RDSA_Asymber----------------------------------');
 
plot(NList,rdsa_ab(:,2));
save(['svm_bank_rdsa_ab_step_',num2str(step),'_x1_',num2str(initial),'.txt'],'rdsa_ab','-ascii')

%% 1RDSA_Perm_DP
algo = 'rdsa_perm';
i = 1;

% Setting constants c1 and c2
c1 = alpha_0*d^3/6;
c2 = (alpha_1 + sigma^2)*d/2;

disp('--------------------------------Starting 1RDSA_Perm_DP----------------------------------');
for N = NList
    str=sprintf('N = %d \n', N); disp(str);
    
    % Setting step size, Probability mass function and eta
    eta_0 = d^(-5/6);
    eta = eta_0/(N^(1/6)).*ones(N,1);
    c = sqrt((2*(f_x_1)*N*eta_0^2)/(L*c1^2*eta_0^6 + N*L*c2));
    
    if step == 0
       gamma = ones(N,1)*min(1/L, 1/(d^2*N)^(2/3));
    elseif step == 1
       gamma = ones(N,1)*min(1/L, c/(N^(2/3)));
    elseif step == 2
       gamma = ones(N,1).*1/L;
    elseif step == 3
       gamma = 1/L.*1./(1:N); % Decreasing step size 
       eta = ones(N,1);
    end
    
    P_R = gamma./sum(gamma);
    
    rdsa_perm(i,2:5) = rsg_svm(algo, N, T, approx_samples, x_1, gamma, P_R, eta, lambda, training_u, training_v, test_u, test_v);
    i = i+1;
end
 disp(rdsa_perm(:,[2,4]));
 disp('--------------------------------End 1RDSA_Perm_DP----------------------------------');
plot(NList,rdsa_perm(:,2));
save(['svm_bank_rdsa_perm_step_',num2str(step),'_x1_',num2str(initial),'.txt'],'rdsa_perm','-ascii')

 %% 1RDSA_KW_DP
rng(ra);
algo = 'rdsa_kw';
i = 1;

% Setting constants c1 and c2
c1 = alpha_0*d^3/6;
c2 = (alpha_1 + sigma^2)*d/2;

 disp('--------------------------------Starting 1RDSA_KW_DP----------------------------------');
for N = NList
    str=sprintf('N = %d \n', N); disp(str);
    
    % Setting step size, Probability mass function and eta
    eta_0 = d^(-5/6);
    eta = eta_0/(N^(1/6)).*ones(N,1);
    c = sqrt((2*(f_x_1)*N*eta_0^2)/(L*c1^2*eta_0^6 + N*L*c2));
    
    if step == 0
       gamma = ones(N,1)*min(1/L, 1/(d^2*N)^(2/3));
    elseif step == 1
       gamma = ones(N,1)*min(1/L, c/(N^(2/3)));
    elseif step == 2
       gamma = ones(N,1).*1/L;
    elseif step == 3
       gamma = 1/L.*1./(1:N); % Decreasing step size 
       eta = ones(N,1);
    end
    
    P_R = gamma./sum(gamma);
    
    rdsa_kw(i,2:5) = rsg_svm(algo, N, T, approx_samples, x_1, gamma, P_R, eta, lambda, training_u, training_v, test_u, test_v);
    i = i+1;
end
 disp(rdsa_kw(:,[2,4]));
 disp('--------------------------------End 1RDSA_KW_DP----------------------------------');
plot(NList,rdsa_kw(:,2));
save(['svm_bank_rdsa_kw_step_',num2str(step),'_x1_',num2str(initial),'.txt'],'rdsa_kw','-ascii')

%% 1RDSA_Lex_DP
% rng(ra);
% algo = 'rdsa_lex';
% i = 1;
% 
% % Setting constants c1 and c2
% c1 = alpha_0*d^3*3^(d-1);
% c2 = (alpha_1 + sigma^2)*d*2;
% 
% disp('--------------------------------Starting 1RDSA_Lex_DP----------------------------------');
% for N = NList
%     str=sprintf('N = %d \n', N); disp(str);
%     
%     % Setting step size, Probability mass function and eta
%     eta_0 = d^(-5/6);
%     eta = eta_0/(N^(1/6)).*ones(N,1);
%     c = sqrt((2*(f_x_1)*N*eta_0^2)/(L*c1^2*eta_0^6 + N*L*c2));
%     
%     if step == 0
%        gamma = ones(N,1)*min(1/L, 1/(d^2*N)^(2/3));
%     elseif step == 1
%        gamma = ones(N,1)*min(1/L, c/(N^(2/3)));
%     elseif step == 2
%        gamma = ones(N,1).*1/L;
%     elseif step == 3
%        gamma = 1/L.*1./(1:N); % Decreasing step size 
%        eta = ones(N,1);
%     end
%     P_R = gamma./sum(gamma);
%     
%     rdsa_lex(i,2:5) = rsg_svm(algo, N, T, approx_samples, x_1, gamma, P_R, eta, lambda, training_u, training_v, test_u, test_v);
%     i = i+1;
% end
%  disp(rdsa_lex(:,[2,4]));
%  disp('--------------------------------End 1RDSA_Lex_DP----------------------------------');
% plot(NList,rdsa_lex(:,2));
% save(['svm_bank_rdsa_lex_step_',num2str(step),'_x1_',num2str(initial),'.txt'],'rdsa_kw','-ascii')


%% Plotting results

% Plotting SNG
figure;
plot(NList,ub(:,2));
hold on
plot(NList,gs(:,2));
plot(NList,spsa(:,2));
plot(NList,rdsa_u(:,2));
plot(NList,rdsa_ab(:,2));
% plot(NList,rdsa_lex(:,2));
plot(NList,rdsa_perm(:,2));
plot(NList,rdsa_kw(:,2));
title('SNG');
% legend('Unbiased','GS','SPSA','RDSA-Unif','RDSA-AsymBer','RDSA-Lex-DP','RDSA-Perm-DP','RDSA-KW-DP');
legend('Unbiased','GS','SPSA','RDSA-Unif','RDSA-AsymBer','RDSA-Perm-DP','RDSA-KW-DP');


%% Saving results to txt file

save(['svm_bank_step_',num2str(step),'_x1_',num2str(initial)]);