# 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}:$ Put your answer here. 

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

$\color{blue}{Answer}:$ Put your answer here.

## 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:}$ Put your anser here (use latex)

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

$\color{blue}{Answer:}$ Put your answer here (use latex)

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

$\color{blue}{Answer}:$ Put your answer here. 

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

$\color{blue}{Answer}:$ Put your answer here.

## 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:}$ Put your anser here (use latex)

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

$\color{blue}{Answer:}$ Put your answer here (use latex)

# 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:}$ Write your answer in this cell.

## 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)

        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)

        # Assemble the matrices for the constraints   
        P = TODO
        q = TODO
        G = TODO
        h = TODO
        A = TODO
        b = TODO

        # 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:}$ Write your answer here.

# 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')
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')
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')
svm2.train_kernel_svm(train, np.reshape(traint[:, 2], (np.shape(train[:, :2])[0], 1)))
output[:, 2] = svm2.classify(test, soft=True).T
```


```python
# Make a decision about which class
# Pick the one with the largest margin
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")
```

## 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?

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