Chapter 9 Support Vector Machine¶

82¶

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import scipy
from scipy import stats
from numpy.random import randn  # Gaussian random number
import cvxopt
from cvxopt import matrix
In [ ]:
a = randn(1)
b = randn(1)
n = 100
X = randn(n, 2)
y = np.sign(a*X[:, 0] + b*X[:, 1] + 0.1*randn(n))
y = y.reshape(-1, 1)  # The shape should be explicit
In [ ]:
def svm_1(X, y, C):
    eps = 0.0001
    n = X.shape[0]
    P = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            P[i, j] = np.dot(X[i, :], X[j, :]) * y[i] * y[j]
    # Use the matrix function in the package to spacify them
    P = matrix(P + np.eye(n)*eps)
    A = matrix(- y.T.astype(np.float64))
    b = matrix(np.array([0]).astype(np.float64))
    h = matrix(np.array([C]*n + [0]*n).reshape(-1, 1).astype(np.float64))
    G = matrix(np.concatenate([np.diag(np.ones(n)), np.diag(-np.ones(n))]))
    q = matrix(np.array([-1]*n).astype(np.float64))
    res = cvxopt.solvers.qp(P, q, A=A, b=b, G=G, h=h)    # Execute the solver
    alpha = np.array(res["x"])  # x is the alpha in the text
    beta = ((alpha*y).T@X).reshape(2, 1)
    index = (eps < alpha[:, 0]) & (alpha[:, 0] < C - eps)
    beta_0 = np.mean(y[index] - X[index, :]@beta)
    return {"beta":beta, "beta_0":beta_0}
In [ ]:
def f(x):
    return - res["beta_0"]/res["beta"][1] - x*res["beta"][0]/res["beta"][1]
In [ ]:
a = randn(1)
b = randn(1)
n = 100
X = randn(n, 2)
y = np.sign(a*X[:, 0] + b*X[:, 1] + 0.1*randn(n))
y = y.reshape(-1, 1)  # The shape should be explicit.
for i in range(n):
    if y[i] == 1:
        plt.scatter(X[i, 0], X[i, 1], c="red")
    else :
        plt.scatter(X[i, 0], X[i, 1], c="blue")
res = svm_1(X, y, C=10)
x_seq = np.arange(-3, 3, 0.5)
plt.plot(x_seq, f(x_seq))
<ipython-input-14-9ae01250edad>:7: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
  P[i, j] = np.dot(X[i, :], X[j, :]) * y[i] * y[j]
     pcost       dcost       gap    pres   dres
 0: -1.6325e+02 -8.9101e+03  3e+04  1e+00  1e-14
 1: -3.0366e+01 -3.1994e+03  5e+03  1e-01  1e-14
 2:  1.4212e+01 -4.3868e+02  6e+02  1e-02  3e-14
 3: -6.3227e+01 -1.6285e+02  1e+02  1e-03  7e-15
 4: -8.2897e+01 -1.3421e+02  5e+01  5e-04  5e-15
 5: -9.4134e+01 -1.2559e+02  3e+01  2e-04  5e-15
 6: -1.0001e+02 -1.1851e+02  2e+01  9e-05  5e-15
 7: -1.0399e+02 -1.1207e+02  8e+00  3e-05  5e-15
 8: -1.0717e+02 -1.0812e+02  1e+00  2e-16  5e-15
 9: -1.0755e+02 -1.0763e+02  7e-02  2e-15  6e-15
10: -1.0758e+02 -1.0759e+02  2e-03  7e-15  5e-15
11: -1.0758e+02 -1.0758e+02  2e-05  2e-15  6e-15
Optimal solution found.
Out[ ]:
[<matplotlib.lines.Line2D at 0x7c2524c3e170>]
No description has been provided for this image

84¶

In [ ]:
def K_linear(x, y):
    return x.T@y
In [ ]:
# The polynomial kernel may be constructed as folows:
def K_poly(x, y):
    return (1 + x.T@y) ** 2

85¶

