%% Distance Metric Learning for Large Margin Nearest Neighbor Classification
%% Introduction
% The k-Nearest Neighbour(kNN) rule classifies each unlabeled example by the majority label of its k-nearest 
% neighbors( based on distance measure) in the training data. Though the
% idea is simple, it yields competetive results. The performance of the kNN
% depends on the way the distances are computed between examples in the
% data. More often when no prior knowledge is available euclidean distance
% between examples. Euclidean distance may not be the right distance metric
% and can not capture all the properties in the data. In this project,
% Mahalanobis distance metric, which is a linear transformation of the input space, that precedes the kNN classification is learnt. The distance metric
% learning is modeled as a problem in convex optimization. 

%% Problem Formulation
% The algorithm is based on the observation that kNN will accurately
% classifiy if its k nearest neighbour belong to the same class label. The algorithm increases the number
% of samples from the training data with this property by learning a linear transformation of the input space that
% precedes kNN classification using Euclidean distances. The linear
% transfomration is derived by minimizing a loss function,  The first term penalizes large distances between
% samples in the same class desired as k-nearest neighbors and the second term penalizes
% small distances between samples with different labels. The euclidean
% distance in transformed space is equivalent to mahalanobis distance in
% the input space. 
%% Basic Notations
% First a few basic formulae and notations are introduced which is used later on in the
% convex optimization. Let $\vec{x}$ be the d-dimensional samples. Let the
% linear transformation be $\textbf{L}$, which is a dxd dimensional matrix.
% The tranformed or projected samples are then $\grave{x}$=
% $\textbf{L}\vec{x}$. The squared euclidean distance in the transformed domain is
% 
% $$ D_{L}(\grave{x_i},\grave{x_j})= ||\textbf{L}(\vec{x_i}-\vec{x_j})||_{2}^{2}$$
% 
% It is common to express squared distance metrics as square matrix as
% follow
% 
% $$ \textbf{M}=\textbf{L}^T\textbf{L} $$
% 
% Any matrix $\textbf{M}$ formed from real valued matrix \textbf{L} is
% guaranteed to be a positive semi-definite matrix. The squared distances interms of \textbf{M} is
% 
% $$ D_{M}(\grave{x_i},\grave{x_j})=
% (\vec{x_i}-\vec{x_j})^{T}M(\vec{x_i}-\vec{x_j}) $$
% 
% For each sample $\vec{x_i}$, a target nighbour $\vec{x_j}$ is defined as one
% of the k closest neighbour of the same class. This is denoted as $j \to i$,
%  $\vec{x_j}$ is a target neighbour to $\vec{x_i}$. For each sample there are k target neighbours(same class). Note that $j \to i$ does not
% imply the reverse is true. For kNN classification to succeed, target neighbors of each input $\vec{x_i}$ should be closer than all
% differently labeled inputs. For each input a perimeter can be imagined around the target neighbours. 
% A diffferently labeled input should not inavde this perimeter. These
% inputs are defined as impostors. The important goal of the learning is to
% redice the number of impostors. A large distance is to be maintained
% between the impostors and the perimeter defined by the target neighbours.
% By ensuring a margin of safety, the algorithm is robust to noise in the
% training data. Hence the name Large Margin Nearest Neighbour(LMNN)
% classification. For an input  $\vec{x_i}$ with label $y_i$
% and target neighbor $\vec{x_j}$
% , an impostor is any input $\vec{x_l}$ with a different label $y_l$ such that
% 
% $$ D_{L}(\grave{x_i},\grave{x_l})<=D_{L}(\grave{x_i},\grave{x_j})+1 $$
% 

