Xin chào các bạn,
Trong bài post này, mình sẽ giải thích về toán học và code Support Vector Machine.
1. Giới thiệu
Support Vector Machine là một giải thuật phân loại có giám sát (supervised learning) được đề xuất bởi nhà khoa học người Liên Xô Vladimir Vapnik và các cộng sự. Đây là một giải thuật cực kì thú vị cho các bạn mới học ML/AI vì cách diễn giải và chứng minh toán học không quá phức tạp nhưng cực kì hiệu quả.
Ban đầu, giải thuật này chỉ được sử dụng cho việc phân loại 2 classes (binary classification). Tuy nhiên, sau đó cộng đồng AI/ML đã cải biến để có thể dùng nó cho multiple-class classification, regression.
2. Giải thích toán học
Các lý thuyết về của SVM bao quanh các định lý cơ bản của toán học như vector, matrix, tối ưu Larange.
Primal form của SVM có dạng như sau:
\[\min_{w, b} \frac{1}{2}||w|| + C \sum_{1}^{n}{\xi_i}\]
\[\text{s.t. } y_i(w^T\phi(x_i) + b) \ge 1 - \xi_i, \forall i\]
\[\xi_i \ge 0, \forall i\]
Tới đây, nếu không vướng phải hàm \(\phi\) thì ta có thể dùng các thư viện tối ưu để tìm ra được \(w\) và \(b\). Hàm \(\phi\) là một hàm mà sẽ chiếu dữ liệu gốc lên một không gian mới (có thể nhiều chiều hơn hoặc ít chiều hơn), tuy nhiên chúng ta không biết hàm này và sẽ phải chọn thủ công để tìm được hàm tối ưu nhất cho từng task.
Có một cách mà có thể bỏ qua được việc biết hàm \(\phi\) là gì, đó là biến đổi hàm tối ưu này sang dạng dual và dùng phương pháp kernelisation (mình sẽ nói rõ hơn ở phần áp dụng nó vào).
Bây giờ, chúng ta sẽ từng bước biến đổi từ bài toán primal sang dual. Ở bài toán primal, ta sẽ áp dụng phương pháp Lagrange và được như sau:
\[L(w, b, [{\xi_i}], [\alpha_i], [\lambda_i]) = \frac{1}{2}||w|| + \sum_{i=1}^{n}\alpha_i(1 - \xi_i - y_i(w^T\phi(x_i) + b)) + \sum_{i=1}^{n}\lambda_i \xi_i\]
Đây là hàm Lagrange, hàm này biến đổi từ một mớ constraints thành một hàm unconstraints. Để tìm được điểm tối ưu, chúng ta phải giải phương trình sau.
\[L^*(w, b, [{\xi_i}], [\alpha_i], [\lambda_i]) = \inf_{w} \sup_{[\alpha_i], [\lambda_i]}\frac{1}{2}||w|| + \sum_{i=1}^{n}\alpha_i(1 - \xi_i - y_i(w^T\phi(x_i) + b)) + \sum_{i=1}^{n}\lambda_i \xi_i\]
Bình thường, để có thể biến đổi từ dạng primal sang dual mà không làm biến đổi kết quả cuối cùng của bài toán, chúng ta phải chứng minh chúng đối ngẫu (strong duality). Tuy nhiên, bài toán này đã được chứng minh là đối ngẫu nếu primal có nghiệm (chắc chắn primal của bài toán này có nghiệm), và mình sẽ chứng minh nó đối ngẫu ở các bài riêng. Tạm thời, chúng ta hãy chấp nhận strong duality và tiếp tục biến đổi.
Với strong duality, ta có hệ quả như sau:
\[L^* = \inf_{w, b, [\xi_i]} \sup_{[\alpha_i], [\lambda_i]} = \sup_{[\alpha_i], [\lambda_i]} \inf_{w, b, [\xi_i]}\]
Áp dụng hệ quả của strong duality, ta được phương trình sau
\[L^*(w, b, [{\xi_i}], [\alpha_i], [\lambda_i]) = \sup_{[\alpha_i], [\lambda_i]} \inf_{w, b, [\xi_i]} \frac{1}{2}||w|| + \sum_{i=1}^{n}\alpha_i(1 - \xi_i - y_i(w^T\phi(x_i) + b)) + \sum_{i=1}^{n}\lambda_i \xi_i\]
Với phương trình trên, ta sẽ tìm \(w\) và \(b\) sao cho hàm Lagrange có giá trị nhỏ nhất. Ta dùng đạo hàm để tìm giá trị cực tiểu.
\[\frac{\partial L}{\partial w} = 0 <=> w - \sum_{i=1}^{n}\alpha_n y_n \phi(x_n) = 0 <=> w = \sum_{i=1}^{n}\alpha_n y_n \phi(x_n)\]
\[\frac{\partial L}{\partial b} = 0 <=> \sum_{i=1}^{n}\alpha_n y_n\]
\[\frac{\partial L}{\partial \xi} = 0 <=> C - \bold{\alpha} - \bold{\lambda} = 0\]
Áp dụng các phương trình trên vào hàm Lagrange, ta được:
\[L([\alpha_i], [\lambda_i]) = \sum_{i=1}^{n} \alpha_i + \frac{1}{2} \sum_{i=1}^{m}\sum_{i=1}^{n} \alpha_m \alpha_n y_m y_n \phi^T(x_m)\phi(x_n)\]
\[\text{s.t. } \alpha_i, \lambda_i \ge 0\]
\[\sum_{i=1}^{n}\alpha_i y_i = 0\]
\[C - \alpha - \lambda = 0\]
Như quan sát, trong làm Lagrange bây giờ không còn tồn tại \(\lambda\) nhưng \(\lambda\) vẫn phải lớn hơn 0 do ràng buộc Lagrange. Ta sẽ dùng nó để biến đổi tiếp như sau:
\[L([\alpha_i], [\lambda_i]) = \sum_{i=1}^{n} \alpha_i + \frac{1}{2} \sum_{i=1}^{m}\sum_{i=1}^{n} \alpha_m \alpha_n y_m y_n \phi^T(x_m)\phi(x_n)\]
\[\text{s.t. } \alpha_i \ge 0\]
\[\sum_{i=1}^{n}\alpha_i y_i = 0\]
\[C - \alpha \ge 0\]
Và cuối cùng, ta cần phải giải phương trình với các constraints như sau để có thể tìm ra nghiệm của bài toán.
\[L^*([\alpha_i], [\lambda_i]) = \sup_{[\alpha_i]} \sum_{i=1}^{n} \alpha_i + \frac{1}{2} \sum_{i=1}^{m}\sum_{i=1}^{n} \alpha_m \alpha_n y_m y_n \phi^T(x_m)\phi(x_n)\]
\[\text{s.t. } 0 \le \alpha_i \le C\]
\[\sum_{i=1}^{n}\alpha_i y_i = 0\]
Nhìn có vẻ phức tạp, nhưng thực chất nó có thể được biểu diễn như bài toán quadratic convex optimisation và có thể được giải dễ dàng. Tuy nhiên, chúng ta vẫn còn bị cấn ở phần hàm \(\phi\) và mục đích biến đổi nãy giờ chỉ để bypass phần này. Để có thể thực hiện việc này, các nhà khoa học đã đề xuất một phương pháp được gọi là kernelisation trick.
Giả sử, hàm \(\phi \in \mathbb{R}^2 \to \mathbb{R}^3\) có dạng như sau:
\[\phi: x \to \phi(x) = \begin{bmatrix} x_1^2 \\ \sqrt{2}x_1x_2 \\ x_2^2 \end{bmatrix}\]
\[\phi^T(x_m) \cdot \phi(x_n) = \begin{bmatrix} x_{m1}^2 \\ \sqrt{2}x_{m1}x_{m2} \\ x_{m2}^2 \end{bmatrix}^T \cdot \begin{bmatrix} x_{n1}^2 \\ \sqrt{2}x_{n1}x_{n2} \\ x_{n2}^2 \end{bmatrix}\]
\[<=> \phi^T(x_m) \cdot \phi(x_n) = x_{m1}^2x_{n1}^2 + 2x_{m1}x_{m2}x_{n1}x_{n2} + x_{m2}^2x_{n2}^2\]
\[<=> \phi^T(x_m) \cdot \phi(x_n) = (x_{m1}x_{n1} + x_{m2}x_{n2})^2 = (x_m^T \cdot x_n)^2\]
Như các bạn thấy, sau khi dot 2 vectors lại thì ta được một kết quả là dot của 2 vectors gốc không liên quan gì đến \(\phi\) mà phụ thuộc vào kết quả của một cách kết hợp của 2 dữ liệu ở chiều gốc và người ta đặt tên cho nó là kernelisation.
Có rất nhiều loại kernel, một vài loại phổ biến là:
\[\text{Polynomial: } k(x_m, x_n) = (x_m^Tx_n + r)^d\]
\[\text{Linear: } k(x_m, x_n) = x_m^Tx_n\]
\[\text{Gaussian: } k(x_m, x_n) = e^{\frac{||x_m - x_n||^2}{2\sigma^2}}\]
\[\text{Laplace: } k(x_m, x_n) = e^{-\alpha||x_m - x_n||}\]
Nhưng mình thấy dùng phổ biến nhất có Gaussian. Việc kernel tốt và phù hợp với bài toán sẽ là quyết định của người sử dụng.
Áp dụng kernelisation vào phương trình Larange ở trên, ta được:
\[L^*([\alpha_i], [\lambda_i]) = \sup_{[\alpha_i]} \sum_{i=1}^{n} \alpha_i + \frac{1}{2} \sum_{i=1}^{m}\sum_{i=1}^{n} \alpha_m \alpha_n y_m y_n k(x_m, x_n)\]
\[\text{s.t. } 0 \le \alpha_i \le C\]
\[\sum_{i=1}^{n}\alpha_i y_i = 0\]
Sau khi dùng quadratic programming tìm ra nghiệm tối ưu cho phương trình trên, ta sẽ được \(\alpha^*\) và sẽ dùng nó để tìm \(w\) và \(b\).
Áp dụng phương trình từ đạo hàm Larange lúc nãy, ta có \(w\) và \(b\) được tính như sau:
\[w^* = \sum_{i=1}^{n}\alpha_i^* y_n \phi(x_n)\]
Dựa vào phương trình trên, ta quan sát thấy \(w^*\) chỉ bị ảnh hưởng khi \(\alpha_i > 0\) và sẽ không bị ảnh hưởng khi \(\alpha_i = 0\). Ta gọi những training data \(x_i\) có \(\alpha > 0\) là support vectors do chỉ các training data này mới có tác động.
Nhớ lại ở đầu mục, chúng ta đã set điều kiện \(w^T\phi(x_{sv}) + b = y_i = \pm 1\). Vì vậy sau khi biết \(w\) và biết được các support vectors, ta sẽ áp dụng để tìm \(b\).
\[b = y_i - {w^*}^T \phi(x_{SV}) = y_i - \sum_{i=1}^{n}\alpha_i^* y_n k(x_n, x_{SV})\]
Chỉ cần chọn một support vector là chúng ta có thể tìm ra \(b\).
Và sau khi tìm ra \(w\) và \(b\), chúng ta đã tìm ra được hyperplane tối ưu nhất để phân tách được 2 class.
3. Code SVM
Ở phần này, mình sẽ code SVM chỉ dùng các thư viện numpy (tính toán ma trận), cvxopt (Tối ưu), matplotlib (visualise)
- Bước 1: Import các thư viện cần thiết
import numpy as np
import cvxopt
import matplotlib.pyplot as plt
- Bước 2: Viết các kernelisation functions
Có nhiều hàm kernelisation. Tuy nhiên, trong nội dung bài post, mình sẽ chỉ dùng gaussian (hay còn được gọi là RBF) và đây cũng là kernel phổ biến nhất khi SVM được dùng.
def linear(x, z):
return np.dot(x, z.T)
def polynomial(x, z, p=5):
return (1 + np.dot(x, z.T))**p
def gaussian(x, z, sigma=0.1):
return np.exp(-np.linalg.norm(x-z, axis=1)**2/(2*(sigma**2)))
- **Bước 3: Viết code training, predict **
class SVM:
def __init__(self, kernel=gaussian, C=1):
self.kernel = kernel
self.C = C
def fit(self, X, y):
self.y = y
self.X = X
m, n = X.shape
# Calculate Kernel
self.K = np.zeros((m, m))
for i in range(m):
self.K[i, :] = self.kernel(X[i, np.newaxis], self.X)
# Solve with cvxopt final QP needs to be reformulated
# to match the input form for cvxopt.solvers.qp
P = cvxopt.matrix(np.outer(y, y) * self.K)
q = cvxopt.matrix(-np.ones((m, 1)))
G = cvxopt.matrix(np.vstack((np.eye(m) * -1, np.eye(m))))
h = cvxopt.matrix(np.hstack((np.zeros(m), np.ones(m) * self.C)))
A = cvxopt.matrix(y, (1, m), "d")
b = cvxopt.matrix(np.zeros(1))
cvxopt.solvers.options["show_progress"] = False
sol = cvxopt.solvers.qp(P, q, G, h, A, b)
self.alphas = np.array(sol["x"])
def predict(self, X):
y_predict = np.zeros((X.shape[0]))
sv = self.get_parameters(self.alphas)
for i in range(X.shape[0]):
y_predict[i] = np.sum(
self.alphas[sv]
* self.y[sv, np.newaxis]
* self.kernel(X[i], self.X[sv])[:, np.newaxis]
)
return np.sign(y_predict + self.b)
def get_parameters(self, alphas):
threshold = 1e-5
sv = ((alphas > threshold) * (alphas < self.C)).flatten()
self.w = np.dot(self.X[sv].T, alphas[sv] * self.y[sv, np.newaxis])
self.b = np.mean(
self.y[sv, np.newaxis]
- self.alphas[sv] * self.y[sv, np.newaxis] * self.K[sv, sv][:, np.newaxis]
)
return sv
- Bước 4: Viết các hàm helpers:
def create_dataset(N, D=2, K=2):
X = np.zeros((N * K, D)) # data matrix (each row = single example)
y = np.zeros(N * K) # class labels
for j in range(K):
ix = range(N * j, N * (j + 1))
r = np.linspace(0.0, 1, N) # radius
t = np.linspace(j * 4, (j + 1) * 4, N) + np.random.randn(N) * 0.2 # theta
X[ix] = np.c_[r * np.sin(t), r * np.cos(t)]
y[ix] = j
# lets visualize the data:
plt.scatter(X[:, 0], X[:, 1], c=y, s=40, cmap=plt.cm.Spectral)
plt.show()
y[y == 0] -= 1
return X, y
def plot_contour(X, y, svm):
# plot the resulting classifier
h = 0.01
x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
points = np.c_[xx.ravel(), yy.ravel()]
Z = svm.predict(points)
Z = Z.reshape(xx.shape)
plt.contourf(xx, yy, Z, cmap=plt.cm.Spectral, alpha=0.8)
# plt the points
plt.scatter(X[:, 0], X[:, 1], c=y, s=40, cmap=plt.cm.Spectral)
plt.show()
- Bước 4: Tạo data giả và train SVM để đánh giá
np.random.seed(2)
X, y = create_dataset(N=1000)
permuted_idxes = np.random.permutation(len(X))
X = X[permuted_idxes]
y = y[permuted_idxes]
val_split = 50
X_train, Y_train = X[:val_split], y[:val_split]
X_val, Y_val = X[val_split:], y[val_split:]
svm = SVM(kernel=gaussian, C=0.5)
svm.fit(X_train, Y_train)
y_pred = svm.predict(X_val)
print(f"Accuracy: {sum(Y_val==y_pred)/Y_val.shape[0]}")
Accuracy: 0.9384615384615385
from sklearn.svm import SVC
clf = SVC()
clf.fit(X_train, Y_train)
pred = clf.predict(X_val)
print(f"Accuracy: {len(pred[pred==Y_val])/(Y_val.shape[0])}")
Accuracy: 0.981025641025641
4. Kết luận
Support Vector Machine là một giải thuật binary classification được hình thành dựa trên các định lý và công thức toán học về vector, matrix, và tối ưu. Giải thuật này đã được ra đời khá lâu và hiện tại đã có rất nhiều giải thuật machine learning khác cho kết quả vượt trội hơn, tuy nhiên các giá trị toán học nền tảng cho người mới bắt đầu là rất, rất nhiều. Vì vậy, ở các trường đại học trong các môn liên quan đến machine learning, SVM vẫn luôn được dạy.
Hy vọng các bạn cảm thấy bài viết hữu ích. Nếu có chỗ nào chưa rõ, hãy email và mình sẽ trả lời. Chúc các bạn học tốt.
References
1. SVM playlist (Youtube) - Prof. Cynthia Rudin
2. Support Vector Machines - THE MATH YOU SHOULD KNOW (Youtube)
3. Understanding the mathematics behind Support Vector Machines