In [ ]:
# Execution
a = 3
b = -1
n = 200
X = randn(n, 2)
y = np.sign(a*X[:, 0] + b*X[:, 1]**2 + 0.3*randn(n))
y = y.reshape(-1, 1)
In [ ]:
def plot_kernel(K, line):  # Specify the line format using ``line"
    res = svm_2(X, y, 1, K)
    alpha = res["alpha"][:, 0]
    beta_0 = res["beta_0"]
    def f(u, v):
        S = beta_0
        for i in range(X.shape[0]):
            S = S + alpha[i]*y[i]*K(X[i, :], [u, v])
        return S[0]
    uu = np.arange(-2, 2, 0.1)
    vv = np.arange(-2, 2, 0.1)
    ww = []
    for v in vv:
        w = []
        for u in uu:
            w.append(f(u, v))
        ww.append(w)
    plt.contour(uu, vv, ww, levels=0, linestyles=line)
In [ ]:
def svm_2(X, y, C, K):
    eps = 0.0001
    n = X.shape[0]
    P = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            P[i, j] = K(X[i, :], X[j, :])*y[i]*y[j]
    # Use the matrix function in the package to spacify them
    P = matrix(P + np.eye(n)*eps)
    A = matrix(-y.T.astype(np.float64))
    b = matrix(np.array([0]).astype(np.float64))
    h = matrix(np.array([C]*n + [0]*n).reshape(-1, 1).astype(np.float64))
    G = matrix(np.concatenate([np.diag(np.ones(n)), np.diag(-np.ones(n))]))
    q = matrix(np.array([-1]*n).astype(np.float64))
    res = cvxopt.solvers.qp(P, q, A=A, b=b, G=G, h=h)
    alpha = np.array(res["x"])  # x is the alpha in the text
    beta = ((alpha*y).T@X).reshape(2, 1)
    index = (eps < alpha[:, 0]) & (alpha[:, 0] < C - eps)
    beta_0 = np.mean(y[index] - X[index, :]@beta)
    return {"alpha":alpha, "beta":beta, "beta_0":beta_0}
    uu = np.arange(-2, 2, 0.1)
    vv = np.arange(-2, 2, 0.1)
    ww = []
    for v in vv:
        w = []
        for u in uu:
            w.append(f(u, v))
        ww.append(w)
    plt.contour(uu, vv, ww, levels=0, linestyles=line)
In [ ]:
for i in range(n):
    if y[i] == 1:
        plt.scatter(X[i, 0], X[i, 1], c="red")
    else:
        plt.scatter(X[i, 0], X[i, 1], c="blue")
plot_kernel(K_poly, line="dashed")
plot_kernel(K_linear, line="solid")
<ipython-input-23-afc9250a7f7c>:7: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
  P[i, j] = K(X[i, :], X[j, :])*y[i]*y[j]
     pcost       dcost       gap    pres   dres
 0: -6.9533e+01 -4.9851e+02  3e+03  3e+00  2e-14
 1: -4.3954e+01 -3.2535e+02  6e+02  5e-01  1e-14
 2: -2.7322e+01 -1.3056e+02  2e+02  1e-01  7e-15
 3: -1.9699e+01 -5.0913e+01  5e+01  4e-02  4e-15
 4: -1.9472e+01 -2.7246e+01  1e+01  8e-03  4e-15
 5: -2.0169e+01 -2.2575e+01  3e+00  2e-03  4e-15
 6: -2.0477e+01 -2.1305e+01  1e+00  3e-04  4e-15
 7: -2.0647e+01 -2.0954e+01  3e-01  6e-05  4e-15
 8: -2.0766e+01 -2.0787e+01  2e-02  4e-06  4e-15
 9: -2.0775e+01 -2.0775e+01  4e-04  6e-08  5e-15