%% Loss Function
% The loss function consists of two terms, one which pulls target
% neighbors closer together, and another which pushes differently labeled examples further apart. 
% The two terms have competing effects, since the first is reduced by shrinking the distances between
% examples while the second is generally reduced by magnifying them. The first term in the loss function penalizes large distances between 
% each input and its target
% neighbors. In terms of the linear transformation $\textbf{L}$ of the input space, the sum of these squared
% distances is given by
% 
% $$ \textbf{E}_{pull}=\Sigma_{j\to i}||\textbf{L}(\vec{x_i}-\vec{x_j})||_{2}^{2} $$
% 
% It does not penalize all the samples with same labels, but only the
% target neighbours. The second term in the loss function penalizes small distances between differently labeled examples.
% To simplify notation,
% we introduce a new indicator variable $y_{il}$ =1 if and only if $y_i =y_l$
% , and $y_{il}$ =0. The second term is
%  
% $$ \textbf{E}_{push}=\Sigma_{i, j\to
% i}\Sigma_l(1-y_{il})[1+||\textbf{L}(\vec{x_i}-\vec{x_j})||_{2}^{2}-||\textbf{L}(\vec{x_i}-\vec{x_l})||_{2}^{2}]_+$$
% 
% where the term $[z]_+= max(z,0)$, denotes the standard hinge loss. The
% final loss function is
% 
% $$ \textbf{E(L)}= (1-u)\textbf{E}_{pull} +u\textbf{E}_{push}$$
% 
% The above can be rewritten in terms of matrix $\textbf{M} as follows,
% 
% $$ \textbf{E(M)}=(1-u)\Sigma_{j\to i}D_{M}(\grave{x_i},\grave{x_j})+ u\Sigma_{i, j\to
% i}\Sigma_l(1-y_{il})[1+ D_{M}(\grave{x_i},\grave{x_j})-D_{M}(\grave{x_i},\grave{x_l})]_+$$
% 
% With this substitution, the loss function is now expressed over positive semidefinite matrices $\textbf{M}\succeq 0$,
% as opposed to real-valued matrices $\textbf{L}$. The loss function is a piecewise linear, convex function of the elements in the matrix
% $\textbf{M}$. The first term(penalizing large distances between target neighbors)
% is linear in the elements of M, while the second term (penalizing impostors) is derived from
% the convex hinge loss. The loss function is brought to a standard form,
% formulated as an semidefinite programming(SDP) by introducing slack variables $e_{ijl}$ for all
% triplets of target neighbours $j\to i$ and impostots $\vec{x_l}$. The
% slack variable $e_{ijl}\geq 0$, is used to measure the amount by which the large margin
% inequality is violated. The SDP is as follows
% 
% $$ \textbf{Minimize}_{M, e_{ijl}} (1-u)\Sigma_{j\to i}(\vec{x_i}-\vec{x_j})^{T}M(\vec{x_i}-\vec{x_j})+ u\Sigma_{i, j\to
% i}\Sigma_l(1-y_{il})e_{ijl}$$
% 
% $$\textbf{Subject to} $$
% 
% $$  (\vec{x_i}-\vec{x_l})^{T}M(\vec{x_i}-\vec{x_l})-(\vec{x_i}-\vec{x_j})^{T}M(\vec{x_i}-\vec{x_j}) \geq
% 1-e_{ijl} $$
% 
% $$ M \succeq 0$$
% 
% $$ e_{ijl} \geq 0$$
%% Code Run through
% The data used here is 4 dimensional data with 429 training samples. The slack variable used above is a combination of 3
% incides leading to a total of 22787 inequalities. For computational sake
% the data is sampled. The nearest neighbour is set to 3.
load bal.mat % xtr yTr xTe yTe
sample=80;
xTr=xTr(:,35:end);
yTr=yTr(35:end);
y=randsample(size(xTr,2),sample); %radom sampling
y=sort(y);
K=3; % nearest neighbour
xTr=xTr(:,y);
yTr=yTr(:,y);
% test_sample=randsample(length(yTe),50);
% test_sample=sort(test_sample);
%% 
% The target neighbours are calculated for each sample and their
% corresponding distances.
[d,n]=size(xTr);
un=unique(yTr);
NN=zeros(K,n);
for c=un
    i=find(yTr==c);
    nn=LSKnn(xTr(:,i),xTr(:,i),2:K+1); % gives the indices of k nearest neighbours
    NN(:,i)=i(nn);
end
Ni=zeros(K,n);
dis=zeros(K,n);
%target neighbours
 for nnid=1:K
        Ni(nnid,:)=sum((xTr-xTr(:,NN(nnid,:))).^2,1)+1; % euclidean distance
        disvector{nnid}=xTr-xTr(:,NN(nnid,:));
 end;
 
    
