Morten Hjorth-Jensen
1376 lines
44 KiB
Plaintext
1376 lines
44 KiB
Plaintext
TITLE: Week 45: Decisions Trees, Random Forests, Bagging and Boosting
|
|
AUTHOR: Morten Hjorth-Jensen {copyright, 1999-present|CC BY-NC} at Department of Physics, University of Oslo & Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University
|
|
DATE: today
|
|
|
|
!split
|
|
===== Overview of week 45 =====
|
|
|
|
* Thursday: Boosting methods, froma AdaBoost to Gradient boosting
|
|
* "Video of lecture":"https://youtu.be/mK48PfCxgYk"
|
|
* Friday: Gradient boosting and discussion of Decision trees and ensemble methods. Wrapping up trees and start discussing Support Vector Machines
|
|
* "Video of lecture":"https://youtu.be/v8eJBFeZKuI"
|
|
|
|
!bblock Videos
|
|
o "Video on Decision trees":"https://www.youtube.com/watch?v=RmajweUFKvM&ab_channel=Simplilearn"
|
|
o "Video on boosting methods by Hastie":"https://www.youtube.com/watch?v=wPqtzj5VZus&ab_channel=H2O.ai".
|
|
o "Video on AdaBoost":"https://www.youtube.com/watch?v=LsK-xG1cLYA"
|
|
o "Video on Gradient boost, part 1, parts 2-4 follows":"https://www.youtube.com/watch?v=3CC4N4z3GJc"
|
|
!eblock
|
|
|
|
!bblock Reading
|
|
o "Hastie et al, chapter 10.1-10.10":"https://github.com/CompPhysics/MachineLearning/blob/master/doc/Textbooks/elementsstat.pdf". Geron's chapters 6 and 7 are also useful.
|
|
!eblock
|
|
|
|
|
|
!split
|
|
===== Brief code reminder from last wekk =====
|
|
|
|
!bc pycod
|
|
%matplotlib inline
|
|
# Common imports
|
|
from IPython.display import Image
|
|
from pydot import graph_from_dot_data
|
|
import pandas as pd
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
from sklearn.model_selection import train_test_split
|
|
from sklearn.tree import export_graphviz
|
|
from sklearn.preprocessing import StandardScaler, OneHotEncoder
|
|
from sklearn.compose import ColumnTransformer
|
|
from pydot import graph_from_dot_data
|
|
from sklearn.datasets import load_breast_cancer
|
|
from sklearn.svm import SVC
|
|
from sklearn.linear_model import LogisticRegression
|
|
from sklearn.tree import DecisionTreeClassifier
|
|
from sklearn.ensemble import RandomForestClassifier
|
|
from sklearn.preprocessing import LabelEncoder
|
|
from sklearn.model_selection import cross_validate
|
|
import scikitplot as skplt
|
|
from sklearn.preprocessing import StandardScaler
|
|
import os
|
|
|
|
# Where to save the figures and data files
|
|
PROJECT_ROOT_DIR = "Results"
|
|
FIGURE_ID = "Results/FigureFiles"
|
|
DATA_ID = "DataFiles/"
|
|
|
|
if not os.path.exists(PROJECT_ROOT_DIR):
|
|
os.mkdir(PROJECT_ROOT_DIR)
|
|
|
|
if not os.path.exists(FIGURE_ID):
|
|
os.makedirs(FIGURE_ID)
|
|
|
|
if not os.path.exists(DATA_ID):
|
|
os.makedirs(DATA_ID)
|
|
|
|
def image_path(fig_id):
|
|
return os.path.join(FIGURE_ID, fig_id)
|
|
|
|
def data_path(dat_id):
|
|
return os.path.join(DATA_ID, dat_id)
|
|
|
|
def save_fig(fig_id):
|
|
plt.savefig(image_path(fig_id) + ".png", format='png')
|
|
|
|
# Load the cancer data
|
|
cancer = load_breast_cancer()
|
|
|
|
X_train, X_test, y_train, y_test = train_test_split(cancer.data,cancer.target,random_state=0)
|
|
print(X_train.shape)
|
|
print(X_test.shape)
|
|
#Scale the data
|
|
scaler = StandardScaler()
|
|
scaler.fit(X_train)
|
|
X_train_scaled = scaler.transform(X_train)
|
|
X_test_scaled = scaler.transform(X_test)
|
|
#define methods
|
|
# Logistic Regression
|
|
logreg = LogisticRegression(solver='lbfgs')
|
|
logreg.fit(X_train_scaled, y_train)
|
|
print("Test set accuracy Logistic Regression with scaled data: {:.2f}".format(logreg.score(X_test_scaled,y_test)))
|
|
# Decision Trees
|
|
deep_tree_clf = DecisionTreeClassifier(max_depth=None)
|
|
deep_tree_clf.fit(X_train_scaled, y_train)
|
|
print("Test set accuracy with Decision Trees and scaled data: {:.2f}".format(deep_tree_clf.score(X_test_scaled,y_test)))
|
|
# Support Vector Machine
|
|
svm = SVC(gamma='auto', C=100)
|
|
svm.fit(X_train_scaled, y_train)
|
|
print("Test set accuracy SVM with scaled data: {:.2f}".format(logreg.score(X_test_scaled,y_test)))
|
|
# Random forests
|
|
#Instantiate the model with 500 trees and entropy as splitting criteria
|
|
Random_Forest_model = RandomForestClassifier(n_estimators=500,criterion="entropy")
|
|
Random_Forest_model.fit(X_train_scaled, y_train)
|
|
print("Test set accuracy with Random Forests and scaled data: {:.2f}".format(Random_Forest_model.score(X_test_scaled,y_test)))
|
|
|
|
|
|
|
|
y_pred = Random_Forest_model.predict(X_test_scaled)
|
|
skplt.metrics.plot_confusion_matrix(y_test, y_pred, normalize=True)
|
|
plt.show()
|
|
y_probas = Random_Forest_model.predict_proba(X_test_scaled)
|
|
skplt.metrics.plot_roc(y_test, y_probas)
|
|
plt.show()
|
|
skplt.metrics.plot_cumulative_gain(y_test, y_probas)
|
|
plt.show()
|
|
!ec
|
|
|
|
|
|
|
|
|
|
|
|
!split
|
|
===== Boosting, a Bird's Eye View =====
|
|
|
|
The basic idea is to combine weak classifiers in order to create a good
|
|
classifier. With a weak classifier we often intend a classifier which
|
|
produces results which are only slightly better than we would get by
|
|
random guesses.
|
|
|
|
This is done by applying in an iterative way a weak (or a standard
|
|
classifier like decision trees) to modify the data. In each iteration
|
|
we emphasize those observations which are misclassified by weighting
|
|
them with a factor.
|
|
|
|
|
|
!split
|
|
===== What is boosting? Additive Modelling/Iterative Fitting =====
|
|
|
|
Boosting is a way of fitting an additive expansion in a set of
|
|
elementary basis functions like for example some simple polynomials.
|
|
Assume for example that we have a function
|
|
!bt
|
|
\[
|
|
f_M(x) = \sum_{i=1}^M \beta_m b(x;\gamma_m),
|
|
\]
|
|
!et
|
|
|
|
where $\beta_m$ are the expansion parameters to be determined in a
|
|
minimization process and $b(x;\gamma_m)$ are some simple functions of
|
|
the multivariable parameter $x$ which is characterized by the
|
|
parameters $\gamma_m$.
|
|
|
|
As an example, consider the Sigmoid function we used in logistic
|
|
regression. In that case, we can translate the function
|
|
$b(x;\gamma_m)$ into the Sigmoid function
|
|
|
|
|
|
!bt
|
|
\[
|
|
\sigma(t) = \frac{1}{1+\exp{(-t)}},
|
|
\]
|
|
!et
|
|
|
|
where $t=\gamma_0+\gamma_1 x$ and the parameters $\gamma_0$ and
|
|
$\gamma_1$ were determined by the Logistic Regression fitting
|
|
algorithm.
|
|
|
|
As another example, consider the cost function we defined for linear regression
|
|
!bt
|
|
\[
|
|
C(\bm{y},\bm{f}) = \frac{1}{n} \sum_{i=0}^{n-1}(y_i-f(x_i))^2.
|
|
\]
|
|
!et
|
|
|
|
In this case the function $f(x)$ was replaced by the design matrix
|
|
$\bm{X}$ and the unknown linear regression parameters $\bm{\beta}$,
|
|
that is $\bm{f}=\bm{X}\bm{\beta}$. In linear regression we can
|
|
simply invert a matrix and obtain the parameters $\beta$ by
|
|
|
|
!bt
|
|
\[
|
|
\bm{\beta}=\left(\bm{X}^T\bm{X}\right)^{-1}\bm{X}^T\bm{y}.
|
|
\]
|
|
!et
|
|
|
|
In iterative fitting or additive modeling, we minimize the cost function with respect to the parameters $\beta_m$ and $\gamma_m$.
|
|
|
|
|
|
!split
|
|
===== Iterative Fitting, Regression and Squared-error Cost Function =====
|
|
|
|
The way we proceed is as follows (here we specialize to the squared-error cost function)
|
|
|
|
o Establish a cost function, here ${\cal C}(\bm{y},\bm{f}) = \frac{1}{n} \sum_{i=0}^{n-1}(y_i-f_M(x_i))^2$ with $f_M(x) = \sum_{i=1}^M \beta_m b(x;\gamma_m)$.
|
|
o Initialize with a guess $f_0(x)$. It could be one or even zero or some random numbers.
|
|
o For $m=1:M$
|
|
o minimize $\sum_{i=0}^{n-1}(y_i-f_{m-1}(x_i)-\beta b(x;\gamma))^2$ wrt $\gamma$ and $\beta$
|
|
o This gives the optimal values $\beta_m$ and $\gamma_m$
|
|
o Determine then the new values $f_m(x)=f_{m-1}(x) +\beta_m b(x;\gamma_m)$
|
|
|
|
We could use any of the algorithms we have discussed till now. If we
|
|
use trees, $\gamma$ parameterizes the split variables and split points
|
|
at the internal nodes, and the predictions at the terminal nodes.
|
|
|
|
|
|
!split
|
|
===== Squared-Error Example and Iterative Fitting =====
|
|
|
|
To better understand what happens, let us develop the steps for the iterative fitting using the above squared error function.
|
|
|
|
For simplicity we assume also that our functions $b(x;\gamma)=1+\gamma x$.
|
|
|
|
This means that for every iteration $m$, we need to optimize
|
|
|
|
!bt
|
|
\[
|
|
(\beta_m,\gamma_m) = \mathrm{argmin}_{\beta,\lambda}\hspace{0.1cm} \sum_{i=0}^{n-1}(y_i-f_{m-1}(x_i)-\beta b(x;\gamma))^2=\sum_{i=0}^{n-1}(y_i-f_{m-1}(x_i)-\beta(1+\gamma x_i))^2.
|
|
\]
|
|
!et
|
|
|
|
We start our iteration by simply setting $f_0(x)=0$.
|
|
Taking the derivatives with respect to $\beta$ and $\gamma$ we obtain
|
|
!bt
|
|
\[
|
|
\frac{\partial {\cal C}}{\partial \beta} = -2\sum_{i}(1+\gamma x_i)(y_i-\beta(1+\gamma x_i))=0,
|
|
\]
|
|
!et
|
|
and
|
|
!bt
|
|
\[
|
|
\frac{\partial {\cal C}}{\partial \gamma} =-2\sum_{i}\beta x_i(y_i-\beta(1+\gamma x_i))=0.
|
|
\]
|
|
!et
|
|
We can then rewrite these equations as (defining $\bm{w}=\bm{e}+\gamma \bm{x})$ with $\bm{e}$ being the unit vector)
|
|
!bt
|
|
\[
|
|
\gamma \bm{w}^T(\bm{y}-\beta\gamma \bm{w})=0,
|
|
\]
|
|
!et
|
|
which gives us $\beta = \bm{w}^T\bm{y}/(\bm{w}^T\bm{w})$. Similarly we have
|
|
!bt
|
|
\[
|
|
\beta\gamma \bm{x}^T(\bm{y}-\beta(1+\gamma \bm{x}))=0,
|
|
\]
|
|
!et
|
|
|
|
which leads to $\gamma =(\bm{x}^T\bm{y}-\beta\bm{x}^T\bm{e})/(\beta\bm{x}^T\bm{x})$. Inserting
|
|
for $\beta$ gives us an equation for $\gamma$. This is a non-linear equation in the unknown $\gamma$ and has to be solved numerically.
|
|
|
|
The solution to these two equations gives us in turn $\beta_1$ and $\gamma_1$ leading to the new expression for $f_1(x)$ as
|
|
$f_1(x) = \beta_1(1+\gamma_1x)$. Doing this $M$ times results in our final estimate for the function $f$.
|
|
|
|
|
|
|
|
!split
|
|
===== Iterative Fitting, Classification and AdaBoost =====
|
|
|
|
Let us consider a binary classification problem with two outcomes $y_i \in \{-1,1\}$ and $i=0,1,2,\dots,n-1$ as our set of
|
|
observations. We define a classification function $G(x)$ which produces a prediction taking one or the other of the two values
|
|
$\{-1,1\}$.
|
|
|
|
The error rate of the training sample is then
|
|
|
|
!bt
|
|
\[
|
|
\mathrm{\overline{err}}=\frac{1}{n} \sum_{i=0}^{n-1} I(y_i\ne G(x_i)).
|
|
\]
|
|
!et
|
|
|
|
The iterative procedure starts with defining a weak classifier whose
|
|
error rate is barely better than random guessing. The iterative
|
|
procedure in boosting is to sequentially apply a weak
|
|
classification algorithm to repeatedly modified versions of the data
|
|
producing a sequence of weak classifiers $G_m(x)$.
|
|
|
|
Here we will express our function $f(x)$ in terms of $G(x)$. That is
|
|
!bt
|
|
\[
|
|
f_M(x) = \sum_{i=1}^M \beta_m b(x;\gamma_m),
|
|
\]
|
|
!et
|
|
will be a function of
|
|
!bt
|
|
\[
|
|
G_M(x) = \mathrm{sign} \sum_{i=1}^M \alpha_m G_m(x).
|
|
\]
|
|
!et
|
|
|
|
|
|
|
|
!split
|
|
===== Adaptive Boosting, AdaBoost =====
|
|
|
|
In our iterative procedure we define thus
|
|
!bt
|
|
\[
|
|
f_m(x) = f_{m-1}(x)+\beta_mG_m(x).
|
|
\]
|
|
!et
|
|
|
|
The simplest possible cost function which leads (also simple from a computational point of view) to the AdaBoost algorithm is the
|
|
exponential cost/loss function defined as
|
|
!bt
|
|
\[
|
|
C(\bm{y},\bm{f}) = \sum_{i=0}^{n-1}\exp{(-y_i(f_{m-1}(x_i)+\beta G(x_i))}.
|
|
\]
|
|
!et
|
|
|
|
We optimize $\beta$ and $G$ for each value of $m=1:M$ as we did in the regression case.
|
|
This is normally done in two steps. Let us however first rewrite the cost function as
|
|
|
|
!bt
|
|
\[
|
|
C(\bm{y},\bm{f}) = \sum_{i=0}^{n-1}w_i^{m}\exp{(-y_i\beta G(x_i))},
|
|
\]
|
|
!et
|
|
where we have defined $w_i^m= \exp{(-y_if_{m-1}(x_i))}$.
|
|
|
|
!split
|
|
===== Building up AdaBoost =====
|
|
|
|
First, for any $\beta > 0$, we optimize $G$ by setting
|
|
!bt
|
|
\[
|
|
G_m(x) = \mathrm{sign} \sum_{i=0}^{n-1} w_i^m I(y_i \ne G_(x_i)),
|
|
\]
|
|
!et
|
|
which is the classifier that minimizes the weighted error rate in predicting $y$.
|
|
|
|
We can do this by rewriting
|
|
!bt
|
|
\[
|
|
\exp{-(\beta)}\sum_{y_i=G(x_i)}w_i^m+\exp{(\beta)}\sum_{y_i\ne G(x_i)}w_i^m,
|
|
\]
|
|
!et
|
|
which can be rewritten as
|
|
!bt
|
|
\[
|
|
(\exp{(\beta)}-\exp{-(\beta)})\sum_{i=0}^{n-1}w_i^mI(y_i\ne G(x_i))+\exp{(-\beta)}\sum_{i=0}^{n-1}w_i^m=0,
|
|
\]
|
|
!et
|
|
which leads to
|
|
!bt
|
|
\[
|
|
\beta_m = \frac{1}{2}\log{\frac{1-\mathrm{\overline{err}}}{\mathrm{\overline{err}}}},
|
|
\]
|
|
!et
|
|
where we have redefined the error as
|
|
!bt
|
|
\[
|
|
\mathrm{\overline{err}}_m=\frac{1}{n}\frac{\sum_{i=0}^{n-1}w_i^mI(y_i\ne G(x_i)}{\sum_{i=0}^{n-1}w_i^m},
|
|
\]
|
|
!et
|
|
which leads to an update of
|
|
!bt
|
|
\[
|
|
f_m(x) = f_{m-1}(x) +\beta_m G_m(x).
|
|
\]
|
|
!et
|
|
This leads to the new weights
|
|
!bt
|
|
\[
|
|
w_i^{m+1} = w_i^m \exp{(-y_i\beta_m G_m(x_i))}
|
|
\]
|
|
!et
|
|
|
|
|
|
!split
|
|
===== Adaptive boosting: AdaBoost, Basic Algorithm =====
|
|
|
|
The algorithm here is rather straightforward. Assume that our weak
|
|
classifier is a decision tree and we consider a binary set of outputs
|
|
with $y_i \in \{-1,1\}$ and $i=0,1,2,\dots,n-1$ as our set of
|
|
observations. Our design matrix is given in terms of the
|
|
feature/predictor vectors
|
|
$\bm{X}=[\bm{x}_0\bm{x}_1\dots\bm{x}_{p-1}]$. Finally, we define also a
|
|
classifier determined by our data via a function $G(x)$. This function tells us how well we are able to classify our outputs/targets $\bm{y}$.
|
|
|
|
We have already defined the misclassification error $\mathrm{err}$ as
|
|
!bt
|
|
\[
|
|
\mathrm{err}=\frac{1}{n}\sum_{i=0}^{n-1}I(y_i\ne G(x_i)),
|
|
\]
|
|
!et
|
|
where the function $I()$ is one if we misclassify and zero if we classify correctly.
|
|
|
|
!split
|
|
===== Basic Steps of AdaBoost =====
|
|
|
|
With the above definitions we are now ready to set up the algorithm for AdaBoost.
|
|
The basic idea is to set up weights which will be used to scale the correctly classified and the misclassified cases.
|
|
o We start by initializing all weights to $w_i = 1/n$, with $i=0,1,2,\dots n-1$. It is easy to see that we must have $\sum_{i=0}^{n-1}w_i = 1$.
|
|
o We rewrite the misclassification error as
|
|
!bt
|
|
\[
|
|
\mathrm{\overline{err}}_m=\frac{\sum_{i=0}^{n-1}w_i^m I(y_i\ne G(x_i))}{\sum_{i=0}^{n-1}w_i},
|
|
\]
|
|
!et
|
|
o Then we start looping over all attempts at classifying, namely we start an iterative process for $m=1:M$, where $M$ is the final number of classifications. Our given classifier could for example be a plain decision tree.
|
|
o Fit then a given classifier to the training set using the weights $w_i$.
|
|
o Compute then $\mathrm{err}$ and figure out which events are classified properly and which are classified wrongly.
|
|
o Define a quantity $\alpha_{m} = \log{(1-\mathrm{\overline{err}}_m)/\mathrm{\overline{err}}_m}$
|
|
o Set the new weights to $w_i = w_i\times \exp{(\alpha_m I(y_i\ne G(x_i)}$.
|
|
o Compute the new classifier $G(x)= \sum_{i=0}^{n-1}\alpha_m I(y_i\ne G(x_i)$.
|
|
|
|
For the iterations with $m \le 2$ the weights are modified
|
|
individually at each steps. The observations which were misclassified
|
|
at iteration $m-1$ have a weight which is larger than those which were
|
|
classified properly. As this proceeds, the observations which were
|
|
difficult to classifiy correctly are given a larger influence. Each
|
|
new classification step $m$ is then forced to concentrate on those
|
|
observations that are missed in the previous iterations.
|
|
|
|
|
|
|
|
!split
|
|
===== AdaBoost Examples =====
|
|
|
|
Using _Scikit-Learn_ it is easy to apply the adaptive boosting algorithm, as done here.
|
|
|
|
!bc pycod
|
|
from sklearn.ensemble import AdaBoostClassifier
|
|
|
|
ada_clf = AdaBoostClassifier(
|
|
DecisionTreeClassifier(max_depth=2), n_estimators=200,
|
|
algorithm="SAMME.R", learning_rate=0.01, random_state=42)
|
|
ada_clf.fit(X_train, y_train)
|
|
y_pred = ada_clf.predict(X_test)
|
|
skplt.metrics.plot_confusion_matrix(y_test, y_pred, normalize=True)
|
|
plt.show()
|
|
y_probas = ada_clf.predict_proba(X_test)
|
|
skplt.metrics.plot_roc(y_test, y_probas)
|
|
plt.show()
|
|
skplt.metrics.plot_cumulative_gain(y_test, y_probas)
|
|
plt.show()
|
|
!ec
|
|
|
|
|
|
|
|
!split
|
|
===== Gradient boosting: Basics with Steepest Descent/Functional Gradient Descent =====
|
|
|
|
Gradient boosting is again a similar technique to Adaptive boosting,
|
|
it combines so-called weak classifiers or regressors into a strong
|
|
method via a series of iterations.
|
|
|
|
In order to understand the method, let us illustrate its basics by
|
|
bringing back the essential steps in linear regression, where our cost
|
|
function was the least squares function.
|
|
|
|
!split
|
|
===== The Squared-Error again! Steepest Descent =====
|
|
|
|
We start again with our cost function ${\cal C}(\bm{y}m\bm{f})=\sum_{i=0}^{n-1}{\cal L}(y_i, f(x_i))$ where we want to minimize
|
|
This means that for every iteration, we need to optimize
|
|
|
|
!bt
|
|
\[
|
|
(\hat{\bm{f}}) = \mathrm{argmin}_{\bm{f}}\hspace{0.1cm} \sum_{i=0}^{n-1}(y_i-f(x_i))^2.
|
|
\]
|
|
!et
|
|
|
|
We define a real function $h_m(x)$ that defines our final function $f_M(x)$ as
|
|
!bt
|
|
\[
|
|
f_M(x) = \sum_{m=0}^M h_m(x).
|
|
\]
|
|
!et
|
|
|
|
In the steepest decent approach we approximate $h_m(x) = -\rho_m g_m(x)$, where $\rho_m$ is a scalar and $g_m(x)$ the gradient defined as
|
|
!bt
|
|
\[
|
|
g_m(x_i) = \left[ \frac{\partial {\cal L}(y_i, f(x_i))}{\partial f(x_i)}\right]_{f(x_i)=f_{m-1}(x_i)}.
|
|
\]
|
|
!et
|
|
|
|
With the new gradient we can update $f_m(x) = f_{m-1}(x) -\rho_m g_m(x)$. Using the above squared-error function we see that
|
|
the gradient is $g_m(x_i) = -2(y_i-f(x_i))$.
|
|
|
|
Choosing $f_0(x)=0$ we obtain $g_m(x) = -2y_i$ and inserting this into the minimization problem for the cost function we have
|
|
!bt
|
|
\[
|
|
(\rho_1) = \mathrm{argmin}_{\rho}\hspace{0.1cm} \sum_{i=0}^{n-1}(y_i+2\rho y_i)^2.
|
|
\]
|
|
!et
|
|
|
|
!split
|
|
===== Steepest Descent Example =====
|
|
|
|
Optimizing with respect to $\rho$ we obtain (taking the derivative) that $\rho_1 = -1/2$. We have then that
|
|
!bt
|
|
\[
|
|
f_1(x) = f_{0}(x) -\rho_1 g_1(x)=-y_i.
|
|
\]
|
|
!et
|
|
We can then proceed and compute
|
|
!bt
|
|
\[
|
|
g_2(x_i) = \left[ \frac{\partial {\cal L}(y_i, f(x_i))}{\partial f(x_i)}\right]_{f(x_i)=f_{1}(x_i)=y_i}=-4y_i,
|
|
\]
|
|
!et
|
|
and find a new value for $\rho_2=-1/2$ and continue till we have reached $m=M$. We can modify the steepest descent method, or steepest boosting, by introducing what is called _gradient boosting_.
|
|
|
|
!split
|
|
===== Gradient Boosting, algorithm =====
|
|
|
|
Steepest descent is however not much used, since it only optimizes $f$ at a fixed set of $n$ points,
|
|
so we do not learn a function that can generalize. However, we can modify the algorithm by
|
|
fitting a weak learner to approximate the negative gradient signal.
|
|
|
|
Suppose we have a cost function $C(f)=\sum_{i=0}^{n-1}L(y_i, f(x_i))$ where $y_i$ is our target and $f(x_i)$ the function which is meant to model $y_i$. The above cost function could be our standard squared-error function
|
|
!bt
|
|
\[
|
|
C(\bm{y},\bm{f})=\sum_{i=0}^{n-1}(y_i-f(x_i))^2.
|
|
\]
|
|
!et
|
|
|
|
The way we proceed in an iterative fashion is to
|
|
o Initialize our estimate $f_0(x)$.
|
|
o For $m=1:M$, we
|
|
o compute the negative gradient vector $\bm{u}_m = -\partial C(\bm{y},\bm{f})/\partial \bm{f}(x)$ at $f(x) = f_{m-1}(x)$;
|
|
o fit the so-called base-learner to the negative gradient $h_m(u_m,x)$;
|
|
o update the estimate $f_m(x) = f_{m-1}(x)+h_m(u_m,x)$;
|
|
o The final estimate is then $f_M(x) = \sum_{m=1}^M h_m(u_m,x)$.
|
|
|
|
|
|
|
|
!split
|
|
===== Gradient Boosting, Examples of Regression =====
|
|
!bc pycod
|
|
import matplotlib.pyplot as plt
|
|
import numpy as np
|
|
from sklearn.model_selection import train_test_split
|
|
from sklearn.ensemble import GradientBoostingRegressor
|
|
import scikitplot as skplt
|
|
from sklearn.metrics import mean_squared_error
|
|
|
|
n = 100
|
|
maxdegree = 6
|
|
|
|
# Make data set.
|
|
x = np.linspace(-3, 3, n).reshape(-1, 1)
|
|
y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)
|
|
|
|
error = np.zeros(maxdegree)
|
|
bias = np.zeros(maxdegree)
|
|
variance = np.zeros(maxdegree)
|
|
polydegree = np.zeros(maxdegree)
|
|
X_train, X_test, y_train, y_test = train_test_split(x, y, test_size=0.2)
|
|
|
|
for degree in range(1,maxdegree):
|
|
model = GradientBoostingRegressor(max_depth=degree, n_estimators=100, learning_rate=1.0)
|
|
model.fit(X_train,y_train)
|
|
y_pred = model.predict(X_test)
|
|
polydegree[degree] = degree
|
|
error[degree] = np.mean( np.mean((y_test - y_pred)**2) )
|
|
bias[degree] = np.mean( (y_test - np.mean(y_pred))**2 )
|
|
variance[degree] = np.mean( np.var(y_pred) )
|
|
print('Max depth:', degree)
|
|
print('Error:', error[degree])
|
|
print('Bias^2:', bias[degree])
|
|
print('Var:', variance[degree])
|
|
print('{} >= {} + {} = {}'.format(error[degree], bias[degree], variance[degree], bias[degree]+variance[degree]))
|
|
|
|
plt.xlim(1,maxdegree-1)
|
|
plt.plot(polydegree, error, label='Error')
|
|
plt.plot(polydegree, bias, label='bias')
|
|
plt.plot(polydegree, variance, label='Variance')
|
|
plt.legend()
|
|
save_fig("gdregression")
|
|
plt.show()
|
|
!ec
|
|
|
|
|
|
!split
|
|
===== Gradient Boosting, Classification Example =====
|
|
!bc pycod
|
|
import matplotlib.pyplot as plt
|
|
import numpy as np
|
|
from sklearn.model_selection import train_test_split
|
|
from sklearn.datasets import load_breast_cancer
|
|
import scikitplot as skplt
|
|
from sklearn.ensemble import GradientBoostingClassifier
|
|
from sklearn.model_selection import cross_validate
|
|
|
|
# Load the data
|
|
cancer = load_breast_cancer()
|
|
|
|
X_train, X_test, y_train, y_test = train_test_split(cancer.data,cancer.target,random_state=0)
|
|
print(X_train.shape)
|
|
print(X_test.shape)
|
|
#now scale the data
|
|
from sklearn.preprocessing import StandardScaler
|
|
scaler = StandardScaler()
|
|
scaler.fit(X_train)
|
|
X_train_scaled = scaler.transform(X_train)
|
|
X_test_scaled = scaler.transform(X_test)
|
|
|
|
gd_clf = GradientBoostingClassifier(max_depth=3, n_estimators=100, learning_rate=1.0)
|
|
gd_clf.fit(X_train_scaled, y_train)
|
|
#Cross validation
|
|
accuracy = cross_validate(gd_clf,X_test_scaled,y_test,cv=10)['test_score']
|
|
print(accuracy)
|
|
print("Test set accuracy with Gradient boosting and scaled data: {:.2f}".format(gd_clf.score(X_test_scaled,y_test)))
|
|
|
|
import scikitplot as skplt
|
|
y_pred = gd_clf.predict(X_test_scaled)
|
|
skplt.metrics.plot_confusion_matrix(y_test, y_pred, normalize=True)
|
|
save_fig("gdclassiffierconfusion")
|
|
plt.show()
|
|
y_probas = gd_clf.predict_proba(X_test_scaled)
|
|
skplt.metrics.plot_roc(y_test, y_probas)
|
|
save_fig("gdclassiffierroc")
|
|
plt.show()
|
|
skplt.metrics.plot_cumulative_gain(y_test, y_probas)
|
|
save_fig("gdclassiffiercgain")
|
|
plt.show()
|
|
!ec
|
|
|
|
|
|
!split
|
|
===== XGBoost: Extreme Gradient Boosting =====
|
|
|
|
|
|
"XGBoost":"https://github.com/dmlc/xgboost" or Extreme Gradient
|
|
Boosting, is an optimized distributed gradient boosting library
|
|
designed to be highly efficient, flexible and portable. It implements
|
|
machine learning algorithms under the Gradient Boosting
|
|
framework. XGBoost provides a parallel tree boosting that solve many
|
|
data science problems in a fast and accurate way. See the "article by Chen and Guestrin":"https://arxiv.org/abs/1603.02754".
|
|
|
|
The authors design and build a highly scalable end-to-end tree
|
|
boosting system. It has a theoretically justified weighted quantile
|
|
sketch for efficient proposal calculation. It introduces a novel sparsity-aware algorithm for parallel tree learning and an effective cache-aware block structure for out-of-core tree learning.
|
|
|
|
It is now the algorithm which wins essentially all ML competitions!!!
|
|
|
|
!split
|
|
===== Regression Case =====
|
|
|
|
!bc pycod
|
|
import matplotlib.pyplot as plt
|
|
import numpy as np
|
|
from sklearn.model_selection import train_test_split
|
|
import xgboost as xgb
|
|
import scikitplot as skplt
|
|
from sklearn.metrics import mean_squared_error
|
|
|
|
n = 100
|
|
maxdegree = 6
|
|
|
|
# Make data set.
|
|
x = np.linspace(-3, 3, n).reshape(-1, 1)
|
|
y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)
|
|
|
|
error = np.zeros(maxdegree)
|
|
bias = np.zeros(maxdegree)
|
|
variance = np.zeros(maxdegree)
|
|
polydegree = np.zeros(maxdegree)
|
|
X_train, X_test, y_train, y_test = train_test_split(x, y, test_size=0.2)
|
|
|
|
for degree in range(maxdegree):
|
|
model = xgb.XGBRegressor(objective ='reg:squarederror', colsaobjective ='reg:squarederror', colsample_bytree = 0.3, learning_rate = 0.1,max_depth = degree, alpha = 10, n_estimators = 200)
|
|
|
|
model.fit(X_train,y_train)
|
|
y_pred = model.predict(X_test)
|
|
polydegree[degree] = degree
|
|
error[degree] = np.mean( np.mean((y_test - y_pred)**2) )
|
|
bias[degree] = np.mean( (y_test - np.mean(y_pred))**2 )
|
|
variance[degree] = np.mean( np.var(y_pred) )
|
|
print('Max depth:', degree)
|
|
print('Error:', error[degree])
|
|
print('Bias^2:', bias[degree])
|
|
print('Var:', variance[degree])
|
|
print('{} >= {} + {} = {}'.format(error[degree], bias[degree], variance[degree], bias[degree]+variance[degree]))
|
|
|
|
plt.xlim(1,maxdegree-1)
|
|
plt.plot(polydegree, error, label='Error')
|
|
plt.plot(polydegree, bias, label='bias')
|
|
plt.plot(polydegree, variance, label='Variance')
|
|
plt.legend()
|
|
plt.show()
|
|
|
|
!ec
|
|
|
|
!split
|
|
===== Xgboost on the Cancer Data =====
|
|
|
|
As you will see from the confusion matrix below, XGBoots does an excellent job on the Wisconsin cancer data and outperforms essentially all agorithms we have discussed till now.
|
|
!bc pycod
|
|
|
|
import matplotlib.pyplot as plt
|
|
import numpy as np
|
|
from sklearn.model_selection import train_test_split
|
|
from sklearn.datasets import load_breast_cancer
|
|
from sklearn.preprocessing import LabelEncoder
|
|
from sklearn.model_selection import cross_validate
|
|
import scikitplot as skplt
|
|
import xgboost as xgb
|
|
# Load the data
|
|
cancer = load_breast_cancer()
|
|
|
|
X_train, X_test, y_train, y_test = train_test_split(cancer.data,cancer.target,random_state=0)
|
|
print(X_train.shape)
|
|
print(X_test.shape)
|
|
#now scale the data
|
|
from sklearn.preprocessing import StandardScaler
|
|
scaler = StandardScaler()
|
|
scaler.fit(X_train)
|
|
X_train_scaled = scaler.transform(X_train)
|
|
X_test_scaled = scaler.transform(X_test)
|
|
|
|
xg_clf = xgb.XGBClassifier()
|
|
xg_clf.fit(X_train_scaled,y_train)
|
|
|
|
y_test = xg_clf.predict(X_test_scaled)
|
|
|
|
print("Test set accuracy with Gradient Boosting and scaled data: {:.2f}".format(xg_clf.score(X_test_scaled,y_test)))
|
|
|
|
import scikitplot as skplt
|
|
y_pred = xg_clf.predict(X_test_scaled)
|
|
skplt.metrics.plot_confusion_matrix(y_test, y_pred, normalize=True)
|
|
save_fig("xdclassiffierconfusion")
|
|
plt.show()
|
|
y_probas = xg_clf.predict_proba(X_test_scaled)
|
|
skplt.metrics.plot_roc(y_test, y_probas)
|
|
save_fig("xdclassiffierroc")
|
|
plt.show()
|
|
skplt.metrics.plot_cumulative_gain(y_test, y_probas)
|
|
save_fig("gdclassiffiercgain")
|
|
plt.show()
|
|
|
|
|
|
xgb.plot_tree(xg_clf,num_trees=0)
|
|
plt.rcParams['figure.figsize'] = [50, 10]
|
|
save_fig("xgtree")
|
|
plt.show()
|
|
|
|
xgb.plot_importance(xg_clf)
|
|
plt.rcParams['figure.figsize'] = [5, 5]
|
|
save_fig("xgparams")
|
|
plt.show()
|
|
|
|
!ec
|
|
|
|
|
|
|
|
!split
|
|
===== Support Vector Machines, overarching aims =====
|
|
|
|
A Support Vector Machine (SVM) is a very powerful and versatile
|
|
Machine Learning method, capable of performing linear or nonlinear
|
|
classification, regression, and even outlier detection. It is one of
|
|
the most popular models in Machine Learning, and anyone interested in
|
|
Machine Learning should have it in their toolbox. SVMs are
|
|
particularly well suited for classification of complex but small-sized or
|
|
medium-sized datasets.
|
|
|
|
The case with two well-separated classes only can be understood in an
|
|
intuitive way in terms of lines in a two-dimensional space separating
|
|
the two classes (see figure below).
|
|
|
|
The basic mathematics behind the SVM is however less familiar to most of us.
|
|
It relies on the definition of hyperplanes and the
|
|
definition of a _margin_ which separates classes (in case of
|
|
classification problems) of variables. It is also used for regression
|
|
problems.
|
|
|
|
With SVMs we distinguish between hard margin and soft margins. The
|
|
latter introduces a so-called softening parameter to be discussed
|
|
below. We distinguish also between linear and non-linear
|
|
approaches. The latter are the most frequent ones since it is rather
|
|
unlikely that we can separate classes easily by say straight lines.
|
|
|
|
!split
|
|
===== Hyperplanes and all that =====
|
|
|
|
The theory behind support vector machines (SVM hereafter) is based on
|
|
the mathematical description of so-called hyperplanes. Let us start
|
|
with a two-dimensional case. This will also allow us to introduce our
|
|
first SVM examples. These will be tailored to the case of two specific
|
|
classes, as displayed in the figure here based on the usage of the petal data.
|
|
|
|
We assume here that our data set can be well separated into two
|
|
domains, where a straight line does the job in the separating the two
|
|
classes. Here the two classes are represented by either squares or
|
|
circles.
|
|
!bc pycod
|
|
from sklearn import datasets
|
|
from sklearn.svm import SVC, LinearSVC
|
|
from sklearn.linear_model import SGDClassifier
|
|
from sklearn.preprocessing import StandardScaler
|
|
import matplotlib
|
|
import matplotlib.pyplot as plt
|
|
plt.rcParams['axes.labelsize'] = 14
|
|
plt.rcParams['xtick.labelsize'] = 12
|
|
plt.rcParams['ytick.labelsize'] = 12
|
|
|
|
|
|
iris = datasets.load_iris()
|
|
X = iris["data"][:, (2, 3)] # petal length, petal width
|
|
y = iris["target"]
|
|
|
|
setosa_or_versicolor = (y == 0) | (y == 1)
|
|
X = X[setosa_or_versicolor]
|
|
y = y[setosa_or_versicolor]
|
|
|
|
|
|
|
|
C = 5
|
|
alpha = 1 / (C * len(X))
|
|
|
|
lin_clf = LinearSVC(loss="hinge", C=C, random_state=42)
|
|
svm_clf = SVC(kernel="linear", C=C)
|
|
sgd_clf = SGDClassifier(loss="hinge", learning_rate="constant", eta0=0.001, alpha=alpha,
|
|
max_iter=100000, random_state=42)
|
|
|
|
scaler = StandardScaler()
|
|
X_scaled = scaler.fit_transform(X)
|
|
|
|
lin_clf.fit(X_scaled, y)
|
|
svm_clf.fit(X_scaled, y)
|
|
sgd_clf.fit(X_scaled, y)
|
|
|
|
print("LinearSVC: ", lin_clf.intercept_, lin_clf.coef_)
|
|
print("SVC: ", svm_clf.intercept_, svm_clf.coef_)
|
|
print("SGDClassifier(alpha={:.5f}):".format(sgd_clf.alpha), sgd_clf.intercept_, sgd_clf.coef_)
|
|
|
|
# Compute the slope and bias of each decision boundary
|
|
w1 = -lin_clf.coef_[0, 0]/lin_clf.coef_[0, 1]
|
|
b1 = -lin_clf.intercept_[0]/lin_clf.coef_[0, 1]
|
|
w2 = -svm_clf.coef_[0, 0]/svm_clf.coef_[0, 1]
|
|
b2 = -svm_clf.intercept_[0]/svm_clf.coef_[0, 1]
|
|
w3 = -sgd_clf.coef_[0, 0]/sgd_clf.coef_[0, 1]
|
|
b3 = -sgd_clf.intercept_[0]/sgd_clf.coef_[0, 1]
|
|
|
|
# Transform the decision boundary lines back to the original scale
|
|
line1 = scaler.inverse_transform([[-10, -10 * w1 + b1], [10, 10 * w1 + b1]])
|
|
line2 = scaler.inverse_transform([[-10, -10 * w2 + b2], [10, 10 * w2 + b2]])
|
|
line3 = scaler.inverse_transform([[-10, -10 * w3 + b3], [10, 10 * w3 + b3]])
|
|
|
|
# Plot all three decision boundaries
|
|
plt.figure(figsize=(11, 4))
|
|
plt.plot(line1[:, 0], line1[:, 1], "k:", label="LinearSVC")
|
|
plt.plot(line2[:, 0], line2[:, 1], "b--", linewidth=2, label="SVC")
|
|
plt.plot(line3[:, 0], line3[:, 1], "r-", label="SGDClassifier")
|
|
plt.plot(X[:, 0][y==1], X[:, 1][y==1], "bs") # label="Iris-Versicolor"
|
|
plt.plot(X[:, 0][y==0], X[:, 1][y==0], "yo") # label="Iris-Setosa"
|
|
plt.xlabel("Petal length", fontsize=14)
|
|
plt.ylabel("Petal width", fontsize=14)
|
|
plt.legend(loc="upper center", fontsize=14)
|
|
plt.axis([0, 5.5, 0, 2])
|
|
|
|
plt.show()
|
|
|
|
|
|
|
|
!ec
|
|
|
|
|
|
|
|
|
|
!split
|
|
===== What is a hyperplane? =====
|
|
|
|
The aim of the SVM algorithm is to find a hyperplane in a
|
|
$p$-dimensional space, where $p$ is the number of features that
|
|
distinctly classifies the data points.
|
|
|
|
In a $p$-dimensional space, a hyperplane is what we call an affine subspace of dimension of $p-1$.
|
|
As an example, in two dimension, a hyperplane is simply as straight line while in three dimensions it is
|
|
a two-dimensional subspace, or stated simply, a plane.
|
|
|
|
In two dimensions, with the variables $x_1$ and $x_2$, the hyperplane is defined as
|
|
!bt
|
|
\[
|
|
b+w_1x_1+w_2x_2=0,
|
|
\]
|
|
!et
|
|
|
|
where $b$ is the intercept and $w_1$ and $w_2$ define the elements of a vector orthogonal to the line
|
|
$b+w_1x_1+w_2x_2=0$.
|
|
In two dimensions we define the vectors $\bm{x} =[x1,x2]$ and $\bm{w}=[w1,w2]$.
|
|
We can then rewrite the above equation as
|
|
|
|
!bt
|
|
\[
|
|
\bm{x}^T\bm{w}+b=0.
|
|
\]
|
|
!et
|
|
|
|
!split
|
|
===== A $p$-dimensional space of features =====
|
|
|
|
We limit ourselves to two classes of outputs $y_i$ and assign these classes the values $y_i = \pm 1$.
|
|
In a $p$-dimensional space of say $p$ features we have a hyperplane defines as
|
|
!bt
|
|
\[
|
|
b+wx_1+w_2x_2+\dots +w_px_p=0.
|
|
\]
|
|
!et
|
|
If we define a
|
|
matrix $\bm{X}=\left[\bm{x}_1,\bm{x}_2,\dots, \bm{x}_p\right]$
|
|
of dimension $n\times p$, where $n$ represents the observations for each feature and each vector $x_i$ is a column vector of the matrix $\bm{X}$,
|
|
!bt
|
|
\[
|
|
\bm{x}_i = \begin{bmatrix} x_{i1} \\ x_{i2} \\ \dots \\ \dots \\ x_{ip} \end{bmatrix}.
|
|
\]
|
|
!et
|
|
If the above condition is not met for a given vector $\bm{x}_i$ we have
|
|
!bt
|
|
\[
|
|
b+w_1x_{i1}+w_2x_{i2}+\dots +w_px_{ip} >0,
|
|
\]
|
|
!et
|
|
if our output $y_i=1$.
|
|
In this case we say that $\bm{x}_i$ lies on one of the sides of the hyperplane and if
|
|
!bt
|
|
\[
|
|
b+w_1x_{i1}+w_2x_{i2}+\dots +w_px_{ip} < 0,
|
|
\]
|
|
!et
|
|
for the class of observations $y_i=-1$,
|
|
then $\bm{x}_i$ lies on the other side.
|
|
|
|
Equivalently, for the two classes of observations we have
|
|
!bt
|
|
\[
|
|
y_i\left(b+w_1x_{i1}+w_2x_{i2}+\dots +w_px_{ip}\right) > 0.
|
|
\]
|
|
!et
|
|
|
|
When we try to separate hyperplanes, if it exists, we can use it to construct a natural classifier: a test observation is assigned a given class depending on which side of the hyperplane it is located.
|
|
|
|
!split
|
|
===== The two-dimensional case =====
|
|
|
|
Let us try to develop our intuition about SVMs by limiting ourselves to a two-dimensional
|
|
plane. To separate the two classes of data points, there are many
|
|
possible lines (hyperplanes if you prefer a more strict naming)
|
|
that could be chosen. Our objective is to find a
|
|
plane that has the maximum margin, i.e the maximum distance between
|
|
data points of both classes. Maximizing the margin distance provides
|
|
some reinforcement so that future data points can be classified with
|
|
more confidence.
|
|
|
|
What a linear classifier attempts to accomplish is to split the
|
|
feature space into two half spaces by placing a hyperplane between the
|
|
data points. This hyperplane will be our decision boundary. All
|
|
points on one side of the plane will belong to class one and all points
|
|
on the other side of the plane will belong to the second class two.
|
|
|
|
Unfortunately there are many ways in which we can place a hyperplane
|
|
to divide the data. Below is an example of two candidate hyperplanes
|
|
for our data sample.
|
|
|
|
!split
|
|
===== Getting into the details =====
|
|
|
|
Let us define the function
|
|
!bt
|
|
\[
|
|
f(x) = \bm{w}^T\bm{x}+b = 0,
|
|
\]
|
|
!et
|
|
as the function that determines the line $L$ that separates two classes (our two features), see the figure here.
|
|
|
|
|
|
Any point defined by $\bm{x}_i$ and $\bm{x}_2$ on the line $L$ will satisfy $\bm{w}^T(\bm{x}_1-\bm{x}_2)=0$.
|
|
|
|
The signed distance $\delta$ from any point defined by a vector $\bm{x}$ and a point $\bm{x}_0$ on the line $L$ is then
|
|
!bt
|
|
\[
|
|
\delta = \frac{1}{\vert\vert \bm{w}\vert\vert}(\bm{w}^T\bm{x}+b).
|
|
\]
|
|
!et
|
|
|
|
!split
|
|
===== First attempt at a minimization approach =====
|
|
|
|
How do we find the parameter $b$ and the vector $\bm{w}$? What we could
|
|
do is to define a cost function which now contains the set of all
|
|
misclassified points $M$ and attempt to minimize this function
|
|
|
|
!bt
|
|
\[
|
|
C(\bm{w},b) = -\sum_{i\in M} y_i(\bm{w}^T\bm{x}_i+b).
|
|
\]
|
|
!et
|
|
|
|
We could now for example define all values $y_i =1$ as misclassified in case we have $\bm{w}^T\bm{x}_i+b < 0$ and the opposite if we have $y_i=-1$. Taking the derivatives gives us
|
|
!bt
|
|
\[
|
|
\frac{\partial C}{\partial b} = -\sum_{i\in M} y_i,
|
|
\]
|
|
!et
|
|
and
|
|
!bt
|
|
\[
|
|
\frac{\partial C}{\partial \bm{w}} = -\sum_{i\in M} y_ix_i.
|
|
\]
|
|
!et
|
|
|
|
!split
|
|
===== Solving the equations =====
|
|
|
|
We can now use the Newton-Raphson method or different variants of the gradient descent family (from plain gradient descent to various stochastic gradient descent approaches) to solve the equations
|
|
!bt
|
|
\[
|
|
b \leftarrow b +\eta \frac{\partial C}{\partial b},
|
|
\]
|
|
!et
|
|
and
|
|
!bt
|
|
\[
|
|
\bm{w} \leftarrow \bm{w} +\eta \frac{\partial C}{\partial \bm{w}},
|
|
\]
|
|
!et
|
|
where $\eta$ is our by now well-known learning rate.
|
|
|
|
|
|
!split
|
|
===== Code Example =====
|
|
|
|
The equations we discussed above can be coded rather easily (the
|
|
framework is similar to what we developed for logistic
|
|
regression). We are going to set up a simple case with two classes only and we want to find a line which separates them the best possible way.
|
|
!bc pycod
|
|
|
|
!ec
|
|
|
|
!split
|
|
===== Problems with the Simpler Approach =====
|
|
|
|
|
|
There are however problems with this approach, although it looks
|
|
pretty straightforward to implement. When running the above code, we see that we can easily end up with many diffeent lines which separate the two classes.
|
|
|
|
|
|
For small
|
|
gaps between the entries, we may also end up needing many iterations
|
|
before the solutions converge and if the data cannot be separated
|
|
properly into two distinct classes, we may not experience a converge
|
|
at all.
|
|
|
|
!split
|
|
===== A better approach =====
|
|
|
|
A better approach is rather to try to define a large margin between
|
|
the two classes (if they are well separated from the beginning).
|
|
|
|
Thus, we wish to find a margin $M$ with $\bm{w}$ normalized to
|
|
$\vert\vert \bm{w}\vert\vert =1$ subject to the condition
|
|
|
|
!bt
|
|
\[
|
|
y_i(\bm{w}^T\bm{x}_i+b) \geq M \hspace{0.1cm}\forall i=1,2,\dots, p.
|
|
\]
|
|
!et
|
|
All points are thus at a signed distance from the decision boundary defined by the line $L$. The parameters $b$ and $w_1$ and $w_2$ define this line.
|
|
|
|
We seek thus the largest value $M$ defined by
|
|
!bt
|
|
\[
|
|
\frac{1}{\vert \vert \bm{w}\vert\vert}y_i(\bm{w}^T\bm{x}_i+b) \geq M \hspace{0.1cm}\forall i=1,2,\dots, n,
|
|
\]
|
|
!et
|
|
or just
|
|
!bt
|
|
\[
|
|
y_i(\bm{w}^T\bm{x}_i+b) \geq M\vert \vert \bm{w}\vert\vert \hspace{0.1cm}\forall i.
|
|
\]
|
|
!et
|
|
If we scale the equation so that $\vert \vert \bm{w}\vert\vert = 1/M$, we have to find the minimum of
|
|
$\bm{w}^T\bm{w}=\vert \vert \bm{w}\vert\vert$ (the norm) subject to the condition
|
|
!bt
|
|
\[
|
|
y_i(\bm{w}^T\bm{x}_i+b) \geq 1 \hspace{0.1cm}\forall i.
|
|
\]
|
|
!et
|
|
|
|
We have thus defined our margin as the invers of the norm of
|
|
$\bm{w}$. We want to minimize the norm in order to have a as large as
|
|
possible margin $M$. Before we proceed, we need to remind ourselves
|
|
about Lagrangian multipliers.
|
|
|
|
!split
|
|
===== A quick Reminder on Lagrangian Multipliers =====
|
|
|
|
Consider a function of three independent variables $f(x,y,z)$ . For the function $f$ to be an
|
|
extreme we have
|
|
!bt
|
|
\[
|
|
df=0.
|
|
\]
|
|
!et
|
|
A necessary and sufficient condition is
|
|
!bt
|
|
\[
|
|
\frac{\partial f}{\partial x} =\frac{\partial f}{\partial y}=\frac{\partial f}{\partial z}=0,
|
|
\]
|
|
!et
|
|
due to
|
|
!bt
|
|
\[
|
|
df = \frac{\partial f}{\partial x}dx+\frac{\partial f}{\partial y}dy+\frac{\partial f}{\partial z}dz.
|
|
\]
|
|
!et
|
|
In many problems the variables $x,y,z$ are often subject to constraints (such as those above for the margin)
|
|
so that they are no longer all independent. It is possible at least in principle to use each
|
|
constraint to eliminate one variable
|
|
and to proceed with a new and smaller set of independent varables.
|
|
|
|
The use of so-called Lagrangian multipliers is an alternative technique when the elimination
|
|
of variables is incovenient or undesirable. Assume that we have an equation of constraint on
|
|
the variables $x,y,z$
|
|
!bt
|
|
\[
|
|
\phi(x,y,z) = 0,
|
|
\]
|
|
!et
|
|
resulting in
|
|
!bt
|
|
\[
|
|
d\phi = \frac{\partial \phi}{\partial x}dx+\frac{\partial \phi}{\partial y}dy+\frac{\partial \phi}{\partial z}dz =0.
|
|
\]
|
|
!et
|
|
Now we cannot set anymore
|
|
!bt
|
|
\[
|
|
\frac{\partial f}{\partial x} =\frac{\partial f}{\partial y}=\frac{\partial f}{\partial z}=0,
|
|
\]
|
|
!et
|
|
if $df=0$ is wanted
|
|
because there are now only two independent variables! Assume $x$ and $y$ are the independent
|
|
variables.
|
|
Then $dz$ is no longer arbitrary.
|
|
|
|
!split
|
|
===== Adding the Multiplier =====
|
|
|
|
However, we can add to
|
|
!bt
|
|
\[
|
|
df = \frac{\partial f}{\partial x}dx+\frac{\partial f}{\partial y}dy+\frac{\partial f}{\partial z}dz,
|
|
\]
|
|
!et
|
|
a multiplum of $d\phi$, viz. $\lambda d\phi$, resulting in
|
|
!bt
|
|
\[
|
|
df+\lambda d\phi = (\frac{\partial f}{\partial z}+\lambda
|
|
\frac{\partial \phi}{\partial x})dx+(\frac{\partial f}{\partial y}+\lambda\frac{\partial \phi}{\partial y})dy+
|
|
(\frac{\partial f}{\partial z}+\lambda\frac{\partial \phi}{\partial z})dz =0.
|
|
\]
|
|
!et
|
|
Our multiplier is chosen so that
|
|
!bt
|
|
\[
|
|
\frac{\partial f}{\partial z}+\lambda\frac{\partial \phi}{\partial z} =0.
|
|
\]
|
|
!et
|
|
|
|
We need to remember that we took $dx$ and $dy$ to be arbitrary and thus we must have
|
|
!bt
|
|
\[
|
|
\frac{\partial f}{\partial x}+\lambda\frac{\partial \phi}{\partial x} =0,
|
|
\]
|
|
!et
|
|
and
|
|
!bt
|
|
\[
|
|
\frac{\partial f}{\partial y}+\lambda\frac{\partial \phi}{\partial y} =0.
|
|
\]
|
|
!et
|
|
When all these equations are satisfied, $df=0$. We have four unknowns, $x,y,z$ and
|
|
$\lambda$. Actually we want only $x,y,z$, $\lambda$ needs not to be determined,
|
|
it is therefore often called
|
|
Lagrange's undetermined multiplier.
|
|
If we have a set of constraints $\phi_k$ we have the equations
|
|
!bt
|
|
\[
|
|
\frac{\partial f}{\partial x_i}+\sum_k\lambda_k\frac{\partial \phi_k}{\partial x_i} =0.
|
|
\]
|
|
!et
|
|
|
|
!split
|
|
===== Setting up the Problem =====
|
|
In order to solve the above problem, we define the following Lagrangian function to be minimized
|
|
!bt
|
|
\[
|
|
{\cal L}(\lambda,b,\bm{w})=\frac{1}{2}\bm{w}^T\bm{w}-\sum_{i=1}^n\lambda_i\left[y_i(\bm{w}^T\bm{x}_i+b)-1\right],
|
|
\]
|
|
!et
|
|
where $\lambda_i$ is a so-called Lagrange multiplier subject to the condition $\lambda_i \geq 0$.
|
|
|
|
Taking the derivatives with respect to $b$ and $\bm{w}$ we obtain
|
|
!bt
|
|
\[
|
|
\frac{\partial {\cal L}}{\partial b} = -\sum_{i} \lambda_iy_i=0,
|
|
\]
|
|
!et
|
|
and
|
|
!bt
|
|
\[
|
|
\frac{\partial {\cal L}}{\partial \bm{w}} = 0 = \bm{w}-\sum_{i} \lambda_iy_i\bm{x}_i.
|
|
\]
|
|
!et
|
|
Inserting these constraints into the equation for ${\cal L}$ we obtain
|
|
!bt
|
|
\[
|
|
{\cal L}=\sum_i\lambda_i-\frac{1}{2}\sum_{ij}^n\lambda_i\lambda_jy_iy_j\bm{x}_i^T\bm{x}_j,
|
|
\]
|
|
!et
|
|
subject to the constraints $\lambda_i\geq 0$ and $\sum_i\lambda_iy_i=0$.
|
|
We must in addition satisfy the "Karush-Kuhn-Tucker":"https://en.wikipedia.org/wiki/Karush%E2%80%93Kuhn%E2%80%93Tucker_conditions" (KKT) condition
|
|
!bt
|
|
\[
|
|
\lambda_i\left[y_i(\bm{w}^T\bm{x}_i+b) -1\right] \hspace{0.1cm}\forall i.
|
|
\]
|
|
!et
|
|
o If $\lambda_i > 0$, then $y_i(\bm{w}^T\bm{x}_i+b)=1$ and we say that $x_i$ is on the boundary.
|
|
o If $y_i(\bm{w}^T\bm{x}_i+b)> 1$, we say $x_i$ is not on the boundary and we set $\lambda_i=0$.
|
|
When $\lambda_i > 0$, the vectors $\bm{x}_i$ are called support vectors. They are the vectors closest to the line (or hyperplane) and define the margin $M$.
|
|
|
|
!split
|
|
===== The problem to solve =====
|
|
|
|
We can rewrite
|
|
!bt
|
|
\[
|
|
{\cal L}=\sum_i\lambda_i-\frac{1}{2}\sum_{ij}^n\lambda_i\lambda_jy_iy_j\bm{x}_i^T\bm{x}_j,
|
|
\]
|
|
!et
|
|
and its constraints in terms of a matrix-vector problem where we minimize w.r.t. $\lambda$ the following problem
|
|
!bt
|
|
\[
|
|
\frac{1}{2} \bm{\lambda}^T\begin{bmatrix} y_1y_1\bm{x}_1^T\bm{x}_1 & y_1y_2\bm{x}_1^T\bm{x}_2 & \dots & \dots & y_1y_n\bm{x}_1^T\bm{x}_n \\
|
|
y_2y_1\bm{x}_2^T\bm{x}_1 & y_2y_2\bm{x}_2^T\bm{x}_2 & \dots & \dots & y_1y_n\bm{x}_2^T\bm{x}_n \\
|
|
\dots & \dots & \dots & \dots & \dots \\
|
|
\dots & \dots & \dots & \dots & \dots \\
|
|
y_ny_1\bm{x}_n^T\bm{x}_1 & y_ny_2\bm{x}_n^T\bm{x}_2 & \dots & \dots & y_ny_n\bm{x}_n^T\bm{x}_n \\
|
|
\end{bmatrix}\bm{\lambda}-\mathbb{1}\bm{\lambda},
|
|
\]
|
|
!et
|
|
subject to $\bm{y}^T\bm{\lambda}=0$. Here we defined the vectors $\bm{\lambda} =[\lambda_1,\lambda_2,\dots,\lambda_n]$ and
|
|
$\bm{y}=[y_1,y_2,\dots,y_n]$.
|
|
|
|
|
|
!split
|
|
===== The last steps =====
|
|
|
|
Solving the above problem, yields the values of $\lambda_i$.
|
|
To find the coefficients of your hyperplane we need simply to compute
|
|
!bt
|
|
\[
|
|
\bm{w}=\sum_{i} \lambda_iy_i\bm{x}_i.
|
|
\]
|
|
!et
|
|
With our vector $\bm{w}$ we can in turn find the value of the intercept $b$ (here in two dimensions) via
|
|
!bt
|
|
\[
|
|
y_i(\bm{w}^T\bm{x}_i+b)=1,
|
|
\]
|
|
!et
|
|
resulting in
|
|
!bt
|
|
\[
|
|
b = \frac{1}{y_i}-\bm{w}^T\bm{x}_i,
|
|
\]
|
|
!et
|
|
or if we write it out in terms of the support vectors only, with $N_s$ being their number, we have
|
|
!bt
|
|
\[
|
|
b = \frac{1}{N_s}\sum_{j\in N_s}\left(y_j-\sum_{i=1}^n\lambda_iy_i\bm{x}_i^T\bm{x}_j\right).
|
|
\]
|
|
!et
|
|
With our hyperplane coefficients we can use our classifier to assign any observation by simply using
|
|
!bt
|
|
\[
|
|
y_i = \mathrm{sign}(\bm{w}^T\bm{x}_i+b).
|
|
\]
|
|
!et
|
|
Below we discuss how to find the optimal values of $\lambda_i$. Before we proceed however, we discuss now the so-called soft classifier.
|
|
|
|
!split
|
|
===== A soft classifier =====
|
|
|
|
Till now, the margin is strictly defined by the support vectors. This defines what is called a hard classifier, that is the margins are well defined.
|
|
|
|
Suppose now that classes overlap in feature space, as shown in the
|
|
figure here. One way to deal with this problem before we define the
|
|
so-called _kernel approach_, is to allow a kind of slack in the sense
|
|
that we allow some points to be on the wrong side of the margin.
|
|
|
|
We introduce thus the so-called _slack_ variables $\bm{\xi} =[\xi_1,x_2,\dots,x_n]$ and
|
|
modify our previous equation
|
|
!bt
|
|
\[
|
|
y_i(\bm{w}^T\bm{x}_i+b)=1,
|
|
\]
|
|
!et
|
|
to
|
|
!bt
|
|
\[
|
|
y_i(\bm{w}^T\bm{x}_i+b)=1-\xi_i,
|
|
\]
|
|
!et
|
|
with the requirement $\xi_i\geq 0$. The total violation is now $\sum_i\xi$.
|
|
The value $\xi_i$ in the constraint the last constraint corresponds to the amount by which the prediction
|
|
$y_i(\bm{w}^T\bm{x}_i+b)=1$ is on the wrong side of its margin. Hence by bounding the sum $\sum_i \xi_i$,
|
|
we bound the total amount by which predictions fall on the wrong side of their margins.
|
|
|
|
Misclassifications occur when $\xi_i > 1$. Thus bounding the total sum by some value $C$ bounds in turn the total number of
|
|
misclassifications.
|
|
|
|
!split
|
|
===== Soft optmization problem =====
|
|
|
|
|
|
This has in turn the consequences that we change our optmization problem to finding the minimum of
|
|
!bt
|
|
\[
|
|
{\cal L}=\frac{1}{2}\bm{w}^T\bm{w}-\sum_{i=1}^n\lambda_i\left[y_i(\bm{w}^T\bm{x}_i+b)-(1-\xi_)\right]+C\sum_{i=1}^n\xi_i-\sum_{i=1}^n\gamma_i\xi_i,
|
|
\]
|
|
!et
|
|
subject to
|
|
!bt
|
|
\[
|
|
y_i(\bm{w}^T\bm{x}_i+b)=1-\xi_i \hspace{0.1cm}\forall i,
|
|
\]
|
|
!et
|
|
with the requirement $\xi_i\geq 0$.
|
|
|
|
Taking the derivatives with respect to $b$ and $\bm{w}$ we obtain
|
|
!bt
|
|
\[
|
|
\frac{\partial {\cal L}}{\partial b} = -\sum_{i} \lambda_iy_i=0,
|
|
\]
|
|
!et
|
|
and
|
|
!bt
|
|
\[
|
|
\frac{\partial {\cal L}}{\partial \bm{w}} = 0 = \bm{w}-\sum_{i} \lambda_iy_i\bm{x}_i,
|
|
\]
|
|
!et
|
|
and
|
|
!bt
|
|
\[
|
|
\lambda_i = C-\gamma_i \hspace{0.1cm}\forall i.
|
|
\]
|
|
!et
|
|
Inserting these constraints into the equation for ${\cal L}$ we obtain the same equation as before
|
|
!bt
|
|
\[
|
|
{\cal L}=\sum_i\lambda_i-\frac{1}{2}\sum_{ij}^n\lambda_i\lambda_jy_iy_j\bm{x}_i^T\bm{x}_j,
|
|
\]
|
|
!et
|
|
but now subject to the constraints $\lambda_i\geq 0$, $\sum_i\lambda_iy_i=0$ and $0\leq\lambda_i \leq C$.
|
|
We must in addition satisfy the Karush-Kuhn-Tucker condition which now reads
|
|
!bt
|
|
\[
|
|
\lambda_i\left[y_i(\bm{w}^T\bm{x}_i+b) -(1-\xi_)\right]=0 \hspace{0.1cm}\forall i,
|
|
\]
|
|
!et
|
|
!bt
|
|
\[
|
|
\gamma_i\xi_i = 0,
|
|
\]
|
|
!et
|
|
and
|
|
!bt
|
|
\[
|
|
y_i(\bm{w}^T\bm{x}_i+b) -(1-\xi_) \geq 0 \hspace{0.1cm}\forall i.
|
|
\]
|
|
!et
|
|
|