【SVM】soft margin SVM 代码以及详解

心已赠人 2023-01-21 14:22 152阅读 0赞

### 文章目录 ###

*  导包
 *  10-------------------------------------------------
 *  11-------------------------------------------------
 *  12-------------------------------------------------
 *  13-------------------------------------------------
 *  14-------------------------------------------------
 *  15-------------------------------------------------
 *  16-------------------------------------------------
 *  17-------------------------------------------------
 *  18-------------------------------------------------
 *  测试部分
 *   *  线性核函数
     *  多项式核函数
     *  高斯核函数
 *  整体代码

# 导包 #

import matplotlib.pyplot as plt
    import numpy as np
    import scipy as sp
    import cvxopt
    from cvxopt import matrix
    from cvxopt import solvers
    from numpy import linalg

# 10------------------------------------------------- #

Make the function `H = make_H_k(X,t, kernel_function, param)` to compute the matrix  H H H for a data matrix `X` and an associated target vector `t` given a kernel function `kernel_function` with parameter `param`.

def make_H_k(X, t, kernel_function, param):
        n_samples, n_feature = X.shape
        X_prime = np.zeros((n_samples, n_samples))
        for i in range(n_samples):
            for j in range(n_samples):
                X_prime[i, j] = kernel_function(X[i], X[j], param)
        H = np.outer(t, t) * X_prime
        return H

\- 在这一个步骤，kernel\_function可以为线性核函数（对应教案中的例子），多项式核函数（对应soft margin SVM），和高斯核函数。这里核函数的输入 x i , x j ∈ X x\_i, x\_j \\in\{X\} xi,xj∈X，

# 11------------------------------------------------- #

Make the function `P, q, G, h, A, b = setup_optimization_soft_k(X, t, C, kernel_function, param)` to compute the necessary matrices to solve the QP problem associated to the soft margin kernelized classification task using the library [cvxopt][]. The function receives in input a data matrix `X`, an associated target vector `t`, a regularization parameter `C` and a kernel function `kernel_function` with parameter `param`.

def setup_optimization_soft_k(X, t, C, kernel_function, param):
        m, n = X.shape
        t = t.reshape(-1, 1).astype(float)
        H = make_H_k(X, t, kernel_function, param)
        P = matrix(H)
        q = matrix(-np.ones((m, 1)))
        G = matrix(np.vstack((-np.eye(m), np.eye(m))))
        h = matrix(np.hstack((np.zeros(m), np.ones(m) * C)))
        A = matrix(t.reshape(1, -1))
        b = matrix(np.zeros(1))
    
        return P, q, G, h, A, b

\- 此处未作更改

# 12------------------------------------------------- #

Make the function `support_vectors, support_targets, support_alphas, support_ids = compute_support_vectors(X,t,alphas)` to compute:

*  `support_vectors` the support vectors
 *  `support_targets` the target vector associated to the support vectors
 *  `support_alphas` the vector of the dual variables  α i \\alpha\_i αi associated to the support vectors
 *  `support_ids` a boolean vector that indicates if a position is associated to a support vector (w.r.t. the original data matrix) (i.e. if the first instance is a support vector and the second one it is not, then `support_ids` = \[True, False\])

def compute_support_vectors(X, t, alphas):
        t = t.reshape(-1, 1).astype(float)
        support_ids = (alphas > 1e-4).flatten()
        support_vectors = X[support_ids]
        support_targets = t[support_ids]
        support_alphas = alphas[support_ids]
        return support_vectors, support_targets, support_alphas, support_ids

\- 此处未作更改

# 13------------------------------------------------- #

Make the function `model = train_slmkc(X, t, C, kernel_function, param)` to compute the necessary elements to represent a kernelized soft large margin classifier trained over an input data matrix `X` with an associated target vector `t`, for a regularization parameter `C` and a kernel function `kernel_function` with parameter `param`.

Note: you may implement `model` as a tuple containing information on the support vectors and the dual variables, i.e. `model = (support_vectors,support_targets,support_alphas, support_ids)`.

