%% ECE 602 Project: Multiple Kernel Learning, Conic Duality, and the SMO Algorithm

%% Problem Formulation
% Classical kernel based classifiers are based on a single kernel, but in
% practice classifiers based on multiple kernels are more commonly used.
% The support vector machine (SVM) based on a conic combination of kernels 
% is defined as support kernel machine (SKM). The optimization of the 
% coefficients of such a combination reduces to a convex optimization problem 
% known as a quadratically constrained quadratic program (QCQP). In this project, 
% for the convenience to plot data, we will generate a dataset with two attributes. 
% The whole dataset is divided into two classes with each data point being labeled 
% by 1 or -1. We will first consider the primal and dual problems of SVM, and 
% then the kernelized dual problem will be solved for comparison.

%% Solution 
% In linear SVM, the objective is 
% 
% $$minimize\quad\frac{1}{2}\parallel \omega\parallel_2^2+C\sum_{i=1}^{n}\xi_i$$
% 
% $$subject\ to \quad y_i(\sum_{j}\omega_j^Tx_{ji}+b)\geq 1-\xi_i \quad\forall i\in {1,\cdots, n}$$
% 
% In SVM, we encourage the sparsity of the vector $\omega$ at the level of
% blocks. The primal problem for the SVM is defined as
% 
% $$minimize\quad\frac{1}{2}(\sum_jd_j\parallel\omega\parallel_2)^2+C\sum_i\xi_i$$
% 
% $$subject\ to\quad y_i(\sum_j\omega_j^Tx_{ji}+b)\geq 1-\xi_i, \quad\xi_i\geq 0$$
% 
% The Lagrangian is
% 
% $$L=\frac{1}{2}(\sum_jd_j\parallel \omega\parallel_2)^2+C\sum_i\xi_i-\sum_i\alpha_i[y_i(\sum_j^Tx_{ji}+b)-1+\xi_i]-\sum_i\beta_i\xi_i$$
% 
% When $\omega_j\neq 0$, $L$ is differentiable, we have
% 
% $$\frac{\partial L}{\partial\omega_j}=ud_j\frac{\omega_j}{\parallel\omega_j\parallel_2}-\sum_i\alpha_iy_ix_{ji}$$
% 
% where $u=\sum_jd_j\parallel\omega_j\parallel_2$.
% At the minimum, $\frac{\partial L}{\partial\omega_j}=0$, we have
% 
% $$ud_j\frac{\omega_j}{\parallel\omega_j\parallel_2}=\sum_i\alpha_iy_ix_{ji}$$
% 
% Multiply both sides by $\omega_j^T$, we have
% 
% $$ud_j\parallel\omega_j\parallel_2=\omega_j^T\sum_i\alpha_iy_ix_{ji}$$
% 
% Take the $l_2$ norm of both sides, we have
% 
% $$ud_j=\parallel\sum_i\alpha_iy_ix_{ji}\parallel_2$$
% 
% Thus, the Lagrange dual function is
% 
% $$g(\alpha, \beta, \gamma)=-\frac{1}{2}\gamma^2+\sum_i\alpha_i$$
% 
% Lagrange dual problem becomes
% 
% $$maximize\quad -\frac{1}{2}\gamma^2+e^T\alpha$$
% 
% $$subject\ to\quad 0\leq\alpha\leq C,\quad\alpha^Ty=0,\quad\parallel\sum_i\alpha_iy_ix_{ji}\parallel_2\leq d_j\gamma$$
% 
% We now kernelize the problem using kernel function and get kernelized dual problem.
% 
% $$minimize\quad\frac{1}{2}\gamma^2-e^T\alpha$$
% 
% $$subject\ to\quad 0\leq\alpha\leq C,\quad\alpha^Ty=0,\quad(\alpha^TD(y)K_jD(y)\alpha)^2\leq\gamma d_j, \forall j$$

%% Data Generation

n = 200;
m = 2;

