%% Multiple Kernel Learning, Conic Duality, and the SMO Algorithm
% Team Member: Qiang Zhao

%% Problem Formulation
% Nowadays, classicial kernel-based classifiers are based on a single
% kernel, the paper considered a conic combinations of kernel matricesfor
% the support vector machine(SVM), and showed that the optimization of the
% coefficients of such a combination reduces to a convex optimization
% problem known as a quadratcally-constrained quadratic program(QCQP).
% As proposed in the paper, I will show a novel dual formulation of the
% QCQP as a second-order cone programming problem, and show how to exploit
% the technique of Moreau-Yosida regularization to yield a formulation to
% which sequwntial minimal optimization(SMO) techniques can be applied. The
% dataset used in this prokect is the breast cancer from the UCI repository. 

%% Proposed solution and its related optimization problem
% The primal problem for the SKM is:
%%
% $min \ \frac{1}{2}(\sum_{j=1}^m d_j\|w_j\|_2)^2 + C\sum_{n=1}^n\xi_i$
%%
% $s.t.\ y_i(\sum_{j=1}^m w_j^\mathrm{T}x_{ji} + b) \geq 1 - \xi_i, \forall i\in \{1,...,n\}$
%%
% The primal problem can be put in the following equivalent form, where
% $K_{k}=\{(u,v) \in R^{k+1}, \|u\|_2 \le v\}$ is the second-order cone of
% order k, then we have
%%
% $min \ \frac{1}{2}u^2 + C\sum_{n=1}^n\xi_i$
%%
% $w.r.t\ u\in R, t\in R^m, b\in R, \xi\in R_+^n, (w_j,t_j)\in K_{k_j}, \forall j$
%%
% $s.t.\  y_i(\sum_{j=1}^m w_j^\mathrm{T}x_{ji} + b) \geq 1 - \xi_i, \forall i$
%%
% $\qquad \  \sum_{j=1}^m d_jt_j \leq u$
%%
% The cone $K_k$ is self-dual, so the conic Lagrangian function corresponding to the problem is L = :
%%
% $L=\frac{1}{2}u^2 + C\sum_{n=1}^n\xi_i - \sum_{n=1}^n\alpha_i(y_i(\sum_{j=1}^m w_j^\mathrm{T}x_{ji} + b) - 1 + \xi_i) - \sum_{n=1}^n\beta_i\xi_i + \gamma(\sum_{j=1}^md_jt_j-u)-\sum_{j=1}^m(\lambda_j^\mathrm{T}w_j + \mu_jt_j)$ 
%%
% with $\alpha_i \in R_+, \beta_i\in R_+, \gamma\in R_+, 
% (\lambda_j,\mu_j)\in K_{k_j}$
%%
% If $w_j \neq 0$, we can take derivative to Lagrangian function, thus we
% have $\frac{\partial L}{\partial w_j}=\sum_{j=1}^m d_j\|w_j\|_2\frac{w_j}{\|w_j\|_2}-\sum_{n=1}^n\alpha_iy_ix_{ji}$
%%
% Finally, after computing derivatives with respect to the primal variables
% and setting them to zero, we readily get the dual function:
% $g(\alpha,\beta,\gamma,\lambda,\mu)=-\frac{1}{2}\gamma^2 +
% \sum_{n=1}^n\alpha_i$, defined by $\alpha_i\geq 0,\ \beta_i\geq 0,\ \gamma\geq 0,\ \|\lambda_j\|_2\leq\mu_j,d_j\gamma - \mu_j=0,\ \sum_{n=1}^n\alpha_i y_i x_{ji} + \lambda_j = 0,\ \sum_{n=1}^n\alpha_i y_i = 0,\ C-\alpha_i-\beta_i=0$
%%
% After elimination of dummy variables, we obtain the following dual
% problem:
%%
% $max -\frac{1}{2}\gamma^2 + \alpha^\mathrm{T}e$
%%
% $w.r.t\ \gamma\in R,\ \alpha\in R^n$
%%
% $s.t.\ 0\leq\alpha\leq C,\ \alpha^\mathrm{T}y=0$
%%
% $\|\sum_{n=1}^n\alpha_i y_i x_{ji}\|_2 \leq d_j\gamma,\ \forall j\in\{1,...,m\}$
%%
% In addition, the Karush-Kuhn-Tucker(KKT) optimality conditions give the
% following complementary slackness equations:
%%
% $(a)\ \alpha_i(y_i(\sum_{j=1}^mw_j^\mathrm{T}x_{ji}+b)-1+\xi_i)=0,\ \forall i$
%%
% $(b)\ (C-\alpha_i)\xi_i=0,\ \forall i$
%%
% $(c)\ \gamma(\sum_{j=1}^md_jt_j-\gamma)=0$
%% Data source
% The breast-cancer data from UCI repository and I applied the PCA on data 
% to reduce the data dimension to 2-D.
%read data
filename = 'breast-cancer_scale.csv';
delimiter = ': ';
formatSpec = '%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%[^\n\r]';
fileID = fopen(filename,'r');
dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true,  'ReturnOnError', false);
fclose(fileID);
ylabel = dataArray{:, 1};
Var2 = dataArray{:, 3};
Var3 = dataArray{:, 5};
Var4 = dataArray{:, 7};
Var5 = dataArray{:, 9};
Var6 = dataArray{:, 11};
Var7 = dataArray{:, 13};
Var8 = dataArray{:, 15};
Var9 = dataArray{:, 17};
Var10 = dataArray{:, 19};
Var11 = dataArray{:, 21};
dataset_original = [Var2,Var3,Var4,Var5,Var6,Var7,Var8,Var9,Var10,Var11];
N=683;
ylabel=ylabel-3;
%Do PCA on original data and reduce to 2-dimensional data
[coeff,score,latent,tsquared]=princomp(dataset_original);
dataset = score(:,1:2);
tic
%% SKM Primal 
% According to strong duality, the optimal value of SKM primal problem
% should be same as the SKM dual problem.
C=2;
cvx_begin
    variables w_primal(2) xi_primal(683) b_primal;
    minimize(w_primal'*w_primal/2 + C*sum(xi_primal));
    subject to
    ylabel.*(dataset*w_primal + b_primal) - 1 + xi_primal >= 0;
    xi_primal >= 0;
cvx_end
t1=toc
%% SKM Dual
e = ones(683,1);
cvx_begin
    variables gamma_dual alpha_dual(683);
    dual variables q
    maximize (-1/2*gamma_dual^2 + alpha_dual'*e);
    subject to
    0 <= alpha_dual <= C;
    q:alpha_dual'*ylabel == 0;
    ayx = 0;
    for i = 1:683
        ayx = ayx + alpha_dual(i)*ylabel(i)*dataset(i,:);
    end
    norm(ayx,2) <= gamma_dual;
cvx_end
t2=toc
%% Multiple Kernel SKM Dual
% In my implementation, I use two kernels.
kernel = ones(683*2,683);
index = zeros(1,1);
index(end+1,:) = floor(1);
index(end+1,:) = 2;
%compute the kernel matrix of the original dataset
for i=1:size(index)-1
    dataset_kernel = dataset(:,index(i,1)+1:index(i+1,1));
    K_tmp = exp(-0.5/i*multiKernel(dataset_kernel',dataset_kernel'));
    kernel(((i-1)*683+1):(i*683),:)=K_tmp;
end
diag_y = diag(ylabel);
%compute optimal dual variables
C=2;
kernel_alpha = ones(683,1);
cvx_begin
    variables kernel_gamma kernel_alpha(683,1);
    minimize(1/2*kernel_gamma^2-sum(kernel_alpha))
    subject to
    0 <= kernel_alpha <= C;
    kernel_alpha'*ylabel==0;
    kernel_alpha'*diag_y*kernel((1):(683),:)*diag_y*kernel_alpha<=kernel_gamma;
    kernel_alpha'*diag_y*kernel((684):(1366),:)*diag_y*kernel_alpha<=kernel_gamma;
cvx_end
t3=toc
%% Demonstration of results 
% Now we can visualize how kernelized SKM works compare to the original
% dual SKM classifier, we can see from the plots below that the kernalized
% SKM classifier have a descriptive contour decision boundary compare to
% the linear boundary when two classes have some samples overlapped.

%plot the decision boundary for skm dual solver
%compute the variables for dual SKM Classifier
w=dataset'*(ylabel.*alpha_dual);
b=-q;
dataset_dual = linspace(min(dataset(:,1)),max(dataset(:,1)),683);
y_dual = -(dataset_dual'*w(1)+b)/w(2);
% plot the dual SKM Classifier
figure()
x1_positive1 = dataset(find(ylabel==1),1);
x1_positive2 = dataset(find(ylabel==1),2);
x1_negative1 = dataset(find(ylabel==-1),1);
x1_negative2 = dataset(find(ylabel==-1),2);
plot(x1_positive1, x1_positive2,'.r','MarkerSize',7);
hold on
plot(x1_negative1, x1_negative2,'.b','MarkerSize',7);
axis([-2,4.5,-2.5,2]);
plot(dataset_dual,y_dual,'k','LineWidth',1);
hold off
legend('Class 1','Class -1','Location','southwest','Decision Hyperplane');
title('Dual SKM Classifier');
% plot the kernelized dual SKM Classifier
figure()
plot(x1_positive1, x1_positive2,'.r','MarkerSize',7);
hold on
plot(x1_negative1, x1_negative2,'.b','MarkerSize',7);
axis([-2,4.5,-2.5,2]);
d = 0.02;
[x1,x2] = meshgrid(min(dataset(:,1)):d:max(dataset(:,1)),min(dataset(:,2)):d:max(dataset(:,2)));
x = [x1(:),x2(:)];
kernel_test = zeros(size(x,1), size(dataset,1));
for i=1:size(index)-1,
    dataset_k = dataset(:,index(i,1) + 1:index(i+1,1));
    x_k = x(:,index(i,1) + 1:index(i+1,1));
    K = 1/2*exp(-0.5/i*multiKernel(x_k',dataset_k'));
    kernel_test = K + kernel_test;
end

prediction_test = kernel_test*(kernel_alpha.*ylabel);
contour(x1,x2,reshape(prediction_test,size(x1)),[0 0],'k');
legend('Class 1','Class -1','Location','southwest','Decision Hyperplane');
title('Kernelized SKM Classifier');
hold off

%% Summary and Conclusions
% From the above results of primal SKM and dual SKM, we can see that the
% strong duality holds and both of the optimiazation problem achieved the
% same optimal value. From the paper, the kernelized algorithm is based on
% applying sequential minimization techniques to a smoothed version of a convex non-smooth optimization problem.
% In this project, I implemented the kernelizd SKM with two kernels on brest
% dataset and visualize the results compare to the linear boundary provided by primal SKM.
% Since the breast dataset have parts of overlapped samples, so it's hard
% to judge which algorithm achieve a beter decision boundary, but we can 
% see that the contour decision boundary is more descriptive regarding to
% the distribution of the data set and thus we can say kernelized SKM have 
% better performance on seperate the non-linear data. On the other hand, the paper aims to
% save the complexity of training large scale dataset on multiple kernels,
% and from the time stamp showed on the code, the kernelized SKM algorithm
% performed almost the same computation time as the dual SKM without
% kernel, we can say that the kernelized SKM algorithm make it possible to
% learn kernels on large scale data and the good scaling with respect to
% the number of basis kernels opens up the possibility of application to
% large-scale feature selection.

