# Classification Problem  (Assignment 5, TAO, Spring 2019)
### Instructor: Dr. Pawan Kumar

## Classification Problem
### Given a set of input vectors corresponding to objects (or featues) decide which of the N classes the object belogs to.  

### Reference (some figures for illustration below are taken from this): 
1. SVM without Tears, https://med.nyu.edu/chibi/sites/default/files/chibi/Final.pdf

## Possible Methods for Classification Problem
1. Perceptron
2. SVM
3. Neural Networks

## SVM: Support Vector Machines

We will briefly describe the idea behind support vector machines for classification problems. We first describe linear SVM used to classify linearly separable data, and then we describe how we can use these algorithm for non-linearly separable data by so called kernels. The kernels are functions that map non-linearly separable data to a space usually higher dimensional space where the data becomes linearly separable. Let us quickly start with linear SVM.

### Linear SVM for two class classification
We recall the separating hyperpplane theorem: If there are two non-intersecting convex set, then there exists a hyperplane that separates the two convex sets. This is the assumption we will make: we assume that the convex hull of the given data leads to two convex sets for two classes, such that a hyperplane exists that separates the convex hulls. 

### Main idea of SVM: 
Not just find a hyperplane (as in perceptrons), but find one that keeps a good (largest possible) gap from the the data samples of each class. This gap is popularly called margins.

### Illustration of problem, and kewords
Consider the dataset of cancer and normal patients, hence it is a two class problem. Let us visualize the data:

Let us notice the following about the given data:
0. There are two classes: blue shaded stars and red shaded circles.
2. The input vector is two dimensional, hence it is of the form $(x_1^1, x_2^1).$
2. Here $x_1^1, x_2^2$ are values of the features corresponding to two gene features: Gene X, Gene Y.
3. Here red line is the linear classifier or hyperplane that separates the given input data.
4. There are two dotted lines: one passes through a blue star point, and another dotted line passes through two red shaded circle points.
5. The distance between the two dotted lines is called gap or margin that we mentioned before.
6. Goal of SVM compared to perceptrons is to maximize this margin.

## Formulation of Optimization Model for Linear SVM
We now assume that the red hyperplane above with maximum margin is given by $$w \cdot x + b,$$
We further assume that the dotted lines above are given by $$w \cdot x + b = -1, \quad w \cdot x + b = +1.$$

For reasons, on why we can choose such hyperplane is shown in slides Lecture 16 of TAO. Since we want to maximize the margin the distance between the dotted lines, we recall the formula for diatance between planes. Let $D$ denote the distance, then 
$$D = 2/ \| w \|.$$
So, to maximize the margin $D,$ we need to minimize $\| w \|.$ For convenience of writing algorithm (for differentiable function), we can say that minimizing $\| w \|$ is equivalent to minimizing $1/2 \| w \|^2.$ Hence 
### Objective function: $\dfrac{1}{2} \| w \|^2$
For our hyperplane to classify correctly, we need points of one class on one side of dotted line, more concretely
$$w \cdot x + b \leq -1,$$
and the we want the samples of another class (red ones) be on the other side of other dotted lines, i.e., 
$$ w \cdot x + b \geq +1.$$
Let us now look what constraints mean in figure:

With this we are all set to write the constraints for our optimization model for SVM. 
### Constraints: 
$$
\begin{align}
&w \cdot x_i + b \leq -1, \quad \text{if}~y_i = -1\\
&w \cdot x_i + b \geq +1, \quad \text{if}~y_i = +1
\end{align}
$$
Hence, objective function with constraints, gives us the full model. The data for which the label $y_i$ is $-1$ satisfies $w \cdot x + b \leq -1,$ and the data for which the lable $y_i$ is $+1$ satisfies $w \cdot x + b \geq +1.$ Hence both these conditions can be combined to get
$$
\begin{align}
y_i (w \cdot x_i + b) \geq 1
\end{align}
$$

## Optimization Model (Primal Form):
$$
\begin{align}
\text{minimize} \quad & \dfrac{1}{2} \| w \|^2 \\
\text{subject to} \quad &w \cdot x_i + b \geq 1, \quad i=1,\dots,m,
\end{align}
$$
where $m$ is the number of samples $x_i,$ and $w \in \mathbb{R}^n.$

$\color{red}{Question:}$ Prove that the primal objective is convex.


$\color{blue}{Answer}:$ 

$$
\begin{align}
\text{First, we need to prove that} \quad \dfrac{1}{2} \|w \|^2 \quad \text{is a convex function} \\
\text{Let x and y be two vectors} \quad \quad \quad \quad  \\ \\
Now, \quad \theta ||x||^2 + (1- \theta)||y||^2 - || \theta .x + (1- \theta).y ||^2 \\ \\
=> \theta ||x||^2 + (1- \theta)||y||^2 - \theta^2.||x||^2 -\theta^2||y||^2 -2\theta(1-\theta)<x|y> \\ \\
\text{using the property that } <x|y> \leq |x|.|y| \quad \quad \quad \\ \\
\text{we can write,} \quad \quad \quad \quad \quad \quad\\
\theta ||x||^2 + (1- \theta)||y||^2 - \theta^2.||x||^2 -\theta^2||y||^2 -2\theta(1-\theta)<x|y>\quad \\ \\ \geq  \theta(1-\theta)||x||^2 + \theta(1-\theta)||y||^2 -2\theta(1-\theta)||x|||y|| \\ \\
= \theta(1-\theta)(||x||-||y||)^2 \geq=0 \qquad \qquad \qquad \\ \\ \\
=> \theta ||x||^2 + (1- \theta)||y||^2 \geq  || \theta .x + (1- \theta).y ||^2 \\
\text{Hence proved that the objective function which is} \quad \dfrac{1}{2}||w||^2 \quad \text{is a convex function}
\end{align}
$$


$\color{red}{Question}:$ Write the primal problem in standard form.

