updating week47
This commit is contained in:
mhjensen
@@ -6,7 +6,7 @@ DATE: today
|
||||
===== Overview of week 47 =====
|
||||
|
||||
* _Thursday_: Support Vector Machines, classification and regression. "Video of Lecture":"https://www.uio.no/studier/emner/matnat/fys/FYS-STK4155/h20/forelesningsvideoer/LectureNovember19.mp4?vrtx=view-as-webpage"
|
||||
* _Friday_: Workshop on project 3 (first lecture), Support Vector Machines (second Lecture)
|
||||
* _Friday_: Workshop on project 3. "Video of Lecture":"https://www.uio.no/studier/emner/matnat/fys/FYS-STK4155/h20/forelesningsvideoer/LectureNovember20.mp4?vrtx=view-as-webpage"
|
||||
|
||||
|
||||
Geron's chapter 5. Chapter 12 (sections 12.1-12.3 are the most relevant ones) of Hastie et al contains also a good discussion.
|
||||
@@ -664,514 +664,3 @@ y_i(\bm{w}^T\bm{x}_i+b) -(1-\xi_) \geq 0 \hspace{0.1cm}\forall i.
|
||||
\]
|
||||
!et
|
||||
|
||||
!split
|
||||
===== Kernels and non-linearity =====
|
||||
|
||||
The cases we have studied till now, were all characterized by two classes
|
||||
with a close to linear separability. The classifiers we have described
|
||||
so far find linear boundaries in our input feature space. It is
|
||||
possible to make our procedure more flexible by exploring the feature
|
||||
space using other basis expansions such as higher-order polynomials,
|
||||
wavelets, splines etc.
|
||||
|
||||
If our feature space is not easy to separate, as shown in the figure
|
||||
here, we can achieve a better separation by introducing more complex
|
||||
basis functions. The ideal would be, as shown in the next figure, to, via a specific transformation to
|
||||
obtain a separation between the classes which is almost linear.
|
||||
|
||||
The change of basis, from $x\rightarrow z=\phi(x)$ leads to the same type of equations to be solved, except that
|
||||
we need to introduce for example a polynomial transformation to a two-dimensional training set.
|
||||
|
||||
!bc pycod
|
||||
import numpy as np
|
||||
import os
|
||||
|
||||
np.random.seed(42)
|
||||
|
||||
# To plot pretty figures
|
||||
import matplotlib
|
||||
import matplotlib.pyplot as plt
|
||||
plt.rcParams['axes.labelsize'] = 14
|
||||
plt.rcParams['xtick.labelsize'] = 12
|
||||
plt.rcParams['ytick.labelsize'] = 12
|
||||
|
||||
|
||||
from sklearn.svm import SVC
|
||||
from sklearn import datasets
|
||||
|
||||
|
||||
|
||||
X1D = np.linspace(-4, 4, 9).reshape(-1, 1)
|
||||
X2D = np.c_[X1D, X1D**2]
|
||||
y = np.array([0, 0, 1, 1, 1, 1, 1, 0, 0])
|
||||
|
||||
plt.figure(figsize=(11, 4))
|
||||
|
||||
plt.subplot(121)
|
||||
plt.grid(True, which='both')
|
||||
plt.axhline(y=0, color='k')
|
||||
plt.plot(X1D[:, 0][y==0], np.zeros(4), "bs")
|
||||
plt.plot(X1D[:, 0][y==1], np.zeros(5), "g^")
|
||||
plt.gca().get_yaxis().set_ticks([])
|
||||
plt.xlabel(r"$x_1$", fontsize=20)
|
||||
plt.axis([-4.5, 4.5, -0.2, 0.2])
|
||||
|
||||
plt.subplot(122)
|
||||
plt.grid(True, which='both')
|
||||
plt.axhline(y=0, color='k')
|
||||
plt.axvline(x=0, color='k')
|
||||
plt.plot(X2D[:, 0][y==0], X2D[:, 1][y==0], "bs")
|
||||
plt.plot(X2D[:, 0][y==1], X2D[:, 1][y==1], "g^")
|
||||
plt.xlabel(r"$x_1$", fontsize=20)
|
||||
plt.ylabel(r"$x_2$", fontsize=20, rotation=0)
|
||||
plt.gca().get_yaxis().set_ticks([0, 4, 8, 12, 16])
|
||||
plt.plot([-4.5, 4.5], [6.5, 6.5], "r--", linewidth=3)
|
||||
plt.axis([-4.5, 4.5, -1, 17])
|
||||
plt.subplots_adjust(right=1)
|
||||
plt.show()
|
||||
|
||||
!ec
|
||||
|
||||
|
||||
|
||||
!split
|
||||
===== The equations =====
|
||||
|
||||
Suppose we define a polynomial transformation of degree two only (we continue to live in a plane with $x_i$ and $y_i$ as variables)
|
||||
!bt
|
||||
\[
|
||||
z = \phi(x_i) =\left(x_i^2, y_i^2, \sqrt{2}x_iy_i\right).
|
||||
\]
|
||||
!et
|
||||
|
||||
With our new basis, the equations we solved earlier are basically the same, that is we have now (without the slack option for simplicity)
|
||||
!bt
|
||||
\[
|
||||
{\cal L}=\sum_i\lambda_i-\frac{1}{2}\sum_{ij}^n\lambda_i\lambda_jy_iy_j\bm{z}_i^T\bm{z}_j,
|
||||
\]
|
||||
!et
|
||||
subject to the constraints $\lambda_i\geq 0$, $\sum_i\lambda_iy_i=0$, and for the support vectors
|
||||
!bt
|
||||
\[
|
||||
y_i(\bm{w}^T\bm{z}_i+b)= 1 \hspace{0.1cm}\forall i,
|
||||
\]
|
||||
!et
|
||||
from which we also find $b$.
|
||||
To compute $\bm{z}_i^T\bm{z}_j$ we define the kernel $K(\bm{x}_i,\bm{x}_j)$ as
|
||||
!bt
|
||||
\[
|
||||
K(\bm{x}_i,\bm{x}_j)=\bm{z}_i^T\bm{z}_j= \phi(\bm{x}_i)^T\phi(\bm{x}_j).
|
||||
\]
|
||||
!et
|
||||
For the above example, the kernel reads
|
||||
!bt
|
||||
\[
|
||||
K(\bm{x}_i,\bm{x}_j)=[x_i^2, y_i^2, \sqrt{2}x_iy_i]^T\begin{bmatrix} x_j^2 \\ y_j^2 \\ \sqrt{2}x_jy_j \end{bmatrix}=x_i^2x_j^2+2x_ix_jy_iy_j+y_i^2y_j^2.
|
||||
\]
|
||||
!et
|
||||
|
||||
We note that this is nothing but the dot product of the two original
|
||||
vectors $(\bm{x}_i^T\bm{x}_j)^2$. Instead of thus computing the
|
||||
product in the Lagrangian of $\bm{z}_i^T\bm{z}_j$ we simply compute
|
||||
the dot product $(\bm{x}_i^T\bm{x}_j)^2$.
|
||||
|
||||
|
||||
This leads to the so-called
|
||||
kernel trick and the result leads to the same as if we went through
|
||||
the trouble of performing the transformation
|
||||
$\phi(\bm{x}_i)^T\phi(\bm{x}_j)$ during the SVM calculations.
|
||||
|
||||
|
||||
!split
|
||||
===== The problem to solve =====
|
||||
Using our definition of the kernel We can rewrite again the Lagrangian
|
||||
!bt
|
||||
\[
|
||||
{\cal L}=\sum_i\lambda_i-\frac{1}{2}\sum_{ij}^n\lambda_i\lambda_jy_iy_j\bm{x}_i^T\bm{z}_j,
|
||||
\]
|
||||
!et
|
||||
subject to the constraints $\lambda_i\geq 0$, $\sum_i\lambda_iy_i=0$ in terms of a convex optimization problem
|
||||
!bt
|
||||
\[
|
||||
\frac{1}{2} \bm{\lambda}^T\begin{bmatrix} y_1y_1K(\bm{x}_1,\bm{x}_1) & y_1y_2K(\bm{x}_1,\bm{x}_2) & \dots & \dots & y_1y_nK(\bm{x}_1,\bm{x}_n) \\
|
||||
y_2y_1K(\bm{x}_2,\bm{x}_1) & y_2y_2(\bm{x}_2,\bm{x}_2) & \dots & \dots & y_1y_nK(\bm{x}_2,\bm{x}_n) \\
|
||||
\dots & \dots & \dots & \dots & \dots \\
|
||||
\dots & \dots & \dots & \dots & \dots \\
|
||||
y_ny_1K(\bm{x}_n,\bm{x}_1) & y_ny_2K(\bm{x}_n\bm{x}_2) & \dots & \dots & y_ny_nK(\bm{x}_n,\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]$.
|
||||
If we add the slack constants this leads to the additional constraint $0\leq \lambda_i \leq C$.
|
||||
|
||||
We can rewrite this (see the solutions below) in terms of a convex optimization problem of the type
|
||||
!bt
|
||||
\begin{align*}
|
||||
&\mathrm{min}_{\lambda}\hspace{0.2cm} \frac{1}{2}\bm{\lambda}^T\bm{P}\bm{\lambda}+\bm{q}^T\bm{\lambda},\\ \nonumber
|
||||
&\mathrm{subject\hspace{0.1cm}to} \hspace{0.2cm} \bm{G}\bm{\lambda} \preceq \bm{h} \hspace{0.2cm} \wedge \bm{A}\bm{\lambda}=f.
|
||||
\end{align*}
|
||||
!et
|
||||
Below we discuss how to solve these equations. Here we note that the matrix $\bm{P}$ has matrix elements $p_{ij}=y_iy_jK(\bm{x}_i,\bm{x}_j)$.
|
||||
Given a kernel $K$ and the targets $y_i$ this matrix is easy to set up. The constraint $\bm{y}^T\bm{\lambda}=0$ leads to $f=0$ and $\bm{A}=\bm{y}$. How to set up the matrix $\bm{G}$ is discussed later. Here note that the inequalities $0\leq \lambda_i \leq C$ can be split up into
|
||||
$0\leq \lambda_i$ and $\lambda_i \leq C$. These two inequalities define then the matrix $\bm{G}$ and the vector $\bm{h}$.
|
||||
|
||||
|
||||
!split
|
||||
===== Different kernels and Mercer's theorem =====
|
||||
|
||||
There are several popular kernels being used. These are
|
||||
o Linear: $K(\bm{x},\bm{y})=\bm{x}^T\bm{y}$,
|
||||
o Polynomial: $K(\bm{x},\bm{y})=(\bm{x}^T\bm{y}+\gamma)^d$,
|
||||
o Gaussian Radial Basis Function: $K(\bm{x},\bm{y})=\exp{\left(-\gamma\vert\vert\bm{x}-\bm{y}\vert\vert^2\right)}$,
|
||||
o Tanh: $K(\bm{x},\bm{y})=\tanh{(\bm{x}^T\bm{y}+\gamma)}$,
|
||||
and many other ones.
|
||||
|
||||
An important theorem for us is "Mercer's
|
||||
theorem":"https://en.wikipedia.org/wiki/Mercer%27s_theorem". The
|
||||
theorem states that if a kernel function $K$ is symmetric, continuous
|
||||
and leads to a positive semi-definite matrix $\bm{P}$ then there
|
||||
exists a function $\phi$ that maps $\bm{x}_i$ and $\bm{x}_j$ into
|
||||
another space (possibly with much higher dimensions) such that
|
||||
|
||||
!bt
|
||||
\[
|
||||
K(\bm{x}_i,\bm{x}_j)=\phi(\bm{x}_i)^T\phi(\bm{x}_j).
|
||||
\]
|
||||
!et
|
||||
|
||||
So you can use $K$ as a kernel since you know $\phi$ exists, even if
|
||||
you don’t know what $\phi$ is.
|
||||
|
||||
Note that some frequently used kernels (such as the Sigmoid kernel)
|
||||
don’t respect all of Mercer’s conditions, yet they generally work well
|
||||
in practice.
|
||||
|
||||
|
||||
!split
|
||||
===== The moons example =====
|
||||
!bc pycod
|
||||
from __future__ import division, print_function, unicode_literals
|
||||
|
||||
import numpy as np
|
||||
np.random.seed(42)
|
||||
|
||||
import matplotlib
|
||||
import matplotlib.pyplot as plt
|
||||
plt.rcParams['axes.labelsize'] = 14
|
||||
plt.rcParams['xtick.labelsize'] = 12
|
||||
plt.rcParams['ytick.labelsize'] = 12
|
||||
|
||||
|
||||
from sklearn.svm import SVC
|
||||
from sklearn import datasets
|
||||
|
||||
|
||||
|
||||
from sklearn.pipeline import Pipeline
|
||||
from sklearn.preprocessing import StandardScaler
|
||||
from sklearn.svm import LinearSVC
|
||||
|
||||
|
||||
from sklearn.datasets import make_moons
|
||||
X, y = make_moons(n_samples=100, noise=0.15, random_state=42)
|
||||
|
||||
def plot_dataset(X, y, axes):
|
||||
plt.plot(X[:, 0][y==0], X[:, 1][y==0], "bs")
|
||||
plt.plot(X[:, 0][y==1], X[:, 1][y==1], "g^")
|
||||
plt.axis(axes)
|
||||
plt.grid(True, which='both')
|
||||
plt.xlabel(r"$x_1$", fontsize=20)
|
||||
plt.ylabel(r"$x_2$", fontsize=20, rotation=0)
|
||||
|
||||
plot_dataset(X, y, [-1.5, 2.5, -1, 1.5])
|
||||
plt.show()
|
||||
|
||||
from sklearn.datasets import make_moons
|
||||
from sklearn.pipeline import Pipeline
|
||||
from sklearn.preprocessing import PolynomialFeatures
|
||||
|
||||
polynomial_svm_clf = Pipeline([
|
||||
("poly_features", PolynomialFeatures(degree=3)),
|
||||
("scaler", StandardScaler()),
|
||||
("svm_clf", LinearSVC(C=10, loss="hinge", random_state=42))
|
||||
])
|
||||
|
||||
polynomial_svm_clf.fit(X, y)
|
||||
|
||||
def plot_predictions(clf, axes):
|
||||
x0s = np.linspace(axes[0], axes[1], 100)
|
||||
x1s = np.linspace(axes[2], axes[3], 100)
|
||||
x0, x1 = np.meshgrid(x0s, x1s)
|
||||
X = np.c_[x0.ravel(), x1.ravel()]
|
||||
y_pred = clf.predict(X).reshape(x0.shape)
|
||||
y_decision = clf.decision_function(X).reshape(x0.shape)
|
||||
plt.contourf(x0, x1, y_pred, cmap=plt.cm.brg, alpha=0.2)
|
||||
plt.contourf(x0, x1, y_decision, cmap=plt.cm.brg, alpha=0.1)
|
||||
|
||||
plot_predictions(polynomial_svm_clf, [-1.5, 2.5, -1, 1.5])
|
||||
plot_dataset(X, y, [-1.5, 2.5, -1, 1.5])
|
||||
|
||||
plt.show()
|
||||
|
||||
|
||||
from sklearn.svm import SVC
|
||||
|
||||
poly_kernel_svm_clf = Pipeline([
|
||||
("scaler", StandardScaler()),
|
||||
("svm_clf", SVC(kernel="poly", degree=3, coef0=1, C=5))
|
||||
])
|
||||
poly_kernel_svm_clf.fit(X, y)
|
||||
|
||||
poly100_kernel_svm_clf = Pipeline([
|
||||
("scaler", StandardScaler()),
|
||||
("svm_clf", SVC(kernel="poly", degree=10, coef0=100, C=5))
|
||||
])
|
||||
poly100_kernel_svm_clf.fit(X, y)
|
||||
|
||||
plt.figure(figsize=(11, 4))
|
||||
|
||||
plt.subplot(121)
|
||||
plot_predictions(poly_kernel_svm_clf, [-1.5, 2.5, -1, 1.5])
|
||||
plot_dataset(X, y, [-1.5, 2.5, -1, 1.5])
|
||||
plt.title(r"$d=3, r=1, C=5$", fontsize=18)
|
||||
|
||||
plt.subplot(122)
|
||||
plot_predictions(poly100_kernel_svm_clf, [-1.5, 2.5, -1, 1.5])
|
||||
plot_dataset(X, y, [-1.5, 2.5, -1, 1.5])
|
||||
plt.title(r"$d=10, r=100, C=5$", fontsize=18)
|
||||
|
||||
plt.show()
|
||||
|
||||
def gaussian_rbf(x, landmark, gamma):
|
||||
return np.exp(-gamma * np.linalg.norm(x - landmark, axis=1)**2)
|
||||
|
||||
gamma = 0.3
|
||||
|
||||
x1s = np.linspace(-4.5, 4.5, 200).reshape(-1, 1)
|
||||
x2s = gaussian_rbf(x1s, -2, gamma)
|
||||
x3s = gaussian_rbf(x1s, 1, gamma)
|
||||
|
||||
XK = np.c_[gaussian_rbf(X1D, -2, gamma), gaussian_rbf(X1D, 1, gamma)]
|
||||
yk = np.array([0, 0, 1, 1, 1, 1, 1, 0, 0])
|
||||
|
||||
plt.figure(figsize=(11, 4))
|
||||
|
||||
plt.subplot(121)
|
||||
plt.grid(True, which='both')
|
||||
plt.axhline(y=0, color='k')
|
||||
plt.scatter(x=[-2, 1], y=[0, 0], s=150, alpha=0.5, c="red")
|
||||
plt.plot(X1D[:, 0][yk==0], np.zeros(4), "bs")
|
||||
plt.plot(X1D[:, 0][yk==1], np.zeros(5), "g^")
|
||||
plt.plot(x1s, x2s, "g--")
|
||||
plt.plot(x1s, x3s, "b:")
|
||||
plt.gca().get_yaxis().set_ticks([0, 0.25, 0.5, 0.75, 1])
|
||||
plt.xlabel(r"$x_1$", fontsize=20)
|
||||
plt.ylabel(r"Similarity", fontsize=14)
|
||||
plt.annotate(r'$\mathbf{x}$',
|
||||
xy=(X1D[3, 0], 0),
|
||||
xytext=(-0.5, 0.20),
|
||||
ha="center",
|
||||
arrowprops=dict(facecolor='black', shrink=0.1),
|
||||
fontsize=18,
|
||||
)
|
||||
plt.text(-2, 0.9, "$x_2$", ha="center", fontsize=20)
|
||||
plt.text(1, 0.9, "$x_3$", ha="center", fontsize=20)
|
||||
plt.axis([-4.5, 4.5, -0.1, 1.1])
|
||||
|
||||
plt.subplot(122)
|
||||
plt.grid(True, which='both')
|
||||
plt.axhline(y=0, color='k')
|
||||
plt.axvline(x=0, color='k')
|
||||
plt.plot(XK[:, 0][yk==0], XK[:, 1][yk==0], "bs")
|
||||
plt.plot(XK[:, 0][yk==1], XK[:, 1][yk==1], "g^")
|
||||
plt.xlabel(r"$x_2$", fontsize=20)
|
||||
plt.ylabel(r"$x_3$ ", fontsize=20, rotation=0)
|
||||
plt.annotate(r'$\phi\left(\mathbf{x}\right)$',
|
||||
xy=(XK[3, 0], XK[3, 1]),
|
||||
xytext=(0.65, 0.50),
|
||||
ha="center",
|
||||
arrowprops=dict(facecolor='black', shrink=0.1),
|
||||
fontsize=18,
|
||||
)
|
||||
plt.plot([-0.1, 1.1], [0.57, -0.1], "r--", linewidth=3)
|
||||
plt.axis([-0.1, 1.1, -0.1, 1.1])
|
||||
|
||||
plt.subplots_adjust(right=1)
|
||||
|
||||
plt.show()
|
||||
|
||||
|
||||
x1_example = X1D[3, 0]
|
||||
for landmark in (-2, 1):
|
||||
k = gaussian_rbf(np.array([[x1_example]]), np.array([[landmark]]), gamma)
|
||||
print("Phi({}, {}) = {}".format(x1_example, landmark, k))
|
||||
|
||||
rbf_kernel_svm_clf = Pipeline([
|
||||
("scaler", StandardScaler()),
|
||||
("svm_clf", SVC(kernel="rbf", gamma=5, C=0.001))
|
||||
])
|
||||
rbf_kernel_svm_clf.fit(X, y)
|
||||
|
||||
|
||||
from sklearn.svm import SVC
|
||||
|
||||
gamma1, gamma2 = 0.1, 5
|
||||
C1, C2 = 0.001, 1000
|
||||
hyperparams = (gamma1, C1), (gamma1, C2), (gamma2, C1), (gamma2, C2)
|
||||
|
||||
svm_clfs = []
|
||||
for gamma, C in hyperparams:
|
||||
rbf_kernel_svm_clf = Pipeline([
|
||||
("scaler", StandardScaler()),
|
||||
("svm_clf", SVC(kernel="rbf", gamma=gamma, C=C))
|
||||
])
|
||||
rbf_kernel_svm_clf.fit(X, y)
|
||||
svm_clfs.append(rbf_kernel_svm_clf)
|
||||
|
||||
plt.figure(figsize=(11, 7))
|
||||
|
||||
for i, svm_clf in enumerate(svm_clfs):
|
||||
plt.subplot(221 + i)
|
||||
plot_predictions(svm_clf, [-1.5, 2.5, -1, 1.5])
|
||||
plot_dataset(X, y, [-1.5, 2.5, -1, 1.5])
|
||||
gamma, C = hyperparams[i]
|
||||
plt.title(r"$\gamma = {}, C = {}$".format(gamma, C), fontsize=16)
|
||||
|
||||
plt.show()
|
||||
|
||||
!ec
|
||||
|
||||
|
||||
|
||||
!split
|
||||
===== Mathematical optimization of convex functions =====
|
||||
|
||||
A mathematical (quadratic) optimization problem, or just optimization problem, has the form
|
||||
!bt
|
||||
\begin{align*}
|
||||
&\mathrm{min}_{\lambda}\hspace{0.2cm} \frac{1}{2}\bm{\lambda}^T\bm{P}\bm{\lambda}+\bm{q}^T\bm{\lambda},\\ \nonumber
|
||||
&\mathrm{subject\hspace{0.1cm}to} \hspace{0.2cm} \bm{G}\bm{\lambda} \preceq \bm{h} \wedge \bm{A}\bm{\lambda}=f.
|
||||
\end{align*}
|
||||
!et
|
||||
subject to some constraints for say a selected set $i=1,2,\dots, n$.
|
||||
In our case we are optimizing with respect to the Lagrangian multipliers $\lambda_i$, and the
|
||||
vector $\bm{\lambda}=[\lambda_1, \lambda_2,\dots, \lambda_n]$ is the optimization variable we are dealing with.
|
||||
|
||||
In our case we are particularly interested in a class of optimization problems called convex optmization problems.
|
||||
In our discussion on gradient descent methods we discussed at length the definition of a convex function.
|
||||
|
||||
Convex optimization problems play a central role in applied mathematics and we recommend strongly "Boyd and Vandenberghe's text on the topics":"http://web.stanford.edu/~boyd/cvxbook/".
|
||||
|
||||
|
||||
|
||||
!split
|
||||
===== How do we solve these problems? =====
|
||||
|
||||
If we use Python as programming language and wish to venture beyond
|
||||
_scikit-learn_, _tensorflow_ and similar software which makes our
|
||||
lives so much easier, we need to dive into the wonderful world of
|
||||
quadratic programming. We can, if we wish, solve the minimization
|
||||
problem using say standard gradient methods or conjugate gradient
|
||||
methods. However, these methods tend to exhibit a rather slow
|
||||
converge. So, welcome to the promised land of quadratic programming.
|
||||
|
||||
The functions we need are contained in the quadratic programming package _CVXOPT_ and we need to import it together with _numpy_ as
|
||||
|
||||
!bc pycod
|
||||
import numpy
|
||||
import cvxopt
|
||||
!ec
|
||||
|
||||
This will make our life much easier. You don't need t write your own optimizer.
|
||||
|
||||
|
||||
!split
|
||||
===== A simple example =====
|
||||
|
||||
We remind ourselves about the general problem we want to solve
|
||||
!bt
|
||||
\begin{align*}
|
||||
&\mathrm{min}_{x}\hspace{0.2cm} \frac{1}{2}\bm{x}^T\bm{P}\bm{x}+\bm{q}^T\bm{x},\\ \nonumber
|
||||
&\mathrm{subject\hspace{0.1cm} to} \hspace{0.2cm} \bm{G}\bm{x} \preceq \bm{h} \wedge \bm{A}\bm{x}=f.
|
||||
\end{align*}
|
||||
!et
|
||||
|
||||
Let us show how to perform the optmization using a simple case. Assume we want to optimize the following problem
|
||||
!bt
|
||||
\begin{align*}
|
||||
&\mathrm{min}_{x}\hspace{0.2cm} \frac{1}{2}x^2+5x+3y \\ \nonumber
|
||||
&\mathrm{subject to} \\ \nonumber
|
||||
&x, y \geq 0 \\ \nonumber
|
||||
&x+3y \geq 15 \\ \nonumber
|
||||
&2x+5y \leq 100 \\ \nonumber
|
||||
&3x+4y \leq 80. \\ \nonumber
|
||||
\end{align*}
|
||||
!et
|
||||
The minimization problem can be rewritten in terms of vectors and matrices as (with $x$ and $y$ being the unknowns)
|
||||
!bt
|
||||
\[
|
||||
\frac{1}{2}\begin{bmatrix} x\\ y \end{bmatrix}^T \begin{bmatrix} 1 & 0\\ 0 & 0 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} + \begin{bmatrix}3\\ 4 \end{bmatrix}^T \begin{bmatrix}x \\ y \end{bmatrix}.
|
||||
\]
|
||||
!et
|
||||
Similarly, we can now set up the inequalities (we need to change $\geq$ to $\leq$ by multiplying with $-1$ on bot sides) as the following matrix-vector equation
|
||||
!bt
|
||||
\[
|
||||
\begin{bmatrix} -1 & 0 \\ 0 & -1 \\ -1 & -3 \\ 2 & 5 \\ 3 & 4\end{bmatrix}\begin{bmatrix} x \\ y\end{bmatrix} \preceq \begin{bmatrix}0 \\ 0\\ -15 \\ 100 \\ 80\end{bmatrix}.
|
||||
\]
|
||||
!et
|
||||
We have collapsed all the inequalities into a single matrix $\bm{G}$. We see also that our matrix
|
||||
!bt
|
||||
\[
|
||||
\bm{P} =\begin{bmatrix} 1 & 0\\ 0 & 0 \end{bmatrix}
|
||||
\]
|
||||
!et
|
||||
is clearly positive semi-definite (all eigenvalues larger or equal zero).
|
||||
Finally, the vector $\bm{h}$ is defined as
|
||||
!bt
|
||||
\[
|
||||
\bm{h} = \begin{bmatrix}0 \\ 0\\ -15 \\ 100 \\ 80\end{bmatrix}.
|
||||
\]
|
||||
!et
|
||||
|
||||
|
||||
Since we don't have any equalities the matrix $\bm{A}$ is set to zero
|
||||
The following code solves the equations for us
|
||||
!bc pycod
|
||||
# Import the necessary packages
|
||||
import numpy
|
||||
from cvxopt import matrix
|
||||
from cvxopt import solvers
|
||||
P = matrix(numpy.diag([1,0]), tc=’d’)
|
||||
q = matrix(numpy.array([3,4]), tc=’d’)
|
||||
G = matrix(numpy.array([[-1,0],[0,-1],[-1,-3],[2,5],[3,4]]), tc=’d’)
|
||||
h = matrix(numpy.array([0,0,-15,100,80]), tc=’d’)
|
||||
# Construct the QP, invoke solver
|
||||
sol = solvers.qp(P,q,G,h)
|
||||
# Extract optimal value and solution
|
||||
sol[’x’]
|
||||
sol[’primal objective’]
|
||||
!ec
|
||||
|
||||
!split
|
||||
===== Back to the more realistic cases =====
|
||||
|
||||
We are now ready to return to our setup of the optmization problem for a more realistic case. Introducing the _slack_ parameter $C$ we have
|
||||
!bt
|
||||
\[
|
||||
\frac{1}{2} \bm{\lambda}^T\begin{bmatrix} y_1y_1K(\bm{x}_1,\bm{x}_1) & y_1y_2K(\bm{x}_1,\bm{x}_2) & \dots & \dots & y_1y_nK(\bm{x}_1,\bm{x}_n) \\
|
||||
y_2y_1K(\bm{x}_2,\bm{x}_1) & y_2y_2K(\bm{x}_2,\bm{x}_2) & \dots & \dots & y_1y_nK(\bm{x}_2,\bm{x}_n) \\
|
||||
\dots & \dots & \dots & \dots & \dots \\
|
||||
\dots & \dots & \dots & \dots & \dots \\
|
||||
y_ny_1K(\bm{x}_n,\bm{x}_1) & y_ny_2K(\bm{x}_n\bm{x}_2) & \dots & \dots & y_ny_nK(\bm{x}_n,\bm{x}_n) \\
|
||||
\end{bmatrix}\bm{\lambda}-\mathbb{I}\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]$.
|
||||
With the slack constants this leads to the additional constraint $0\leq \lambda_i \leq C$.
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user