def train_slmkc(X, t, C, kernel_function, param):
        P, q, G, h, A, b = setup_optimization_soft_k(X, t, C, kernel_function, param)
    
        # Setting solver parameters (change default to decrease tolerance)
        solvers.options['show_progress'] = False
        solvers.options['abstol'] = 1e-10
        solvers.options['reltol'] = 1e-10
        solvers.options['feastol'] = 1e-10
    
        # Run solver
        solution_dict = solvers.qp(P, q, G, h, A, b)
        alphas = np.array(solution_dict['x'])
    
        support_vectors, support_targets, support_alphas, support_ids = compute_support_vectors(X, t, alphas)
    
        model = (support_vectors, support_targets, support_alphas, support_ids)
    
        return model

\- 此处未作更改

# 14------------------------------------------------- #

Given the kernelized soft large margin discriminant function of the form:  
 y ( x ) = ∑ n ∈ S V α n t n K ( x n , x ) + b y(x) = \\sum\_\{n \\in SV\} \\alpha\_n t\_n K(x\_n,x) + b y(x)=n∈SV∑αntnK(xn,x)\+b

and given the encoded `model` to represent the necessary elements for a trained kernelized soft large margin classifier

make the function `b = compute_bias(kernel_function, param, C, model)` to compute the offset/bias term `b`, for a regularization parameter `C` and a kernel function `kernel_function` with parameter `param`.

def compute_bias(kernel_function, param, C, model):
        support_vectors, support_targets, support_alphas, support_ids = model
        if param == None:
            w = ((support_targets * support_alphas).T @ support_vectors).reshape(-1, 1)
            b = np.mean(support_targets - np.dot(support_vectors, w))
        else:
            w = None
            b = 0
            for n in range(len(support_targets)):
                b += support_targets[n]
                b -= np.sum(
                    support_alphas * support_targets * kernel_function(support_ids[n], support_alphas, param))
            b /= len(support_alphas)
        return w, b

\- 此处b的计算为重点，对于线性核函数，与教案中计算公式相同。对于非线性核函数或者高斯核函数，应该对上方公式移项，得下方公式计算b值，w值为None，上方因此没有给出w得计算公式。

b = y ( x ) − ∑ n ∈ S V α n t n K ( x n , x ) b = y(x) - \\sum\_\{n \\in SV\} \\alpha\_n t\_n K(x\_n,x) b=y(x)−n∈SV∑αntnK(xn,x)

# 15------------------------------------------------- #

Make the function `scores = score_slmkc(X_test, kernel_function, param, C, model)` to compute the scores for the kernelized soft large margin discriminant function over the data matrix `X_test`, for a regularization parameter `C`, a kernel function `kernel_function` with parameter `param` and an encoded model `model`.

def score_slmkc(X_test, kernel_function, param, C, model):
        w, b = compute_bias(kernel_function, param, C, model)
        if param == None:
            scores = np.dot(X_test, w) + b
            return scores.flatten()
        else:
            support_vectors, support_targets, support_alphas, support_ids = model
            X = X_test
            y_predict = np.zeros(len(X))
            for i in range(len(X)):
                s = 0
                for a, sv_y, sv in zip(support_alphas, support_targets, support_vectors):
                    s += a * sv_y * kernel_function(X[i], sv, param)
                y_predict[i] = s
            return y_predict + b

分为线性与非线性情况，线性下，与教案相同，非线性下，如上述代码所示

# 16------------------------------------------------- #

Make the function `preds = test_slmkc(X_test, kernel_function, param, C, model)` to predict the class label for the kernelized soft large margin discriminant function over the data matrix `X_test`, for a regularization parameter `C`, a kernel function `kernel_function` with parameter `param` and an encoded model `model`.

def test_slmkc(X_test, kernel_function, param, C, model):
        preds = np.sign(score_slmkc(X_test, kernel_function, param, C, model)).astype(int)
        return preds.flatten()

# this function is provided to you, you do not need to alter it in any way.
    class MySoftLargeMarginKernelClassfier(object):
        def __init__(self, C, kernel_function, param):
            self.kernel_function = kernel_function
            self.param = param
            self.C = C
            self.support_vectors = None
            self.support_targets = None
            self.support_alphas = None
            self.support_ids = None
    
        def fit(self, X_train, y_train):
            self.model = train_slmkc(X_train, y_train, self.C, self.kernel_function, self.param)
            self.support_vectors, self.support_targets, self.support_alphas, self.support_ids = self.model
            return self
    
        def predict(self, X_test):
            return test_slmkc(X_test, self.kernel_function, self.param, self.C, self.model)
    
        def decision_function(self, X_test):
            return score_slmkc(X_test, self.kernel_function, self.param, self.C, self.model)
    
        def predict_proba(self, X_test):
            scores = self.decision_function(X_test)
            probs = 1 / (1 + np.exp(-scores))
            return probs