% Generate 100 points uniformly distributed in the unit disk
r1 = sqrt(rand(n/2,1)); % Radius
t1 = 2*pi*rand(n/2,1);  % Angle
data1 = [r1.*cos(t1), r1.*sin(t1)]; % Points

% Generate 100 points uniformly distributed in the annulus
r2 = sqrt(3*rand(n/2,1)+1); % Radius
t2 = 2*pi*rand(n/2,1);      % Angle
data2 = [r2.*cos(t2), r2.*sin(t2)]; % points

% Put the data in one matrix, and make a vector of classifications
X = [data1;data2];
Y = ones(n,1);
Y(1:n/2) = -1;

rp = randperm(n);
X = X(rp,:);
Y = Y(rp,:);

%% Primal Problem of SKM
% The optimal values of the primal and dual problems should be equal according
% to strong duality. We can use this property to verify the correctness of the 
% dual probem.
C = 1;

cvx_begin
    variables w(m) xi(n) b;  % w: 34*1, xi: 300*1, b: 1*1
%     dual variable alpha;
    minimize(1/2*w'*w+C*sum(xi));
    subject to
        Y.*(X*w+b)-1+xi >= 0;
        xi >= 0;
cvx_end

%% Primal Problem of SKM Which Encourages the Group Sparsity

cvx_begin
    variables w2(m) xi2(n) b2;  % w: 34*1, xi: 300*1, b: 1*1
    SUM = norm(w2(1),2);
    for j=2:m
        SUM = SUM+norm(w2(j),2);
    end
    
    minimize(0.5*SUM^2+C*sum(xi2));
    subject to
        Y.*(X*w2+b2)-1+xi2 >= 0;
        xi2 >= 0;
cvx_end

%% Dual Problem of SKM

e = ones(n,1);
cvx_begin
    variables gama alpha(n);
    dual variables eq ie;
    maximize(-1/2*gama^2+alpha'*e);
    subject to
        ie : alpha >= 0;
        alpha <= C;
        eq : alpha'*Y == 0;
        SUM = alpha(1)*Y(1)*X(1,:);
        for i = 2:n
            SUM = SUM+alpha(i)*Y(i)*X(i,:);
        end
        norm(SUM,2) <= gama;
cvx_end

%%
% As can be seen from the optimization results, the optimal values of
% primal problems are equal to that of the dual problem, which proves that
% the implementation is correct.
% 
% Calculate the variables in the primal problem.
w3 = X'*(Y.*alpha);
b3 = -eq;

xt = linspace(min(X(:,1)),max(X(:,1)),n);
yt = -(xt'*w3(1)+b3)/w3(2);

%% Kernelized Dual Problem of SKM
% Generate multiple kernels.
% 
% function d=sqdist(a,b)
% 
% aa = sum(a.*a,1); 
% 
% bb = sum(b.*b,1);
% 
% ab = a'*b;
% 
% d = abs(repmat(aa',[1 size(bb,2)]) + repmat(bb,[size(aa,2) 1]) - 2*ab);

num_kernels=2;
kernel_list=zeros(n*num_kernels,n);
arg_list = [1:num_kernels];
idx_list = zeros(1,1);
for i=1:num_kernels-1,
    idx_list(end+1,:) = floor(m/num_kernels)*i;
end
idx_list(end+1,:) = m;

for i=1:size(idx_list)-1,
    X_temp = X(:,idx_list(i,1)+1:idx_list(i+1,1));
    K = exp(-0.5/i*sqdist(X_temp',X_temp'));
    kernel_list(((i-1)*n+1):(i*n),:)=K;
end

Diganol_Y=diag(Y);

% Use cvx to solve the dual variables.
alpha2=zeros(n,1);
cvx_begin
    variables gamma2 alpha2(n,1);
    minimize(0.5*gamma2^2-sum(alpha2))
    subject to
        alpha2>=0;
        alpha2<=C;
        alpha2'*Y==0;
        for i=1:num_kernels,
            alpha2'*Diganol_Y*kernel_list(((i-1)*n+1):(i*n),:)*...
                Diganol_Y*alpha2<=gamma2;
        end 
cvx_end

% Calculate primal variables.
w_phi = 0;  % calculate w*phi
for i=1:n,
    for j=1:size(idx_list)-1,
        X_temp = X(:,idx_list(j,1)+1:idx_list(j+1,1));
        w_phi = w_phi+1/num_kernels*alpha2(i,1)*Y(i,1)*exp(-0.5/j*sqdist(X_temp',X_temp(i,:)'));
    end
end

idx_bias = find((alpha2>0.01)&(alpha2<0.99));  % calculate bias
b_k = Y(idx_bias,:)-w_phi(idx_bias,:);
bias = mean(b_k);
predict = w_phi + bias;
preidct_label = sign(predict);

% Calculate the percentage of correct prediction.
corr=sum(preidct_label==Y)/size(Y,1);

%% Results
% Plot the data points and the decision boundary for a linear SVM classifier.

figure(1)
plot(X(find(Y==-1),1), X(find(Y==-1),2),'.r','MarkerSize',10);
hold on
plot(X(find(Y==1),1), X(find(Y==1),2),'.b','MarkerSize',10);
axis equal
plot(xt,yt,'k','LineWidth',1);
hold off
legend('Class A','Class B','Decision Boundary');
title('Linear SVM Classifier');

%% 
% As can be seen from the figure above, when the data has complex distribution, 
% the linear classifier fails to classify two classes. We need to map the data 
% to a higher dimensional space and then apply the linear classifier in the high 
% dimensional space. However, this process can be simplified by a kernel function.
% The result is shown below.

figure(2)
plot(X(Y==1,1),X(Y==1,2),'.b','MarkerSize',10);
hold on
plot(X(Y==-1,1),X(Y==-1,2),'.r','MarkerSize',10);

d = 0.02;
[x1Grid,x2Grid] = meshgrid(min(X(:,1)):d:max(X(:,1)),...
    min(X(:,2)):d:max(X(:,2)));
xGrid = [x1Grid(:),x2Grid(:)];
d = 0.02;
[x1Grid,x2Grid] = meshgrid(min(X(:,1)):d:max(X(:,1)),...
    min(X(:,2)):d:max(X(:,2)));
xGrid = [x1Grid(:),x2Grid(:)];

kernel_train_test = zeros(size(xGrid,1), size(X,1));
for i=1:size(idx_list)-1,
    X_temp = X(:,idx_list(i,1)+1:idx_list(i+1,1));
    Grid_temp = xGrid(:,idx_list(i,1)+1:idx_list(i+1,1));
    K = 1/num_kernels*exp(-0.5/i*sqdist(Grid_temp',X_temp'));
    kernel_train_test = K + kernel_train_test;
end

prediction_test=kernel_train_test*(alpha2.*Y)+bias;
contour(x1Grid,x2Grid,reshape(prediction_test,size(x1Grid)),[0 0],'k');
legend('Class A','Class B','Decision Boundary');
title('SKM Classifier with Multiple Kernels');
hold off

%% 
% As shown in the figure, SKM yields a better classification boundary 
% for nonlinear cases.

%% Conclusions
% In the original paper[1], F. R. Bach et al. compared the computational complexity of SKM problem 
% based on sequential minimization optimization (SMO) and Mosek. In this
% project, we generated our own data instead of using real-life data in the original paper. 
% We did not take running time into consideration, but compared the performance of linear SVM with that of SKM with multiple kernels.  
% Using a conic combination of kernels, the classification problem is
% considered as QCQP, and thus can be solved by CVX. As can be seen from the results:
% 
% * The optimal values of the primal and dual problems are identical, which means we have solved both problems correctly.
% 
% * SKM presents a better performance for nonlinear cases.

%% Reference
% [1]Bach, Francis R., Gert RG Lanckriet, and Michael I. Jordan. "Multiple kernel learning, conic duality, and the SMO algorithm." Proceedings of the twenty-first international conference on Machine learning. ACM, 2004.
 