10: -2.0775e+01 -2.0775e+01  8e-06  9e-10  5e-15
Optimal solution found.
     pcost       dcost       gap    pres   dres
 0: -8.5617e+01 -4.6975e+02  2e+03  3e+00  6e-15
 1: -5.7101e+01 -3.0106e+02  4e+02  3e-01  3e-15
 2: -4.7762e+01 -8.7320e+01  5e+01  2e-02  7e-15
 3: -5.4901e+01 -6.4934e+01  1e+01  4e-03  3e-15
 4: -5.6952e+01 -6.1729e+01  5e+00  2e-03  3e-15
 5: -5.7962e+01 -6.0085e+01  2e+00  5e-04  4e-15
 6: -5.8414e+01 -5.9299e+01  9e-01  7e-05  3e-15
 7: -5.8722e+01 -5.8943e+01  2e-01  6e-06  3e-15
 8: -5.8819e+01 -5.8830e+01  1e-02  3e-07  3e-15
 9: -5.8824e+01 -5.8824e+01  3e-04  1e-09  3e-15
10: -5.8824e+01 -5.8824e+01  3e-06  1e-11  3e-15
Optimal solution found.
No description has been provided for this image

86¶

In [ ]:
import sklearn
from sklearn import svm
In [ ]:
x = randn(200, 2)
x[0:100, ] = x[0:100, ] + 2
x[100:150, ] = x[100:150, ] - 2
y = np.concatenate(([1 for i in range(150)], [2 for i in range(50)]))
train = np.random.choice(200, 100, replace=False)
test = list(set(range(200)) - set(train))
res_svm = svm.SVC(kernel="rbf", gamma=1, C=100)  # SVM without tuning
res_svm.fit(x[train, ], y[train])  # Execution
Out[ ]:
SVC(C=100, gamma=1)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
SVC(C=100, gamma=1)
In [ ]:
res_svm.predict(x[test, ])      # Prediction of the test data
Out[ ]:
array([1, 1, 1, 1, 1, 1, 1, 2, 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, 1, 1, 1, 2, 1, 1, 2, 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, 1,
       1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 1,
       1, 1, 2, 2, 2, 1, 2, 2, 1, 2, 2, 2])
In [ ]:
import mlxtend
from mlxtend.plotting import plot_decision_regions
In [ ]:
plot_decision_regions(x, y, clf=res_svm)
Out[ ]:
<Axes: >
No description has been provided for this image
In [ ]:
from sklearn.model_selection import GridSearchCV
In [ ]:
grid = {"C": [0.1, 1, 10, 100, 1000],  "gamma" : [0.5, 1, 2, 3, 4]}
tune = GridSearchCV(svm.SVC(), grid, cv=10)
tune.fit(x[train, ], y[train])
Out[ ]:
GridSearchCV(cv=10, estimator=SVC(),
             param_grid={'C': [0.1, 1, 10, 100, 1000],
                         'gamma': [0.5, 1, 2, 3, 4]})
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
GridSearchCV(cv=10, estimator=SVC(),
             param_grid={'C': [0.1, 1, 10, 100, 1000],
                         'gamma': [0.5, 1, 2, 3, 4]})
SVC()
SVC()

87¶

In [ ]:
from sklearn.datasets import load_iris
In [ ]:
iris = load_iris()
iris.target_names
x = iris.data
y = iris.target
train = np.random.choice(150, 120, replace=False)
test = np.ones(150, dtype=bool)
test[train] = False
iris_svm = svm.SVC(kernel="rbf", gamma=1, C=10)
iris_svm.fit(x[train, ], y[train])
Out[ ]:
SVC(C=10, gamma=1)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
SVC(C=10, gamma=1)
In [ ]:
# the function table_count in Chapter 3
def table_count(m, u, v):
    n = u.shape[0]
    count = np.zeros([m, m])
    for i in range(n):
        count[int(u[i]), int(v[i])] += 1
    return (count)
In [ ]:
y_pre = iris_svm.predict(x[test, ])
table_count(3, y[test], y_pre)
Out[ ]:
array([[11.,  0.,  0.],
       [ 0., 12.,  0.],
       [ 0.,  1.,  6.]])