\- 此处未作改变

# 17------------------------------------------------- #

Make the function `linear_kernel_function(x1,x2,param=None)` corresponding to the linear kernel defined as:  
 K ( x 1 , x 2 ) = x 1 T x 2 K(x\_1,x\_2) = x\_1^T x\_2 K(x1,x2)=x1Tx2  
where  x 1 x\_1 x1 and  x 2 x\_2 x2 are vectors.

def linear_kernel_function(x1, x2, param):
        x_prime = x1.T.dot(x2)
        return x_prime

此处为线性核函数

# 18------------------------------------------------- #

Make the function `polynomial_kernel_function(x1,x2,param)` corresponding to the polynomial kernel defined by the parameter  r = r= r=`param`:  
 K ( x 1 , x 2 ) = ( 1 + x 1 T x 2 ) r K(x\_1,x\_2) = (1 + x\_1^T x\_2)^r K(x1,x2)=(1\+x1Tx2)r  
where  x 1 x\_1 x1 and  x 2 x\_2 x2 are vectors.

def polynomial_kernel_function(x1, x2, param):
        result = np.power(np.dot(x1.T, x2) + 1, param)
        return result

def gaussian_kernel(x, y, sigma=5.0):
        return np.exp(-linalg.norm(x - y) ** 2 / (2 * (sigma ** 2)))

此处为多项式核函数与高斯核函数，高斯核函数是我扩充的内容

# 测试部分 #

def plot(X, y, support_vectors=None, predict_func=None):
        cmap = 'bwr'
        res = 100
        x_min, x_max = X[:, 0].min(), X[:, 0].max()
        y_min, y_max = X[:, 1].min(), X[:, 1].max()
        delta = max((x_max - x_min) / 5, (y_max - y_min) / 5)
        x_min -= delta
        y_min -= delta
        x_max += delta
        y_max += delta
        xx, yy = np.meshgrid(np.linspace(x_min, x_max, res), np.linspace(y_min, y_max, res))
        if predict_func is None:
            Z = np.ones(xx.shape) * .5
        else:
            Z = predict_func(np.hstack([xx.reshape(-1, 1), yy.reshape(-1, 1)]))
            if Z.ndim == 1:
                Z = Z.reshape(xx.shape)
            else:
                Z = Z[:, -1].reshape(xx.shape)
        if support_vectors is None:
            edge_colors = 'w'
            sizes = 50
        else:
            edge_colors = ['k' if s else 'w' for s in support_vectors]
            sizes = [100 if s else 50 for s in support_vectors]
        plt.figure(figsize=(5, 5))
        plt.contourf(xx, yy, Z, levels=20, cmap=cmap, alpha=.5)
        plt.scatter(X[:, 0], X[:, 1], c=y, s=sizes, edgecolors=edge_colors, cmap=cmap)
        plt.grid()
        plt.axis('equal')
        plt.show()
    
    
    # Just run the following code, do not modify it
    
    from sklearn.datasets import make_classification
    
    X, y = make_classification(n_samples=100, n_features=2, n_informative=2, n_redundant=0, n_repeated=0, n_classes=2,
                               n_clusters_per_class=2, weights=None, flip_y=0.0, class_sep=1.0, random_state=5)
    y[y == 0] = -1

## 线性核函数 ##

est = MySoftLargeMarginKernelClassfier(C=1e2, kernel_function=linear_kernel_function, param=None).fit(X, y)
    plot(X, y, support_vectors=est.support_ids, predict_func=est.predict_proba)

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L21pcmFjbGVvYQ_size_16_color_FFFFFF_t_70]

## 多项式核函数 ##

est = MySoftLargeMarginKernelClassfier(C=1e-2, kernel_function=polynomial_kernel_function, param=5).fit(X, y)
    plot(X, y, support_vectors=est.support_ids, predict_func=est.predict_proba)

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L21pcmFjbGVvYQ_size_16_color_FFFFFF_t_70 1]

## 高斯核函数 ##