$\color{blue}{Answer}:$ 
$$
\begin{align}
\text{The standard form of a primal problem is given by:}\\
\text{minimise} \quad f(x) \qquad \qquad \qquad \\
\text{subject to} \quad g_i(x) \leq 0 \quad i=1.......m\\
h_j(x) = 0 \quad j=1.........p \\
\text{Now,}\quad y_i(wx_i+b) \geq 1, \quad i=1.....m \quad \text{can be written as}\\
1-y_i(wx_i+b) \leq 0 \quad for \quad  i=1.......m \\ \\
\text{Thus, the optimisation problem can be formulated as:} \\ \\
\text{minimise} \quad \dfrac{1}{2}||w||^2 \qquad \qquad \qquad \qquad \qquad \\
\text{subject to} \quad 1-y_i(wx_i+b) \leq 0 \quad i=1.......m\\
\end{align}
$$


## Optimization Model (Dual Form)
The dual form was derived in lecture 16:
$$
\begin{align*}
&\text{maximize} \quad \sum_{i=1}^m{\lambda_i} - \dfrac{1}{2} \sum_{i=1}^m \lambda_i \lambda_j y_i y_j (x_i \cdot x_j) \\
&\text{subject to} \quad \lambda_i \geq 0, \quad \sum_{i=1}^m{\lambda_i y_i} = 0, \quad i = 1, \dots, m
\end{align*}, 
$$
where $\lambda_i$ is the Lagrange multiplier. We claim that strong duality holds. 

$\color{red}{Question:}$ Show the derivation of dual.

$\color{blue}{Answer:}$ 

$$
\begin{align*}
\text{The Lagrangian can be written as:} \quad \quad \quad \quad \quad\\
\text{minimise} L(w,b,\lambda) = \dfrac{1}{2} ||w||^2 + \sum_{i=1}^m\lambda_i(1-y_i(w^Tx_i+b)) \\
\text{subject to} \quad \lambda_i \geq 0 \quad \text{for all i=1,2.....m} \quad \quad \quad \quad \quad \\
\text{Now, to minimise we will have to differentiate with respect to w as well as with respect to b}\\
=> \dfrac{dL}{dw}=0= w- \sum_i^m(\lambda_i*y_i*x_i) \quad \quad \quad \quad\\
\text{Also,} \dfrac{dL}{db}=0= - \sum_{i=1}^m \lambda_i y_i=0 \quad \quad \quad \quad\\
\text{Substituting the value of } w=\sum_{i=1}^m(\lambda_i*y_i*x_i) and \sum_{i=1}^m(\lambda_i*y_i)=0, we get \\
\text{L}(\lambda,b)=\dfrac{1}{2}w^T\sum_{i=1}^m\lambda_i*y_i*x_i-\sum_{i=1}^m\lambda_i[y_i(w^TX_i+b)-1]\\
\text{L}(\lambda,b)=\dfrac{1}{2}w^T\sum_{i=1}^m\lambda_i*y_i*x_i-w^T\sum_{i=1}^m\lambda_i*y_i*x_i-b\sum_{i=1}^m\lambda_i*y_i+\sum_{i=1}^m \lambda_i \\
\text{L}(\lambda,b)=\sum_{i=1}^m\lambda_i-b\sum_{i=1}^m\lambda_i*y_i-\dfrac{1}{2}\sum_{i=1}^m\sum_{j=1}^m\lambda_i*\lambda_j*y_i*y_j*(X_i^T*X_j) \\
But, \sum_{i=1}^m\lambda_i*y_i=0 \quad \quad \quad \quad \quad\\
\text{So, the final equation will be:} \quad \quad \quad \quad \\
\text{maximise}\quad \text{L}(\lambda)=\sum_{i=1}^m\lambda_i-\dfrac{1}{2}\sum_{i=1}^m\sum_{j=1}^m\lambda_i*\lambda_j*y_i*y_j*(X_i^T*X_j) \\
\text{subject to}\quad \lambda_i \geq 0 \quad \text{for all i=1,2......m} \quad \text{and} \sum_{i=1}^m \lambda_i*y_i=0
\end{align*}
$$



$\color{red}{Question:}$ Prove that strong duality holds.

$\color{blue}{Answer:}$ 

$$
\begin{align*}
\text{our original problem was to minimise} \quad \dfrac{1}{2}||x||^2 \quad \quad \quad\\
\text{subject to} \quad y_i(w^T*x + b) \geq 1 \quad \text{for i=1........m} \quad \quad \quad\\
\text{Now, to prove the strong duality, we would use slater conditions according to which,given}\quad \quad \quad \\
\text{minimise} f_0(x) \quad \quad \quad \quad \quad \quad\quad \quad \quad \\
\text{subject to} \quad f_i(x) \leq 0 \quad  i=1.....m  \quad \quad \quad\quad \quad \quad\\
Ax=b \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad\\
\text{where all the functions are convex and strict inequality holds for a feasible point such that} \quad \quad \quad\\
f_i(x)  <  0 \quad \text{for i =1.......m} \quad \quad \quad \quad \quad \quad\\ \\
\text{Then, strong duality holds} \quad \quad \quad \quad\quad \quad \quad\\ \\
\text{in case of svm, the support vectors are the ones which follow equality whereas all the other points follow strict inequality} \quad \quad \quad\\ \\
\text{this means that there exists a point such that} \quad \quad \quad y_i(wx_i+b) > 1 \quad \quad \quad \quad \quad \quad\\ \\ 
\text{also, all the objective and constraint functions(affine) are convex} \quad \quad \quad \quad \quad \quad\\ \\ 
\text{thus, slater conditions and strong duality hold} \quad \quad \quad \quad \quad \quad \quad \quad \quad\\
\end{align*}
$$

$\color{red}{Question:}$ Prove that the dual objective is concave.

