|
| 1 | +import numpy as np |
| 2 | +from numpy import linalg |
| 3 | +import cvxopt |
| 4 | +import cvxopt.solvers |
| 5 | + |
| 6 | +## define kenrel functions |
| 7 | +def linear_kernel(x1, x2): |
| 8 | + return np.dot(x1, x2) |
| 9 | + |
| 10 | +def polynomial_kernel(x, y, p=3): |
| 11 | + return (1 + np.dot(x, y)) ** p |
| 12 | + |
| 13 | +def gaussian_kernel(x, y, sigma=5.0): |
| 14 | + return np.exp(-linalg.norm(x-y)**2 / (2 * (sigma ** 2))) |
| 15 | +## end define kernel functions |
| 16 | + |
| 17 | +class SVM(object): |
| 18 | + """ |
| 19 | + Suppoet vector classification by quadratic programming |
| 20 | + """ |
| 21 | + |
| 22 | + def __init__(self, kernel=linear_kernel, C=None): |
| 23 | + """ |
| 24 | +
|
| 25 | + :param kernel: kernel types, should be in the kernel function list above |
| 26 | + :param C: |
| 27 | + """ |
| 28 | + self.kernel = kernel |
| 29 | + self.C = C |
| 30 | + if self.C is not None: self.C = float(self.C) |
| 31 | + |
| 32 | + def fit(self, X, y): |
| 33 | + n_samples, n_features = X.shape |
| 34 | + |
| 35 | + # Gram matrix |
| 36 | + K = np.zeros((n_samples, n_samples)) |
| 37 | + for i in range(n_samples): |
| 38 | + for j in range(n_samples): |
| 39 | + K[i,j] = self.kernel(X[i], X[j]) |
| 40 | + |
| 41 | + P = cvxopt.matrix(np.outer(y,y) * K) |
| 42 | + q = cvxopt.matrix(np.ones(n_samples) * -1) |
| 43 | + A = cvxopt.matrix(y, (1,n_samples)) |
| 44 | + b = cvxopt.matrix(0.0) |
| 45 | + |
| 46 | + if self.C is None: |
| 47 | + G = cvxopt.matrix(np.diag(np.ones(n_samples) * -1)) |
| 48 | + h = cvxopt.matrix(np.zeros(n_samples)) |
| 49 | + else: |
| 50 | + tmp1 = np.diag(np.ones(n_samples) * -1) |
| 51 | + tmp2 = np.identity(n_samples) |
| 52 | + G = cvxopt.matrix(np.vstack((tmp1, tmp2))) |
| 53 | + tmp1 = np.zeros(n_samples) |
| 54 | + tmp2 = np.ones(n_samples) * self.C |
| 55 | + h = cvxopt.matrix(np.hstack((tmp1, tmp2))) |
| 56 | + |
| 57 | + # solve QP problem, DOC: http://cvxopt.org/userguide/coneprog.html?highlight=qp#cvxopt.solvers.qp |
| 58 | + solution = cvxopt.solvers.qp(P, q, G, h, A, b) |
| 59 | + |
| 60 | + # Lagrange multipliers |
| 61 | + a = np.ravel(solution['x']) |
| 62 | + |
| 63 | + # Support vectors have non zero lagrange multipliers |
| 64 | + sv = a > 1e-5 |
| 65 | + ind = np.arange(len(a))[sv] |
| 66 | + self.a = a[sv] |
| 67 | + self.sv = X[sv] |
| 68 | + self.sv_y = y[sv] |
| 69 | + print "%d support vectors out of %d points" % (len(self.a), n_samples) |
| 70 | + |
| 71 | + # Intercept |
| 72 | + self.b = 0 |
| 73 | + for n in range(len(self.a)): |
| 74 | + self.b += self.sv_y[n] |
| 75 | + self.b -= np.sum(self.a * self.sv_y * K[ind[n],sv]) |
| 76 | + self.b /= len(self.a) |
| 77 | + |
| 78 | + # Weight vector |
| 79 | + if self.kernel == linear_kernel: |
| 80 | + self.w = np.zeros(n_features) |
| 81 | + for n in range(len(self.a)): |
| 82 | + self.w += self.a[n] * self.sv_y[n] * self.sv[n] |
| 83 | + else: |
| 84 | + self.w = None |
| 85 | + |
| 86 | + def project(self, X): |
| 87 | + if self.w is not None: |
| 88 | + return np.dot(X, self.w) + self.b |
| 89 | + else: |
| 90 | + y_predict = np.zeros(len(X)) |
| 91 | + for i in range(len(X)): |
| 92 | + s = 0 |
| 93 | + for a, sv_y, sv in zip(self.a, self.sv_y, self.sv): |
| 94 | + s += a * sv_y * self.kernel(X[i], sv) |
| 95 | + y_predict[i] = s |
| 96 | + return y_predict + self.b |
| 97 | + |
| 98 | + def predict(self, X): |
| 99 | + return np.sign(self.project(X)) |
| 100 | + |
0 commit comments