est = MySoftLargeMarginKernelClassfier(C=1e-2, kernel_function=gaussian_kernel, param=5).fit(X, y)
    plot(X, y, support_vectors=est.support_ids, predict_func=est.predict_proba)

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L21pcmFjbGVvYQ_size_16_color_FFFFFF_t_70 2]

# 整体代码 #

import matplotlib.pyplot as plt
    import numpy as np
    import scipy as sp
    import cvxopt
    from cvxopt import matrix
    from cvxopt import solvers
    from numpy import linalg
    
    
    # 10-------------------------------------------------
    def make_H_k(X, t, kernel_function, param):
        n_samples, n_feature = X.shape
        X_prime = np.zeros((n_samples, n_samples))
        for i in range(n_samples):
            for j in range(n_samples):
                X_prime[i, j] = kernel_function(X[i], X[j], param)
        H = np.outer(t, t) * X_prime
        return H
    
    
    # 11-------------------------------------------------
    def setup_optimization_soft_k(X, t, C, kernel_function, param):
        m, n = X.shape
        t = t.reshape(-1, 1).astype(float)
        H = make_H_k(X, t, kernel_function, param)
        P = matrix(H)
        q = matrix(-np.ones((m, 1)))
        G = matrix(np.vstack((-np.eye(m), np.eye(m))))
        h = matrix(np.hstack((np.zeros(m), np.ones(m) * C)))
        A = matrix(t.reshape(1, -1))
        b = matrix(np.zeros(1))
    
        return P, q, G, h, A, b
    
    
    # 12-------------------------------------------------
    def compute_support_vectors(X, t, alphas):
        t = t.reshape(-1, 1).astype(float)
        support_ids = (alphas > 1e-4).flatten()
        support_vectors = X[support_ids]
        support_targets = t[support_ids]
        support_alphas = alphas[support_ids]
        return support_vectors, support_targets, support_alphas, support_ids
    
    
    # 13-------------------------------------------------
    def train_slmkc(X, t, C, kernel_function, param):
        P, q, G, h, A, b = setup_optimization_soft_k(X, t, C, kernel_function, param)
    
        # Setting solver parameters (change default to decrease tolerance)
        solvers.options['show_progress'] = False
        solvers.options['abstol'] = 1e-10
        solvers.options['reltol'] = 1e-10
        solvers.options['feastol'] = 1e-10
    
        # Run solver
        solution_dict = solvers.qp(P, q, G, h, A, b)
        alphas = np.array(solution_dict['x'])
    
        support_vectors, support_targets, support_alphas, support_ids = compute_support_vectors(X, t, alphas)
    
        model = (support_vectors, support_targets, support_alphas, support_ids)
    
        return model
    
    
    # 14-------------------------------------------------
    def compute_bias(kernel_function, param, C, model):
        support_vectors, support_targets, support_alphas, support_ids = model
        if param == None:
            w = ((support_targets * support_alphas).T @ support_vectors).reshape(-1, 1)
            b = np.mean(support_targets - np.dot(support_vectors, w))
        else:
            w = None
            b = 0
            for n in range(len(support_targets)):
                b += support_targets[n]
                b -= np.sum(
                    support_alphas * support_targets * kernel_function(support_ids[n], support_vectors, param))
            b /= len(support_alphas)
        return w, b
    
    
    # 15-------------------------------------------------
    def score_slmkc(X_test, kernel_function, param, C, model):
        w, b = compute_bias(kernel_function, param, C, model)
        if param == None:
            scores = np.dot(X_test, w) + b
            return scores.flatten()
        else:
            support_vectors, support_targets, support_alphas, support_ids = model
            X = X_test
            y_predict = np.zeros(len(X))
            for i in range(len(X)):
                s = 0
                for a, sv_y, sv in zip(support_alphas, support_targets, support_vectors):
                    s += a * sv_y * kernel_function(X[i], sv, param)
                y_predict[i] = s
            return y_predict + b
    
    
    # 16-------------------------------------------------
    def test_slmkc(X_test, kernel_function, param, C, model):
        preds = np.sign(score_slmkc(X_test, kernel_function, param, C, model)).astype(int)
        return preds.flatten()
    
    
    # this function is provided to you, you do not need to alter it in any way.
    
    class MySoftLargeMarginKernelClassfier(object):
        def __init__(self, C, kernel_function, param):
            self.kernel_function = kernel_function
            self.param = param
            self.C = C
            self.support_vectors = None
            self.support_targets = None
            self.support_alphas = None
            self.support_ids = None
    
        def fit(self, X_train, y_train):
            self.model = train_slmkc(X_train, y_train, self.C, self.kernel_function, self.param)
            self.support_vectors, self.support_targets, self.support_alphas, self.support_ids = self.model
            return self
    
        def predict(self, X_test):
            return test_slmkc(X_test, self.kernel_function, self.param, self.C, self.model)
    
        def decision_function(self, X_test):
            return score_slmkc(X_test, self.kernel_function, self.param, self.C, self.model)
    
        def predict_proba(self, X_test):
            scores = self.decision_function(X_test)
            probs = 1 / (1 + np.exp(-scores))
            return probs
    
    
    # 17-------------------------------------------------
    def linear_kernel_function(x1, x2, param):
        x_prime = x1.T.dot(x2)
        return x_prime
    
    
    # 18-------------------------------------------------
    def polynomial_kernel_function(x1, x2, param):
        result = np.power(np.dot(x1.T, x2) + 1, param)
        return result
    
    
    def gaussian_kernel(x, y, sigma=5.0):
        return np.exp(-linalg.norm(x - y) ** 2 / (2 * (sigma ** 2)))
    
    
    # 测试部分
    def plot(X, y, support_vectors=None, predict_func=None):
        cmap = 'bwr'
        res = 100
        x_min, x_max = X[:, 0].min(), X[:, 0].max()
        y_min, y_max = X[:, 1].min(), X[:, 1].max()
        delta = max((x_max - x_min) / 5, (y_max - y_min) / 5)
        x_min -= delta
        y_min -= delta
        x_max += delta
        y_max += delta
        xx, yy = np.meshgrid(np.linspace(x_min, x_max, res), np.linspace(y_min, y_max, res))
        if predict_func is None:
            Z = np.ones(xx.shape) * .5
        else:
            Z = predict_func(np.hstack([xx.reshape(-1, 1), yy.reshape(-1, 1)]))
            if Z.ndim == 1:
                Z = Z.reshape(xx.shape)
            else:
                Z = Z[:, -1].reshape(xx.shape)
        if support_vectors is None:
            edge_colors = 'w'
            sizes = 50
        else:
            edge_colors = ['k' if s else 'w' for s in support_vectors]
            sizes = [100 if s else 50 for s in support_vectors]
        plt.figure(figsize=(5, 5))
        plt.contourf(xx, yy, Z, levels=20, cmap=cmap, alpha=.5)
        plt.scatter(X[:, 0], X[:, 1], c=y, s=sizes, edgecolors=edge_colors, cmap=cmap)
        plt.grid()
        plt.axis('equal')
        plt.show()
    
    
    # Just run the following code, do not modify it
    
    from sklearn.datasets import make_classification
    
    X, y = make_classification(n_samples=100, n_features=2, n_informative=2, n_redundant=0, n_repeated=0, n_classes=2,
                               n_clusters_per_class=2, weights=None, flip_y=0.0, class_sep=1.0, random_state=5)
    y[y == 0] = -1
    
    est = MySoftLargeMarginKernelClassfier(C=1e2, kernel_function=linear_kernel_function, param=None).fit(X, y)
    plot(X, y, support_vectors=est.support_ids, predict_func=est.predict_proba)
    
    est = MySoftLargeMarginKernelClassfier(C=1e-2, kernel_function=polynomial_kernel_function, param=5).fit(X, y)
    plot(X, y, support_vectors=est.support_ids, predict_func=est.predict_proba)
    
    est = MySoftLargeMarginKernelClassfier(C=1e-2, kernel_function=gaussian_kernel, param=5).fit(X, y)
    plot(X, y, support_vectors=est.support_ids, predict_func=est.predict_proba)

[cvxopt]: http://cvxopt.org/
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L21pcmFjbGVvYQ_size_16_color_FFFFFF_t_70]: /images/20221020/bec3df85e2164d07980c58ef15b55ebf.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L21pcmFjbGVvYQ_size_16_color_FFFFFF_t_70 1]: /images/20221020/3ef3e209d4c242bb88947e74e074c27d.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L21pcmFjbGVvYQ_size_16_color_FFFFFF_t_70 2]: /images/20221020/948862e69e1f4ad491fd2bcdbe55f662.png