$\color{blue}{Answer}:$ 
$$
\begin{align}
\text{the dual objective is given by:}\\
\sum_{i=1}^m{\lambda_i} - \dfrac{1}{2} \sum_{i=1}^m \lambda_i \lambda_j y_i y_j (x_i \cdot x_j)\\
\text{this can be written in the form of a quadratic equation like this:}\\
\lambda^T.1 - \dfrac{1}{2}\lambda^T A \lambda \\
=>-(\dfrac{1}{2}\lambda^T A \lambda - \lambda^T.1) \\ \\
\text{Here,}A_{ij}=y_i.y_j.x_i^T.x_j \quad \quad\\
\text{This is basically a negative of convex quadratic function which is given by} \dfrac{1}{2}X^TPX+q^TX \\
\text{We Know that the negative of a convex function is a concave function}\\
\text{Hence Proved} \\
\end{align}
$$

$\color{red}{Question}:$ Write the dual problem in standard form.

$\color{blue}{Answer}:$ 
$$
\begin{align}
\text{The standard form is given by:}\\
\text{minimise} \quad f(x) \\
\text{subject to} \quad g_i(x) \leq 0 \quad i=1.......m\\ \\ \\
\text{Now, the dual problem is,}\\
\text{maximize} \quad \sum_{i=1}^m{\lambda_i} - \dfrac{1}{2} \sum_{i=1}^m \lambda_i \lambda_j y_i y_j (x_i \cdot x_j) \\
\text{subject to} \quad \lambda_i \geq 0, \quad \sum_{i=1}^m{\lambda_i y_i} = 0, \quad i = 1, \dots, m\\ \\ \\
\text{This can be written as} \\
\text{minimise} \quad \dfrac{1}{2} \sum_{i=1}^m \lambda_i \lambda_j y_i y_j (x_i \cdot x_j) -\sum_{i=1}^m{\lambda_i}\\
\text{subject to} \quad \quad -\lambda_i \leq 0 \quad \sum_{i=1}^m{\lambda_i y_i} = 0, \quad i = 1, \dots, m\\
\end{align}
$$

## Soft Margin SVM
In a variant of soft margin SVM, we assume that some data samples may be outliers or noise, and this prevents the data from being linearly separable. For example, see the figure below

In the figure, we see that 

- We believe that two red and one blue sample is noisy or outliers.
- We now want to take into account that real life data is noisy, we decide to allow for some of the noisy data in the margin.
- Let $\xi_i$ denotes how far a data sample is from the middle plane (margin is the area between dotted line).
- For example, one of the noisy red data point in 0.6 away from middle red plane. 
- We introduce this slack variable $\xi_i \geq 0$ for each data sample $x_i.$

## Optimization Model: Primal Soft-Margin
We can then write the primal soft-margin optimization model as follows:
$$
\begin{align*}
&\text{minimize} \quad \dfrac{1}{2} \| w \|^2 + C \sum_{i=1}^m \xi_i \\
&\text{subject to} \quad y_i (w \cdot x_i + b) \geq 1 - \xi_i, \quad \xi_i \geq 0, \quad i = 1, \dots, m.
\end{align*}
$$

## Optimization Model: Dual Soft-Margin
We can also write the dual form of soft-margin SVM as follows:
$$
\begin{align*}
\text{Maximize} \quad &\sum_{i=1}^m \lambda_i - \dfrac{1}{2} \sum_{i,j=1}^m \lambda_i \lambda_j \: y_i y_j \: x_i \cdot x_j \\
\text{subject to} \quad &0 \leq \lambda_i \leq C, \quad  i=1, \dots, m, \\ 
&\sum_{i=1}^m \lambda_i y_i = 0.	 
\end{align*}
$$

$\color{red}{Question:}$ Show the derivation of dual.