%% 
% The impostors are calcuated for each sample $\vec{x_i}$ and target
% neighbour $\vec{x_j}$, the samples that violate the perimeter defined by
% the target neighbours.
for c=un
  
   i=find(yTr==c);
   index=find(yTr~=c);
   X1=xTr(:,i);
   X2=xTr(:,index);
   distmat=zeros(4,length(i),length(index));
   for a=1:length(i)
       for b=1:length(index)
       
       distmat(:,a,b)=X1(:,a)-X2(:,b);
       end
   end
   distancevector{c}=distmat;
end
for i=1:sample
distmat=[];
index=find(yTr~=yTr(i));
dist=distance(xTr(:,i),xTr(:,index));
id=find(dist<=Ni(3,i)+1); %violating the perimeter
% for j=1:length(index)
%     distmat(:,j)=xTr(:,i)-xTr(:,j);
% end
distancevector2{i}=index(id)';
end
%%
% The total number of impostors are calculated for each sample
total=0;
for i=1:sample
    len(i)=size(distancevector2{i},1);

end
total=sum(len);
num_imp=zeros(1,length(len)+1);
num_imp(2:end)=len(1:end);
cc=cumsum(num_imp);
%% CVX Package
% The cvx is run for variables $\textbf{M(d,d)}$ and $e_{ijl}$
cvx_begin
 variable M(d,d) symmetric
 variable e(total*3) nonnegative
 f_x=0;
 for i=1:sample
     for j=1:K
         f_x=f_x+(xTr(:,i)-xTr(:,NN(j,i)))'*M*(xTr(:,i)-xTr(:,NN(j,i)));
     end
 end
 for i=1:sample
     
     for l=1:len(i)
         for j=1:K
             f_x=f_x+e((cc(i))+(l-1)*3+j);
         end
     end
 end
 minimize f_x
 subject to
 for i=1:sample
     
      for l=1:len(i)
         for j=1:K
             (xTr(:,i)-xTr(:,distancevector2{i}(l)))'*M*(xTr(:,i)-xTr(:,distancevector2{i}(l)))-(xTr(:,i)-xTr(:,NN(j,i)))'*M*(xTr(:,i)-xTr(:,NN(j,i)))>=1-e((cc(i))+(l-1)*3+j);
         end
     end
 end
 M>=0
 e>=0
 cvx_end
%% Results
% The projection matrix $\textbf{L}$ is calculated for the proposed
% technique and the results are compared with four algorithms. 1)kNN before LMNN, 2)
% kNN after LMNN, 3) kNN after PCA. We see that
% kNN performs well after LMNN compared to the other two methods
L_lmnn=sqrt(M);
L_knn=[];
[L_pca,score,latent,tsquared,explained]= pca(xTr','Algorithm','eig','Centered',false);
L_pca=L_pca';
Pca_err=knncl(L_pca,xTr,yTr,xTe,yTe,1,'train',0);
knn_err=knncl(L_knn,xTr,yTr,xTe,yTe,1,'train',0);
lmnn_err=knncl(L_lmnn,xTr,yTr,xTe,yTe,1,'train',0);
fprintf('Bal data set\n');
fprintf('\n\nTesting error LMNN: %2.2f%%\n',100*lmnn_err);
fprintf('Testing error kNN: %2.2f%%\n',100*knn_err);
fprintf('Testing error PCA: %2.2f%%\n',100.*Pca_err);
%% Conclusion and Discussion
% A new framework for large margin nearest neighbor (LMNN)
% classification is introduced. From labeled training examples, a Mahalanobis distance
% metric is learned for kNN classification. The required optimization is formulated as a semidefinite programming. The framework makes no parametric assumptions about the structure or
% distribution of the data, though it does not scale well for large data sets. The accuracy of kNN classification
% is imporved significantly by learning a metric in this way.

%% References
% 1)Kilian Q. Weinberger and Lawrence K. Saul. 2009. Distance Metric Learning for Large Margin Nearest Neighbor Classification. J. Mach. Learn. Res. 10 (June 2009), 207-244.


 