$\color{blue}{Answer:}$ 
$$
\begin{align}
\text{We would have to first convert it into lagrangian} \\
=>L(w,b,\lambda,\mu)= \dfrac{1}{2}||w||^2+C\sum_{i=1}^m\xi_i -\sum_{i=1}^m\lambda_i[(y_i(wx_i+b)-1+\xi_i]-\sum_{i=1}^m\mu_i\xi_i\\
=> \text{Now, to find the minimimum we would equate the gradient as zero with respect to w,b,}\xi \\
=> \dfrac{dL}{dw}=w-\sum_{i=1}^m \lambda_i*y_i*x_i=0 \\
=> w=\sum_{i=1}^m \lambda_i*y_i*x_i \\
\text{also,} \quad \quad \quad \quad \quad \qquad\\
=>\dfrac{dL}{db}=\sum_{i=1}^m \lambda_i*y_i=0 \\
\text{also,} \quad \quad \quad \quad \quad \qquad\\
=>\dfrac{dL}{d\xi_i}=C-\lambda_i-\mu_i=0 \quad \quad\\
\text{as the Lagrangian can be expected in the form} \quad L= ......\sum_{i=1}^mC\xi_i-\sum_{i=1}^m\lambda_i\xi_i-\sum_{i=1}^m\mu_i\xi_i \\
=> C= \lambda_i + \mu_i \quad \quad\\
\text{Thus,} \lambda_i \leq C \quad as \quad \mu_i\geq 0\\ \\ \\
\text{Substituting the value of } w=\sum_{i=1}^m(\lambda_i*y_i*x_i) and \sum_{i=1}^m(\lambda_i*y_i)=0, we get \\
\text{L}(\lambda,b)=\dfrac{1}{2}w^T\sum_{i=1}^m\lambda_i*y_i*x_i-\sum_{i=1}^m\lambda_i[y_i(w^TX_i+b)-1]+ \sum_{i=1}^mC\xi_i-\sum_{i=1}^m\lambda_i\xi_i-\sum_{i=1}^m\mu_i\xi_i\\
\text{But we know that}\sum_{i=1}^mC\xi_i=\sum_{i=1}^m\lambda_i\xi_i+\sum_{i=1}^m\mu_i\xi_i\\
\text{L}(\lambda,b)=\dfrac{1}{2}w^T\sum_{i=1}^m\lambda_i*y_i*x_i-w^T\sum_{i=1}^m\lambda_i*y_i*x_i-b\sum_{i=1}^m\lambda_i*y_i+\sum_{i=1}^m \lambda_i \\ \\ \\
\text{L}(\lambda,b)=\sum_{i=1}^m\lambda_i-b\sum_{i=1}^m\lambda_i*y_i-\dfrac{1}{2}\sum_{i=1}^m\sum_{j=1}^m\lambda_i*\lambda_j*y_i*y_j*(X_i^T*X_j) \\
But, \sum_{i=1}^m\lambda_i*y_i=0 \quad \quad \quad \quad \quad\\
\text{So, the final equation will be:} \quad \quad \quad \quad \\
\text{maximise}\quad \text{L}(\lambda)=\sum_{i=1}^m\lambda_i-\dfrac{1}{2}\sum_{i=1}^m\sum_{j=1}^m\lambda_i*\lambda_j*y_i*y_j*(X_i^T*X_j) \\
\text{subject to}\quad \lambda_i \geq 0 ,\quad \lambda_i \leq C\quad \text{for all i=1,2......m} \quad \text{and} \sum_{i=1}^m \lambda_i*y_i=0
\end{align}
$$


$\color{red}{Question:}$ List advantages of dual over primal.

$\color{blue}{Answer:}$ 
$$
\begin{align}
1.\quad \quad \quad\text{ The most significant benefit from solving the dual comes when you are using} \\
\text{the "Kernel Trick" to classify data that is not linearly separable in the original feature space.} \\ \\
2.\text{If D(number of features) is too large, then computing this mapping (and solving for W using the} \\
\text{primal formulation) becomes too expensive. The "trick" part of the the kernel trick comes}\\
\text{here:  we can avoid computing this mapping altogether. Suppose we require only the dot-product between }\\
\text{two data points in the mapped feature space. Now suppose you have a function (unsurprisingly named }\\
\text{the "Kernel Function") which directly computes this dot-product between points in R^D from}\\
\text{their corresponding points in R^d.This is almost the same as the original dual formulation, except}\\
\text{that you compute the kernel function instead of the ordinary dot-product. Note that this is not}\\
\text{possible in the primal formulation, where it would be necessary to explicitly compute}\\
\text{the mapping for each data point.}\quad \quad \quad \quad\\ \\
3.\quad \quad \quad \text{for Data with a large number of features, say, P, the time complexity of the computation}\\
\text{would be of the order of NP where N is the number of samples and P is the number of features}\\
\text{But if we solve it using the dual formulation, then it would be of the order of}N^2.\\
\text{As only the dot products of samples with each other are needed, it would only be of the order of N^2}\\
\text{This is efficient for those problems where P>>>N, i.e. the number of features are way more than the number of samples.}
\end{align}
$$

# Kernels in SVM
## Non-Linear Classifiers
- For nonlinear data, we may map the data to a higher dimensional feature space where it is separable. See the figure below:

Such non-linear transformation can be implemented more effectively using the dual formulation. 
- If we solve the dual form of linear SVM, then the predictions is given by
$$
\begin{align*}
f(x) &= \text{sign}(w \cdot x + b) \\
w &= \sum_{i=1}^m \alpha_i y_i x_i 
\end{align*}
$$
If we assume that we did some transform $\Phi,$ then the classifier is given by
$$
\begin{align*}
f(x) &= \text{sign} (w \cdot \Phi(x) + b) \\
w &= \sum_{i=1}^m \alpha_i y_i \Phi(x_i)
\end{align*}
$$
If we substitute $w$ in $f(x),$ we observe that
$$
\begin{align*}
f(x) = \text{sign} \left ( \sum_{i=1}^m \alpha_i y_i \, \Phi(x_i) \cdot \Phi(x) + b \right) = \text{sign} \left( \sum_{i=1}^m \alpha_i y_i \, K(x_i, x) + b \right)
\end{align*}
$$
Note that doing dot products such as $\Phi(x_i) \cdot \Phi(x),$ if $\Phi(x)$ is a long vector! An important observation is to define this dot product or $K(x,z)$ such that dot products happen in input space rather than the feature space. We can see this with following example:
$$
\begin{align*}
K(x \cdot z) &= (x \cdot z)^2 = \left( \begin{bmatrix}
x_{(1)} \\ x_{(2)} 
\end{bmatrix} \cdot \begin{bmatrix}
z_{(1)} \\ z_{(2)}
\end{bmatrix} \right)^2 = (x_{(1)} z_{(1)} + x_{(2)} z_{(2)})^2 \\
&= x_{(1)}^2 z_{(1)}^2 + 2x_{(1)} z_{(1)} x_{(2)} z_{(2)} + x_{(2)}^2 z_{(2)}^2 = \begin{bmatrix}
x_{(1)}^2 \\ \sqrt{2} x_{(1)} x_{(2)} \\ x_{(2)}^2 
\end{bmatrix} \cdot \begin{bmatrix}
z_{(1)}^2 \\ \sqrt{2} z_{(1)} z_{(2)} \\ z_{(2)}^2 
\end{bmatrix}  \\
&= \Phi(x) \cdot \Phi(z)
\end{align*}
$$

$\color{red}{Question:}$ Let the kernel be defined by $K(x,z) = (x \cdot z)^3.$ Define $\Phi(x).$ Assuming that one multiplications is 1 FLOP, and one addition is 1 FLOP, then how many flops you need to compute $K(x \cdot z)$ in input space versus feature space?

$\color{blue}{Answer:}$ 
$$
\begin{align}
K(x.z)=(x.z)^3 \text{can be expaned as} \\
=>(x_1*z_1+x_2*z_2)^3= x_1^3*z_1^3+x_2^3*z_2^3+3*x_1^2*z_1^2*x_2*z_2+3*x_1*z_1*x_2^2*z_2^2 \\ \\
\begin{bmatrix}
\ x_1^3 \\
\ x_2^3 \\
\sqrt[2]{3}*x_1^2*x_2 \\
\sqrt[2]{3}*x_1*x_2^2 \\
\end{bmatrix}
\begin{bmatrix}
\ z_1^3 \\
\ z_2^3 \\
\sqrt[2]{3}*z_1^2*z_2 \\
\sqrt[2]{3}*z_1*z_2^2 \\
\end{bmatrix}
\text{So, in the feature space there are 4 multiplications and 3 additions involved which is equal to 7 flops}\\
\text{the number of multiplications are 4 as it is a 4 dimensional vector now and to add those 4, 3 additions are required}\\
\text{whereas, in the input space the number of multiplications are 4 and one addition is involved, so, overall 5}\\
\text{in the input space, the number of multiplications are 4 to first calculate the dot dproduct which requires 2}\\
\text{and to raise the power by three, two more multiplications are required}\\ \\ \\
\text{So, answer is 7 Flops in feature space and 5 in input space}\\
\end{align}
$$



## Optimization Model: Dual Soft Margin Kernel SVM
We can now write the dual form of soft-margin Kernel SVM as follows:
$$
\begin{align*}
\text{Maximize} \quad &\sum_{i=1}^m \lambda_i - \dfrac{1}{2} \sum_{i, \, j=1}^m \lambda_i \lambda_j \: y_i y_j \: \Phi(x_i) \cdot \Phi(x_j) \\
\text{subject to} \quad &0 \leq \lambda_i \leq C, \quad  i=1, \dots, m, \\ 
&\sum_{i=1}^m \lambda_i y_i = 0.	 
\end{align*}  
$$

## Solver for Optimization Problem: Quadratic Programming
We aspire to solve the above optimization problem using existing quadraric programming library. But we have a problem: the standard libraries use the standard form of the quadratic optimization problem that looks like the following:
$$
\begin{align*}
\text{minimize} \quad &\dfrac{1}{2} x^T P x + q^T x, \\ 
\text{subject to} \quad &Gx \leq h, \\
&Ax = b
\end{align*}
$$

# Dual Soft-Margin Kernel SVM in Standard QP: Assemble Matrices Vectors
To put the dual Kernel SVM in standard form, we need to set
- matrix $P$
- vector $x$
- vector $q$
- vector $h$
- vector $b$
- matrix $G$
- matrix $A$

### Matrix $P$
Let $$K(x_i, x_j) = \Phi(x_i) \cdot \Phi(x_j),$$ and set $(i,j)$ entry of matrix $P$ as $$P_{ij} = y_iy_j K(x_i,x_j)$$

## Vector $x$
Set $$x = \begin{bmatrix}
\lambda_1 \\
\lambda_2 \\
\vdots \\
\lambda_m
\end{bmatrix}
$$

## Vector $q$
Set $q \in \mathbb{R}^m$
$$ q = 
\begin{bmatrix}
-1 \\ -1 \\ \vdots \\ -1
\end{bmatrix}
$$

## Matrix $A$
Set the matrix (in fact vector) $A$ as 
$$
A = [y_1, y_2, \dots, y_m]
$$

## Vector $b$
In fact vector $b$ is a scalar here: $$b = 0$$

## Matrix $G$
$$
\begin{align*}
G = \begin{bmatrix}
1 & 0 & \dots & 0 \\
0 & 1 & \dots & 0 \\
\vdots & \ddots & \dots & \vdots \\
0 & 0 & \dots & 1 \\ \hline
-1 & 0 & \dots & 0 \\
\vdots & \ddots & \dots & \vdots \\
0 & 0 & \dots& -1
\end{bmatrix}
\end{align*}
$$

## Vector $h$
Set $h$ as 
$$
h = \begin{bmatrix}
C \\
C \\
\vdots \\
C \\
0 \\
\vdots \\
0
\end{bmatrix}
$$

# Implementation of Kernel SVM
We are all set to try out coding the classifier using Kernel SVM. We will first import some libraries. Some of these libraries may not be available in your system. You may install them as follows:
- conda install numpy
- conda install -c conda-forge cvxopt
- sudo apt-get install python-scipy python-matplotlib

Try google search, if these does not work.


```python
import pylab as pl
import cvxopt as cvxopt
from cvxopt import solvers
import numpy as np
```

We will now define a class: svm 
This class will have the following functions:
- __init__: where we will define initial default parameters
- *construct_kernel*: here we will define some kernels such as polynomial and RBF (radial basis or Gaussian kernel)
- *train_kernel_svm*: Here we will train, i.e, we will call a quadratic programming solver from cvxopt
- *classify*: Here we will test our classifier

$\color{red}{Question:}$ Fill the TODO below. 


```python
class svm:

    def __init__(self, kernel='linear', C=None, sigma=1., degree=1., threshold=1e-5):
        self.kernel = kernel
        if self.kernel == 'linear':
            self.degree = 1.
            self.kernel = 'poly'
            
        self.C = C
        self.sigma = sigma
        self.threshold = threshold
        self.degree = degree
        

    def construct_kernel(self, X):
        self.K = np.dot(X, X.T)
        print(np.diag(self.K)*np.ones((1, self.N)))

        if self.kernel == 'poly':
            self.K = (1. + 1./self.sigma*self.K)**self.degree

        elif self.kernel == 'rbf':
            self.xsquared = (np.diag(self.K)*np.ones((1, self.N))).T
            b = np.ones((self.N, 1))
            self.K -= 0.5*(np.dot(self.xsquared, b.T) +
                           np.dot(b, self.xsquared.T))
            self.K = np.exp(self.K/(2.*self.sigma**2))

    def train_kernel_svm(self, X, targets):
        self.N = np.shape(X)[0]
        self.construct_kernel(X)
        print(self.N)

        # Assemble the matrices for the constraints 
        P=np.zeros((self.N,self.N))
        for i in range(self.N):
            for j in range(self.N):
                P[i][j]=targets[i]*targets[j]*self.K[i][j]
                
        q = -np.ones((self.N,1))
        G=np.zeros((2*self.N,self.N))
        for i in range(self.N):
            for j in range(self.N):
                G[i][j]=1
        for i in range(self.N):
            for j in range(self.N):
                G[self.N+i][j]=-1
        h=np.zeros((2*self.N,1))
        for i in range(self.N):
            if self.C==None:
                h[i]=0
            else:
                h[i]=self.C
        
        A = np.zeros((1,self.N))
        for i in range(self.N):
            A[0][i]=targets[i]
        b = np.zeros(1)

        # Call the the quadratic solver of cvxopt library.
        sol = cvxopt.solvers.qp(cvxopt.matrix(P), cvxopt.matrix(q), cvxopt.matrix(
            G), cvxopt.matrix(h), cvxopt.matrix(A), cvxopt.matrix(b))

        # Get the Lagrange multipliers out of the solution dictionary
        lambdas = np.array(sol['x'])

        # Find the (indices of the) support vectors, which are the vectors with non-zero Lagrange multipliers
        self.sv = np.where(lambdas > self.threshold)[0]
        self.nsupport = len(self.sv)
        print ("Number of support vectors = ", self.nsupport)

        # Keep the data corresponding to the support vectors
        self.X = X[self.sv, :]
        self.lambdas = lambdas[self.sv]
        self.targets = targets[self.sv]

        self.b = np.sum(self.targets)
        for n in range(self.nsupport):
            self.b -= np.sum(self.lambdas*self.targets *
                             np.reshape(self.K[self.sv[n], self.sv], (self.nsupport, 1)))
        self.b /= len(self.lambdas)

        if self.kernel == 'poly':
            def classify(Y, soft=False):
                K = (1. + 1./self.sigma*np.dot(Y, self.X.T))**self.degree

                self.y = np.zeros((np.shape(Y)[0], 1))
                for j in range(np.shape(Y)[0]):
                    for i in range(self.nsupport):
                        self.y[j] += self.lambdas[i]*self.targets[i]*K[j, i]
                    self.y[j] += self.b

                if soft:
                    return self.y
                else:
                    return np.sign(self.y)

        elif self.kernel == 'rbf':
            def classify(Y, soft=False):
                K = np.dot(Y, self.X.T)
                c = (1./self.sigma * np.sum(Y**2, axis=1)
                     * np.ones((1, np.shape(Y)[0]))).T
                c = np.dot(c, np.ones((1, np.shape(K)[1])))
                aa = np.dot(self.xsquared[self.sv],
                            np.ones((1, np.shape(K)[0]))).T
                K = K - 0.5*c - 0.5*aa
                K = np.exp(K/(2.*self.sigma**2))

                self.y = np.zeros((np.shape(Y)[0], 1))
                for j in range(np.shape(Y)[0]):
                    for i in range(self.nsupport):
                        self.y[j] += self.lambdas[i]*self.targets[i]*K[j, i]
                    self.y[j] += self.b

                if soft:
                    return self.y
                else:
                    return np.sign(self.y)
        else:
            print ("Error: Invalid kernel")
            return

        self.classify = classify

```

$\color{red}{Question:}$ How $b$ was computed? 

$\color{blue}{Answer:}$
$$
\begin{align}
\text{only the support vectors which lie on the separating hyperplanes are considered in both hard margin}\\
\text{and soft margin classifiers. the condition for support vector lying on separating hyperplane is same}\\
=>1= y_i*(w^T*x_i+b)=y_i(\sum_{j\in SV}^m \lambda_j*y_j*x_j^T*x_i+b) \\
\text{After putting}\quad  w =\sum_{j\in SV}^m \lambda_j*y_j*x_j^T \\
\text{Now, if we multiply by}\quad y_i\quad \text{on both sides, we get}\\
b= y_i-\sum_{j\in SV}^m \lambda_j*y_j*x_j^T*x_i \quad \text{after putting}\quad y_i^2=1 \\
\text{Here, Sv is the set of the support vectors, ie. all the vectors which lie on the separating hyperplane}\\
\text{All these vectors have a value of } \lambda_i > 0 \\
\end{align}
$$

# Test the Classifier
In the following, we will now test our classifier.


```python
from importlib import reload
import pylab as pl
import numpy as np

iris = np.loadtxt('iris_proc.data', delimiter=',')
imax = np.concatenate((iris.max(axis=0)*np.ones((1, 5)),
                       iris.min(axis=0)*np.ones((1, 5))), axis=0).max(axis=0)
target = -np.ones((np.shape(iris)[0], 3), dtype=float)
indices = np.where(iris[:, 4] == 0)
target[indices, 0] = 1.
indices = np.where(iris[:, 4] == 1)
target[indices, 1] = 1.
indices = np.where(iris[:, 4] == 2)
target[indices, 2] = 1.

train = iris[::2, 0:4]
traint = target[::2]
test = iris[1::2, 0:4]
testt = target[1::2]

```


```python
# Training the machines
output = np.zeros((np.shape(test)[0], 3))

# Train for the first set of train data
#svm0 = svm(kernel='linear')
#svm0 = svm(kernel='linear')
#svm0 = svm.svm(kernel='poly',C=0.1,degree=1)
svm0 = svm(kernel='rbf',C=0.2)
svm0.train_kernel_svm(train, np.reshape(traint[:, 0], (np.shape(train[:, :2])[0], 1)))
output[:, 0] = svm0.classify(test, soft=True).T

# Train for the second set of train data
#svm1 = svm(kernel='linear')
#svm1 = svm(kernel='linear')
#svm1 = svm(kernel='poly',degree=3)
svm1 = svm(kernel='rbf',C=0.2)
svm1.train_kernel_svm(train, np.reshape(traint[:, 1], (np.shape(train[:, :2])[0], 1)))
output[:, 1] = svm1.classify(test, soft=True).T

# Train for the third set of train data
#svm2 = svm(kernel='linear')
#svm2 = svm(kernel='linear')
#svm2 = svm(kernel='poly',C=0.1,degree=1)
svm2 = svm(kernel='rbf',C=0.2)
svm2.train_kernel_svm(train, np.reshape(traint[:, 2], (np.shape(train[:, :2])[0], 1)))
output[:, 2] = svm2.classify(test, soft=True).T
```

    [[ 40.26  34.06  39.96  34.77  29.77  45.14  34.01  51.12  46.22  49.91
       43.65  35.16  38.25  39.28  40.6   35.25  46.11  35.88  44.23  30.09
       39.03  31.33  44.22  43.05  44.07  83.29  83.48  73.5   75.23  74.82
       42.25  57.84  54.42  62.86  65.78  71.33  72.2   69.55  79.08  66.91
       51.66  57.58  60.66  78.84  58.86  57.81  57.84  57.98  60.23  42.47
       92.83  98.63  89.73  53.4   88.02  82.5   89.9   73.25  84.74 118.95
       95.63 116.02  92.68  72.56  84.57 103.42  85.    77.29  88.37  71.28
       91.62  70.55  94.52  74.55  84.45]]
    75
         pcost       dcost       gap    pres   dres
     0: -1.0648e-01 -1.5200e+01  2e+01  4e-15  2e-13
     1: -1.0777e-01 -3.4955e-01  2e-01  4e-15  7e-14
     2: -1.6785e-01 -2.0094e-01  3e-02  9e-16  2e-14
     3: -1.9911e-01 -1.9957e-01  5e-04  5e-15  1e-14
     4: -1.9943e-01 -1.9943e-01  5e-06  9e-15  9e-15
     5: -1.9943e-01 -1.9943e-01  5e-08  2e-15  4e-14
    Optimal solution found.
    Number of support vectors =  41
    [[ 40.26  34.06  39.96  34.77  29.77  45.14  34.01  51.12  46.22  49.91
       43.65  35.16  38.25  39.28  40.6   35.25  46.11  35.88  44.23  30.09
       39.03  31.33  44.22  43.05  44.07  83.29  83.48  73.5   75.23  74.82
       42.25  57.84  54.42  62.86  65.78  71.33  72.2   69.55  79.08  66.91
       51.66  57.58  60.66  78.84  58.86  57.81  57.84  57.98  60.23  42.47
       92.83  98.63  89.73  53.4   88.02  82.5   89.9   73.25  84.74 118.95
       95.63 116.02  92.68  72.56  84.57 103.42  85.    77.29  88.37  71.28
       91.62  70.55  94.52  74.55  84.45]]
    75
         pcost       dcost       gap    pres   dres
     0: -1.0667e-01 -1.5200e+01  2e+01  3e-16  3e-14
     1: -1.0796e-01 -3.5000e-01  2e-01  3e-16  6e-14
     2: -1.6844e-01 -2.0150e-01  3e-02  5e-16  8e-15
     3: -1.9968e-01 -2.0014e-01  5e-04  2e-16  7e-16
     4: -2.0000e-01 -2.0000e-01  5e-06  5e-16  2e-15
     5: -2.0000e-01 -2.0000e-01  5e-08  5e-16  2e-15
    Optimal solution found.
    Number of support vectors =  41
    [[ 40.26  34.06  39.96  34.77  29.77  45.14  34.01  51.12  46.22  49.91
       43.65  35.16  38.25  39.28  40.6   35.25  46.11  35.88  44.23  30.09
       39.03  31.33  44.22  43.05  44.07  83.29  83.48  73.5   75.23  74.82
       42.25  57.84  54.42  62.86  65.78  71.33  72.2   69.55  79.08  66.91
       51.66  57.58  60.66  78.84  58.86  57.81  57.84  57.98  60.23  42.47
       92.83  98.63  89.73  53.4   88.02  82.5   89.9   73.25  84.74 118.95
       95.63 116.02  92.68  72.56  84.57 103.42  85.    77.29  88.37  71.28
       91.62  70.55  94.52  74.55  84.45]]
    75
         pcost       dcost       gap    pres   dres
     0: -1.0667e-01 -1.5200e+01  2e+01  5e-16  2e-14
     1: -1.0796e-01 -3.5000e-01  2e-01  9e-16  3e-14
     2: -1.6844e-01 -2.0150e-01  3e-02  4e-16  3e-15
     3: -1.9968e-01 -2.0014e-01  5e-04  1e-16  2e-15
     4: -2.0000e-01 -2.0000e-01  5e-06  1e-15  5e-15
     5: -2.0000e-01 -2.0000e-01  5e-08  6e-16  6e-15
    Optimal solution found.
    Number of support vectors =  36



```python
# Make a decision about which class
# Pick the one with the largest margin
# print(np.shape(output))
print(output)
bestclass = np.argmax(output, axis=1)
print (bestclass)
print (iris[1::2, 4])
print("Misclassified locations:")
err = np.where(bestclass != iris[1::2, 4])[0]
print (err)
print (float(np.shape(testt)[0] - len(err)) /
       (np.shape(testt)[0]), "test accuracy")
```

    [[ 7.64303951 -0.46552986 -0.25061284]
     [ 7.48693593 -0.46488703 -0.24813203]
     [ 7.67708623 -0.45432574 -0.25545977]
     [ 8.1743781  -0.46696161 -0.25565272]
     [ 7.76839503 -0.46511391 -0.2530793 ]
     [ 7.95981174 -0.46462939 -0.25534217]
     [ 6.87078897 -0.46968904 -0.22783511]
     [ 6.94896151 -0.45629642 -0.22962122]
     [ 8.27977622 -0.46880181 -0.25344547]
     [ 8.23512701 -0.46707207 -0.25121927]
     [ 8.17294454 -0.46517706 -0.25405172]
     [ 7.49153053 -0.44991041 -0.26724505]
     [ 7.46865702 -0.4572944  -0.25877561]
     [ 8.2313703  -0.46588391 -0.2561192 ]
     [ 7.64379742 -0.46265932 -0.25380329]
     [ 7.86538329 -0.45861076 -0.25906144]
     [ 7.70640681 -0.46504815 -0.23550261]
     [ 8.04541424 -0.47270719 -0.2463286 ]
     [ 7.76839503 -0.46511391 -0.2530793 ]
     [ 8.1807916  -0.46604668 -0.25660821]
     [ 5.65368955 -0.45416122 -0.22616312]
     [ 7.71306536 -0.4565502  -0.26009382]
     [ 7.53754139 -0.46467014 -0.24992285]
     [ 7.69558131 -0.46897349 -0.24573056]
     [ 8.14660407 -0.46909753 -0.25298448]
     [-6.54884928 -0.3337532  -0.26253334]
     [-4.9465713  -0.26484916 -0.33963469]
     [-6.55158023 -0.28670356 -0.31214133]
     [-1.01225988 -0.30207626 -0.32574881]
     [-4.42702905 -0.26829406 -0.34080633]
     [-6.43772957 -0.28420675 -0.31729282]
     [-6.99566848 -0.32402663 -0.2715709 ]
     [-5.95746076 -0.33510489 -0.26109994]
     [-5.49785455 -0.25752072 -0.34845727]
     [-4.79960332 -0.25493147 -0.35424097]
     [-5.88868462 -0.26914773 -0.33519777]
     [-6.73337845 -0.31225386 -0.28378721]
     [-6.24093049 -0.32645904 -0.27004825]
     [-6.20220868 -0.40553471 -0.18683062]
     [-3.26618242 -0.26193575 -0.35909514]
     [-3.72506837 -0.25998992 -0.35440965]
     [-6.6430455  -0.37346987 -0.21942633]
     [-6.39576807 -0.32917847 -0.26817574]
     [-6.02061222 -0.30237439 -0.29413809]
     [-5.19557434 -0.26052274 -0.34536979]
     [-6.94327207 -0.31640611 -0.28015998]
     [-1.21295795 -0.29840027 -0.32769048]
     [-5.9187223  -0.26792369 -0.33657864]
     [-6.53515956 -0.28965678 -0.30999093]
     [-5.91103596 -0.26283135 -0.34204693]
     [-6.33072335 -0.38781704 -0.20481698]
     [-5.53458038 -0.45109787 -0.14027958]
     [-0.87472304 -0.48644853 -0.10473664]
     [-2.10602575 -0.48854247 -0.10253606]
     [-1.96382451 -0.50194287 -0.08911034]
     [-6.03787289 -0.43198339 -0.1596656 ]
     [-5.97229251 -0.37989815 -0.21287827]
     [-5.41119459 -0.46379517 -0.1276423 ]
     [-0.17846919 -0.47624144 -0.11504338]
     [-5.89522627 -0.35579222 -0.23722866]
     [-6.19055144 -0.36956166 -0.22391251]
     [-6.77327667 -0.37958797 -0.21334217]
     [-2.99894501 -0.48680924 -0.10423938]
     [-6.91193825 -0.37513357 -0.21828157]
     [-3.50909892 -0.47034506 -0.12078213]
     [-0.35432803 -0.47252709 -0.11875164]
     [-6.62190302 -0.38068478 -0.212053  ]
     [-1.50003462 -0.48992119 -0.10119726]
     [-5.69114221 -0.4488064  -0.14267131]
     [-4.87030309 -0.47149253 -0.11972763]
     [-4.93797036 -0.45966456 -0.13178586]
     [-3.65877783 -0.50361732 -0.08734993]
     [-5.28166672 -0.46202677 -0.12936965]
     [-6.08885147 -0.43756384 -0.15417067]
     [-6.58656055 -0.38849529 -0.2043457 ]]
    [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 1 1 1 1 1 2 2 1 1 1 2
     2 2 1 1 2 2 2 1 2 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
     2]
    [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
     0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
     1. 1. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.
     2. 2. 2.]
    Misclassified locations:
    [25 31 32 36 37 38 41 42 43 45]
    0.8666666666666667 test accuracy


## Further Questions

$\color{red}{Question}:$ The IRIS dataset has three classes. Explain by observing the code above how the two class SVM was modified for multiclass classification.

$\color{red}{Question}:$ Write mathematical expressions for the kernels defined above.

$\color{red}{Question}:$ Play with different Kernels. Which kernels (polynomial, RBF, or polynomial) give the best test accuracy?

Question1 from further questions:
$$
\begin{align}
\text{For this the output column in the training data was made into a 3 column form}\\
\text{So, the columns were created for three classes and all the samples were checked}\\
\text{for these 3 classes, the columns contained a one for true class label and -1 otherwise}\\
\text{So, binary classification was done 3 times using the training data and an N*3 size matrix was formed} \\
\text{After this, for each sample, the class corresponding to highest value was given as output}\\ \\ \\
\text{So, basically, in this case one to rest classification was done where the class with the highest value }\\
\text{is given as output for each sample of the data}\\
\end{align}
$$



Question2 from further questions
$$
\begin{align}
\text{Mathematical expressions for kernels} \\ \\
\text{The RBF kernel is given by:}\\
K(x,x')=exp (-\dfrac{||x-x'||^2}{2}) \\
=>exp(-\dfrac{1}{2}||x-x'||^2)=exp(X^TX'-\dfrac{1}{2}||X||^2-\dfrac{1}{2}||X'||^2)\\ \\
=> exp(X^TX')exp(-\dfrac{1}{2}||X||^2)exp(-\dfrac{1}{2}||X'||^2) \\ \\ \\
\text{the Polynomial kernel is given by}\\ 
=> K(x,y)=(x^Ty+c)^d \text{for a d degree polynomial} \\ 
\end{align}
$$

Question3 further questions:
$$
\begin{align}
\text{The RBF kernel gave an accuracy of 86.666 %} \\
\end{align}
$$


```python

```

c34cd2e9067dd064f7d48ecedad65ebf967f93a1
ipynb
Jupyter Notebook
Assignment-5-Question.ipynb
keshavbnsl102/TAO-SVM-assignment
d451f4a9a6e942ea8c6b2b9425c6e3c200cd225d
Apache-2.0
Assignment-5-Question.ipynb
keshavbnsl102/TAO-SVM-assignment
d451f4a9a6e942ea8c6b2b9425c6e3c200cd225d
Apache-2.0
Assignment-5-Question.ipynb
keshavbnsl102/TAO-SVM-assignment
d451f4a9a6e942ea8c6b2b9425c6e3c200cd225d
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES