第一節：線性支持向量機介紹(1)

支持向量機(Support Vector Machine，SVM)主要被設計來進行「二元分類」問題。

– 他希望有自動進行非線性特徵工程的能力。

進行特徵工程的邏輯與多層感知機不同，多層感知機是模仿神經元進行運算，而SVM是試著找出一個「kernel trick」將低維平面投影至高維平面，用以在高維空間中找出一個「線性可分割」的狀態

F01

SVM的求解需要非常多數學基礎，我們需要花兩節課的時間帶大家介紹整個過程。

– 我們首先要先介紹線性SVM，他是在1963年被創造出來做做線性分割的工作(與邏輯斯回歸相同)，但「kernel trick」是在1992年後才發展出以此方法來解決原始特徵空間「線性不可分割」的問題

第一節：線性支持向量機介紹(2)

讓我們先從「線性SVM」開始介紹，他想做的事情和邏輯斯回歸想法不同，但殊途同歸(結果仍有小差異)。

– 在一個2維平面中，他希望找到一條線能完美的分割紅點與藍點

x1 = c(0, 2, 2, 3)
x2 = c(0, 2, 0, 0)
y = c(1, 1, -1, -1)
plot(x1, x2, col = y + 3, pch = 19, cex = 3)

第一節：線性支持向量機介紹(3)

– 我們當然一眼就能看出有好多條線都能輕易完美分割這四個點，舉例來說，\(x_2 = -0.5 + 0.5 \times x_1\)能幫我們輕易的切開紅藍點

plot(x1, x2, col =  y + 3, pch = 19, cex = 3)
abline(a = -0.5, b = 0.5, lwd = 2, lty = 1)

– 但這樣是不夠的，儘管我們不清楚原因，但你應該覺得\(x_2 = -1 + 1 \times x_1\)這條線切得更好吧！

plot(x1, x2, col =  y + 3, pch = 19, cex = 3)
abline(a = -1, b = 1, lwd = 2, lty = 1)

第一節：線性支持向量機介紹(4)

至於最好的原因呢? 我們認為是因為它離最接近的點最遠

plot(x1, x2, col =  y + 3, pch = 19, cex = 3)
abline(a = -1, b = 1, lwd = 2, lty = 1)
abline(a = 0, b = 1, lwd = 2, lty = 2)
abline(a = -2, b = 1, lwd = 2, lty = 2)

所以我們現在重新定義一下問題，

我們想要找到一條線，能完美切割這些點(約束條件)
這條線離最接近的點要最遠(取極大值)

第一節：線性支持向量機介紹(5)

在這樣的條件之下，我們重新定義我們這條線的式子為\(w_0 + w_1x_1 + w_2x_2 = 0\)為主要粗線

– 在這樣的條線下，\(w_0 = 1\)；\(w_1 = -1\)；\(w_2 = 1\)，這能滿足這個例子上我們要的解

當我們把式子改寫後，我們就能計算每個點到為\(w_0 + w_1x_1 + w_2x_2 = 0\)的距離，而距離的計算可以簡化為下列函數(試著自己推導吧)：

distance.func = function (x1, x2, w0, w1, w2) {
  dist = 1/sqrt(w1^2 + w2^2) * abs(w0 + w1 * x1 + w2 * x2)
  return(dist)
}

distance.func(0, 0, w0 = 1, w1 = -1, w2 = 1)

## [1] 0.7071068

distance.func(2, 2, w0 = 1, w1 = -1, w2 = 1)

## [1] 0.7071068

distance.func(2, 0, w0 = 1, w1 = -1, w2 = 1)

## [1] 0.7071068

distance.func(3, 0, w0 = 1, w1 = -1, w2 = 1)

## [1] 1.414214

第一節：線性支持向量機介紹(6)

但有絕對值很討厭，幸好紅點剛好y為-1，而藍點剛好y為1，所以我們改寫一下方程式

distance.func = function (x1, x2, y, w0, w1, w2) {
  dist = 1/sqrt(w1^2 + w2^2) * y * (w0 + w1 * x1 + w2 * x2)
  return(dist)
}

distance.func(0, 0, 1, w0 = 1, w1 = -1, w2 = 1)

## [1] 0.7071068

distance.func(2, 2, 1, w0 = 1, w1 = -1, w2 = 1)

## [1] 0.7071068

distance.func(2, 0, -1, w0 = 1, w1 = -1, w2 = 1)

## [1] 0.7071068

distance.func(3, 0, -1, w0 = 1, w1 = -1, w2 = 1)

## [1] 1.414214

我們把這個方程式寫下來，並令\(w\)向量為\(\begin{pmatrix} w_1 & w_2 \end{pmatrix}\)，\(x_i\)向量第\(i\)個人的座標向量為\(\begin{pmatrix} x_{1,i} & x_{2,i} \end{pmatrix}\)，而\(b\)為\(w0\)，距離公式可以寫成這樣：

\[ \begin{align} \text{distance}(x_i, b, w) & = \frac{1}{||w||} y_i (w^Tx_i + b) \end{align} \]

第一節：線性支持向量機介紹(7)

到目前為止我們已經能夠求每個點到分割線的距離了，並且也取得距離公式了，我們現在的目標就是希望能夠最大化這個距離\(\text{distance}(x_i, b, w)\)，但並不是對每個點都有興趣，是希望找到距離最小的點，並且把他的距離最大化。

– 所有點的距離可以表達成這樣：

\[ \begin{align} \text{distance}(x_1, b, w) & = \frac{1}{||w||} y_i (w^Tx_1 + b) \\ \text{distance}(x_2, b, w) & = \frac{1}{||w||} y_i (w^Tx_2 + b) \\ \cdots \\ \text{distance}(x_n, b, w) & = \frac{1}{||w||} y_i (w^Tx_n + b) \end{align} \]

– 因此，距離最小的點就可以表達成這樣：

\[ \begin{align} min_{i=1, \cdots , n} \text{distance}(x_i, b, w) = min \frac{1}{||w||} y_i (w^Tx_i + b) \end{align} \]

– 我們希望的是，最大化最小距離的點的距離：

\[ \begin{align} max \left[ min \left[ \frac{1}{||w||} y_i (w^Tx_i + b) \right] \right] \end{align} \]

第一節：線性支持向量機介紹(8)

我們這裡先確定一下你知道這件事情，那就是剛剛我們計算的那組\(w_0 = 1\)；\(w_1 = -1\)；\(w_2 = 1\)，如果同時乘上3倍那答案也會一樣，因此縮放是沒有意義的：

distance.func(0, 0, 1, w0 = 1, w1 = -1, w2 = 1)

## [1] 0.7071068

distance.func(0, 0, 1, w0 = 3, w1 = -3, w2 = 3)

## [1] 0.7071068

– 既然縮放是沒有意義的，我們就可以強迫剛剛圖旁邊的兩條虛線為\(w0 + w_1x_1 + w_2x_2 = 1\)與\(w_0 + w_1x_1 + w_2x_2 = -1\)。

這樣的狀況之下，那剛好在虛線上的點就會是最小的點，而剛好在虛線上的點自然可以被寫成這樣：

\[ \begin{align} w^Tx_i + b & = 1 \text{ or } -1 \\ y(w^Tx_i + b) & = 1 \end{align} \]

既然這是最小的點，那我們當然可以要求虛線中沒有任何點，所以我們可以把需要求解式子改寫一下(約束條件)：

\[ min \left[ y(w^Tx_i + b) \right] = 1 \]

在這樣的條件之下，那既然最小的點帶入距離公式\(\text{distance}(x_i, b, w) = \frac{1}{||w||} y_i (w^Tx_i + b)\)的右半部份\(y_i (w^Tx_i + b)\)為1，那我們又希望距離最大化，這時候的距離是不是就等於\(\frac{1}{||w||} \times 1\)了。

– 我們再重新整理一下求解式子：

\[ \begin{align} & max \frac{1}{||w||} \\ \text{subject to } & min \left[ y(w^Tx_i + b) \right] = 1 \end{align} \]

第一節：線性支持向量機介紹(9)

這個式子還是很難解，其中最討厭的地方是他強迫\(min \left[ y(w^Tx_i + b) \right] = 1\)，但要滿足這個條件，我們可以把式子改寫成這個條件：

\[ \begin{align} & max \frac{1}{||w||} \\ \text{subject to } & y(w^Tx_i + b) \geq 1 \end{align} \]

這跟剛剛的\(y(w^Tx_i + b) \geq 1\)相較於\(min \left[ y(w^Tx_i + b) \right] = 1\)顯然較為寬鬆。

– 假設真的有最小的數值真的大於1為一個數字如1.2，那這樣\(w\)與\(b\)必然就可以同時除以這個數字，導致\(\frac{1}{||w||}\)可以來的更大。

– 因此這個改寫是合法的，當保證\(\frac{1}{||w||}\)最大化的前提之下，我們可以保證條件\(y(w^Tx_i + b) \geq 1\)可以求得與條件\(min \left[ y(w^Tx_i + b) \right] = 1\)相同的解。

因為求最小值比較簡單，我們可以再把上面的式子反過來，從而改寫成：

\[ \begin{align} & min \frac{1}{2}w^Tw \\ \text{subject to } & y(w^Tx_i + b) \geq 1 \end{align} \]

– 在這樣的前提之下，我們試著把不同\(w\)在\(b=1\)時的求解\(loss\)給畫出來：

第一節：線性支持向量機介紹(10)

很遺憾的是，我們沒有辦法用梯度下降法求解，因為這不是一個平滑問題，這是一個線性代數問題，這邊我們先不求他如何解出來

\[ \begin{align} & min \frac{1}{2}w^Tw \\ \text{subject to } & y(w^Tx_i + b) \geq 1 \end{align} \]

我們必須改用二次規劃求解，這個部分我們讓套件幫助我們渡過難關，但我們需要把我們的問題轉化。

F02

\[ \begin{align} & min \text{ } \frac{1}{2}b^TDb - d^Tb \\ \text{ subject to } & A^Tb \geq b_0 \end{align} \]

第一節：線性支持向量機介紹(11)

在這樣比較之下，我們類比這兩個式子(注意上下兩式同樣符號所代表的意義可能不同)：

\[ \begin{align} & min \text{ } \frac{1}{2}w^Tw & \text{ subject to } y(w^Tx_i + b) \geq 1 \\ & min \text{ } \frac{1}{2}b^TDb - d^Tb & \text{ subject to } A^Tb \geq b_0 \end{align} \]

如果我們能強迫其中的向量\(b\)、矩陣\(D\)、向量\(d\)、矩陣\(A\)以及向量\(b_0\)為特定數值，這樣上下兩式就會相等，我們就能求出向量\(b\)的數值：

\[ \begin{align} b = \begin{pmatrix} b \\ w_1 \\ w_2 \end{pmatrix} \space D = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \space d = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix} \space A = \begin{pmatrix} y_1 & y_2 & y_3 & y_4 \\ y_1x_{1,1} & y_2x_{1,2} & y_3x_{1,3} & y_4x_{1,4} \\ y_1x_{2,1} & y_2x_{2,2} & y_3x_{2,3} & y_4x_{2,4} \end{pmatrix} \space b_0 = \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \end{pmatrix} \end{align} \]

library(quadprog)
n.sample = 4
n.weight = 3 
small.value = 5e-6
D.matrix = matrix(small.value, nrow = n.weight, ncol = n.weight)
diag(D.matrix) = 1
D.matrix[1,1] = small.value
A.matrix = rbind(rep(1, n.sample)*y, x1*y, x2*y)
fit = solve.QP(Dmat = D.matrix, dvec = rep(0, n.weight), Amat = A.matrix, bvec = rep(1, n.sample))
fit$solution

## [1]  1 -1  1

第一節：線性支持向量機介紹(12)

– 當答案解出來之後，我們就能利用這個解來畫線了

COEF = fit$solution
A0 = -COEF[1]/COEF[3]
A1 = A0 + 1/COEF[3]
A2 = A0 - 1/COEF[3]
B = -COEF[2]/COEF[3]

plot(x1, x2, col =  y + 3, pch = 19, cex = 3)
abline(a = A0, b = B, lwd = 2, lty = 1)
abline(a = A1, b = B, lwd = 2, lty = 2)
abline(a = A2, b = B, lwd = 2, lty = 2)

– 我們把它寫成一個完整的求解函數

mysvm = function (x1, x2, y) {
  
  require(quadprog)
  
  n.sample = length(x1)
  n.weight = 3
  small.value = 5e-6
  D.matrix = matrix(small.value, nrow = n.weight, ncol = n.weight);diag(D.matrix) = 1; D.matrix[1,1] = small.value
  A.matrix = rbind(rep(1, n.sample)*y, x1*y, x2*y)
  fit = solve.QP(Dmat = D.matrix, dvec = rep(0, n.weight), Amat = A.matrix, bvec = rep(1, n.sample))
  
  COEF = fit$solution
  A0 = -COEF[1]/COEF[3]
  A1 = A0 + 1/COEF[3]
  A2 = A0 - 1/COEF[3]
  B = -COEF[2]/COEF[3]
  
  plot(x1, x2, col =  y + 3, pch = 19)
  abline(a = A0, b = B, lwd = 2, lty = 1)
  abline(a = A1, b = B, lwd = 2, lty = 2)
  abline(a = A2, b = B, lwd = 2, lty = 2)
  
}

– 讓我們試試不同的樣本看看：

mysvm(x1 = c(0, 2, 2, 3),
      x2 = c(0, 2, 0, 0),
      y = c(1, 1, -1, -1))

mysvm(x1 = c(0, 2, 3, 4, 5),
      x2 = c(0, 2, 0, 0, 3),
      y = c(1, 1, -1, -1, -1))

練習1：試著更了解線性SVM

你可能聽到這裡為止只知道公式被推導出來了，但並不知道能拿來做什麼。

– 讓我們先用最簡單的範例資料來進行實驗。

data(iris)
sub.iris <- iris[1:100,]
X1 <- sub.iris[,1]
X2 <- sub.iris[,2]
Y <- as.integer(sub.iris[,5]) * 2 - 3

mysvm(x1 = X1, x2 = X2, y = Y)

請比較線性SVM與邏輯斯回歸的差異

glm.model <- glm((Y>0)~X1+X2, family = "binomial")
COEF <- glm.model$coefficients
A <- -COEF[1]/COEF[3]
B <- -COEF[2]/COEF[3]

plot(X1, X2, col = Y + 3, pch = 19)
abline(a = A, b = B, col = "darkgreen")

幾乎一樣，那我們能不能試試看這個樣本呢?

data(iris)
sub.iris <- iris[51:150,]
X1 <- sub.iris[,1]
X2 <- sub.iris[,2]
Y <- as.integer(sub.iris[,5]) * 2 - 5

練習1答案

這是用邏輯思迴歸得到的解：

glm.model <- glm((Y>0)~X1+X2, family = "binomial")
COEF <- glm.model$coefficients
A <- -COEF[1]/COEF[3]
B <- -COEF[2]/COEF[3]

plot(X1, X2, col = Y + 3, pch = 19)
abline(a = A, b = B, col = "darkgreen")

你應該會發現沒有辦法用「mysvm」做出這個解。

mysvm(x1 = X1, x2 = X2, y = Y)

– 為什麼呢? 因為目前我們的推導過程中要求所有的點都必須滿足\(y(w^Tx_i + b) \geq 1\)這個條件，因此如果找不到完美分割線，那就無解！

– 學了半天你有沒有覺得SVM很難用呀?

第二節：對偶問題(1)

學會基本的線性SVM後，我們仍然無法做「非線性分割」，回頭看看我們利用二次規劃所解的最佳化問題，我們若想要進行「非線性分割」，仍需要手動擴充維度數目，這樣SVM並相較於其他方法仍然沒有任何優勢

\[ \begin{align} & min \frac{1}{2}w^Tw \\ \text{subject to } \space \space \space & y(w^Tx_i + b) \geq 1 \space \space \space \text{ for all } i \end{align} \]

在進入下一個部份之前，我們要先轉換一下上述這個「有條件的」最佳化問題，把他的求解的式子加入一個「拉格朗日乘數(Lagrange multiplier)」，讓問題成為一個「無條件的」最佳化問題。

– 我們先把下面那個條件做一些簡單的轉換：

\[ \begin{align} 1 - y(w^Tx_i + b) \leq 0 \end{align} \]

假定有個數字\(\alpha_i \geq 0\)，那這樣這個\(\alpha_i (1 - y(w^Tx_i + b))\)就必然會小於等於0，

– 那這樣我們可以把優化式子改寫成這樣：

\[ \begin{align} L(b, w,\alpha_i) & = \frac{1}{2}w^Tw + \sum \limits_{i=1}^{n} \alpha_i (1 - y(w^Tx_i + b)) \end{align} \]

第二節：對偶問題(2)

對於新的式子\(L(b, w,\alpha_i)\)，我們的優化目標是找到一組最大的\(\alpha_i\)，再同時找一組\(w\)與\(b\)最小化該式。

\[ \begin{align} L(b, w,\alpha_i) & = \frac{1}{2}w^Tw + \sum \limits_{i=1}^{n} \alpha_i (1 - y(w^Tx_i + b)) \end{align} \]

– 需要注意的是在「滿足條件\(y(w^Tx_i + b) \geq 1\)的狀況下」，由於拉格朗日乘數\(\alpha_i \geq 0\)後面的\(1 - y(w^Tx_i + b) \leq 0\)，因此這個數值必然為負數，我們的優化目標改寫成這樣：

\[ \begin{align} \min \limits_{w, \space b} \left( \max \limits_{\alpha_i \geq 0} \left( L(b, w,\alpha_i) \right) \right) \end{align} \]

這個式子有個特點，那就是由於\(1 - y(w^Tx_i + b) \leq 0\)，因此如果你想要找到一組\(\alpha_i\)最大化該式，那只有滿足條件\(y(w^Tx_i + b) = 1\)的點(這個點會剛好在虛線上)會擁有一個\(\alpha_i > 0\)，剩下的所有不在虛線上的點因為\(y(w^Tx_i + b) > 1\)，所以只能有一個\(\alpha_i = 0\)以保持目標成立。

這是什麼意思呢?那就是求解的時候並不需要考慮所有不在虛線上的點(\(\alpha_i = 0\))，因此等等我們求出來的所有\(\alpha_i > 0\)的點就是畫出線的主要依據，我們把這些「剛好在虛線上的點」稱為「支持向量(support vector)」。

第二節：對偶問題(3)

再讓我們轉換一下我們的問題，現在假定有一組特定的\(\hat{\alpha_i}\)，那下面這個式子一定會滿足下列關係：

\[ \begin{align} \min \limits_{w, \space b} \left( \max \limits_{\alpha_i \geq 0} \left( L(b, w,\alpha_i) \right) \right) \geq \min \limits_{w, \space b} \left( L(b, w,\hat{\alpha_i}) \right) \end{align} \]

– 這應該很容易理解，畢竟在本來的式子中我們最大化了\(\alpha_i\)。

– 順這樣推導下去，那我先「隨便」找一組\(\hat{\alpha_i}\)，最小化該式後再來找一組最大的\(\hat{\alpha_i}\)，那這樣應該也不可能比本來的式子來的更大

\[ \begin{align} \min \limits_{w, \space b} \left( \max \limits_{\alpha_i \geq 0} \left( L(b, w,\alpha_i) \right) \right) \geq \max \limits_{\hat{\alpha_i} \geq 0} \left( \min \limits_{w, \space b} \left( L(b, w,\hat{\alpha_i}) \right) \right) \end{align} \]

這是什麼意思呢，我們把新式\(\max \limits_{\hat{\alpha_i} \geq 0} \left( \min \limits_{w, \space b} \left( L(b, w,\hat{\alpha_i}) \right) \right)\)稱作原式\(\min \limits_{w, \space b} \left( \max \limits_{\alpha_i \geq 0} \left( L(b, w,\alpha_i) \right) \right)\)的對偶問題。

– 如果我們真的能夠解出新式的答案，那它必然會比原式來的更小或相等，但我們不是要最小化原式嗎?因此我們可以將我們的優化目標改成新的式子，而求出來的答案會一模一樣。

– 現在，我們可以直接來解這個式子：

\[ \begin{align} \max \limits_{\hat{\alpha_i} \geq 0} \left( \min \limits_{w, \space b} \left( L(b, w,\hat{\alpha_i}) \right) \right) \end{align} \]

第二節：對偶問題(4)

式子轉換之後有甚麼好處呢?那就是我們將想要求的\(w\)與\(b\)放到了裡面，這樣就可以先用微分求\(w\)與\(b\)的極值的位置有甚麼數學特性：

\[ \begin{align} L(b, w,\alpha_i) & = \frac{1}{2}w^Tw + \sum \limits_{i=1}^{n} \alpha_i (1 - y(w^Tx_i + b)) \\ & = \frac{1}{2}w^Tw + \sum \limits_{i=1}^{n} \left( \alpha_i - \alpha_i yw^Tx_i - \alpha_i y b \right) \\\\ \frac{\partial}{\partial b} L(b, w,\hat{\alpha_i}) & = - \sum \limits_{i=1}^{n} \alpha_i y = 0 \\ \frac{\partial}{\partial w} L(b, w,\hat{\alpha_i}) & = w - \sum \limits_{i=1}^{n} \alpha_i yx_i = 0 \end{align} \]

透過微分方程我們獲得兩條重要的關係式，那就是：

\[ \begin{align} \sum \limits_{i=1}^{n} \alpha_i y & = 0 \space \space \space \space w = \sum \limits_{i=1}^{n} \alpha_i yx_i \end{align} \]

這樣我們可以重整一下原來的式子：

\[ \begin{align} L(b, w,\alpha_i) & = \frac{1}{2}w^Tw + \sum \limits_{i=1}^{n} \left( \alpha_i - \alpha_i yw^Tx_i - \alpha_i y b \right) \\ & = \frac{1}{2}w^Tw + \sum \limits_{i=1}^{n} \alpha_i - \sum \limits_{i=1}^{n} \alpha_i yw^Tx_i - \sum \limits_{i=1}^{n} \alpha_i y b \\ & = \frac{1}{2}w^Tw + \sum \limits_{i=1}^{n} \alpha_i - w^T \sum \limits_{i=1}^{n} \alpha_i yx_i - b \sum \limits_{i=1}^{n} \alpha_i y\\ \left[ \because \sum \limits_{i=1}^{n} \alpha_i y = 0 \right] & = \frac{1}{2}w^Tw + \sum \limits_{i=1}^{n} \alpha_i - w^T \sum \limits_{i=1}^{n} \alpha_i yx_i \\ \left[ \because w = \sum \limits_{i=1}^{n} \alpha_i yx_i \right] & = \frac{1}{2}w^Tw + \sum \limits_{i=1}^{n} \alpha_i - w^T w \\ & = - \frac{1}{2}w^Tw + \sum \limits_{i=1}^{n} \alpha_i\\ & = - \frac{1}{2} \sum \limits_{i=1}^{n} \sum \limits_{j=1}^{n} \alpha_i \alpha_j y_i y_j x_i x_j + \sum \limits_{i=1}^{n} \alpha_i \end{align} \]

第二節：對偶問題(5)

根據剛剛的結果，我們再來看看我們的求解目標，其中裡面的式子我們已經解掉了，可以轉換成下面這樣：

\[ \begin{align} \max \limits_{\hat{\alpha_i} \geq 0} \left( \min \limits_{w, \space b} \left( L(b, w,\alpha_i) \right) \right) & = \max \limits_{\hat{\alpha_i} \geq 0} \left( - \frac{1}{2} \sum \limits_{i=1}^{n} \sum \limits_{j=1}^{n} \alpha_i \alpha_j y_i y_j x_i x_j + \sum \limits_{i=1}^{n} \alpha_i \right) \end{align} \]

通常我們是比較想求最小值的，我們可以把這個問題轉換一下，加上一個原來的條件\(\sum \limits_{i=1}^{n} \alpha_i y = 0\)，我們又可以把這個式子改成二次規劃的形式：

\[ \begin{align} min & \space \space \space \frac{1}{2} \sum \limits_{i=1}^{n} \sum \limits_{j=1}^{n} \alpha_i \alpha_j y_i y_j x_i x_j - \sum \limits_{i=1}^{n} \alpha_i \\ \text{subject to } & \space \space \space \sum \limits_{i=1}^{n} \alpha_i y_i = 0 \space \space \text{and} \space \alpha_i \geq 0 \space \space \space \text{ for all } i \end{align} \]

再跟原來的二次規劃問題比較一下：

\[ \begin{align} min & \space \space \space \space \space \frac{1}{2}b^TDb - d^Tb \\ \text{ subject to} & \space \space \space \space \space A^Tb \geq b_0 \end{align} \]

讓我們來類比一下裡面的每個矩陣：

\[ \begin{align} b = \begin{pmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{pmatrix} \space D =YY^T \otimes (XX^T) = \begin{pmatrix} y_1y_1x_1x_1 & y_1y_2x_1x_2 & y_1y_3x_1x_3 & y_1y_4x_1x_4 \\ y_2y_1x_2x_1 & y_2y_2x_2x_2 & y_2y_3x_2x_3 & y_2y_4x_2x_4 \\ y_3y_1x_3x_1 & y_3y_2x_3x_2 & y_3y_3x_3x_3 & y_3y_4x_3x_4 \\ y_4y_1x_4x_1 & y_4y_2x_4x_2 & y_4y_3x_4x_3 & y_4y_4x_4x_4 \end{pmatrix} \space d = \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \end{pmatrix} \\ \space A = \begin{pmatrix} y_1 & 1 & 0 & 0 & 0 \\ y_2 & 0 & 1 & 0 & 0 \\ y_3 & 0 & 0 & 1 & 0 \\ y_4 & 0 & 0 & 0 & 1 \end{pmatrix} \space b_0 = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \end{align} \]

第二節：對偶問題(6)

我們最終能獲得這樣的語法，注意條件必須從\(A^Tb\)中的第一個數條件是相等\(\sum \limits_{i=1}^{n} \alpha_i y_i = 0\)，這需要下參數【meq = 1：the first meq constraints are treated as equality constraints】：

#set.seed(0)
x1 = c(0, 2, 2, 3)
x2 = c(0, 2, 0, 0)
y = c(1, 1, -1, -1)

library(quadprog)

n.sample = 4
small.value = 5e-6

X = cbind(x1, x2)

D.matrix = (y%*%t(y))*(X%*%t(X))
D.matrix = D.matrix + small.value*diag(n.sample)
A.matrix = t(rbind(matrix(y, ncol = n.sample), diag(n.sample)))

fit = solve.QP(Dmat = D.matrix, dvec = rep(1, n.sample), Amat = A.matrix, bvec = rep(0, n.sample + 1), meq = 1, factorized = FALSE)
qpsol <- fit$solution
print(qpsol)

## [1] 0.4999981 0.4999981 0.9999963 0.0000000

聰明的你是否發現，答案是求出來了，但我們求得的是各點的「拉格朗日乘數(Lagrange multiplier)」，我們並沒有求得我們最想要的W向量

– 不要忘記我們剛剛在解極值時給定的條件\(w = \sum \limits_{i=1}^{n} \alpha_i yx_i\)

– 也不要忘記最開始原始SVM給的條件\(y(w^Tx_i + b) \geq 1\)，其中大於1的部份我們不管，而等於1的部分就是支持向量(\(\alpha_i > 0\)），因此利用支持向量獲得條件\(y(w^Tx_i + b) = 1\)（在做式子消去時要特別注意\(y \in (1, \space -1)\)，所以\(y = \frac {1} {y}\)），那這樣我們就可以獲得下面兩個關係式：

\[ \begin{align} w & = \sum \limits_{i=1}^{n} \alpha_i yx_i \\ b & = yi - w^Tx_i \space \space \space \text{ for all } \alpha_i > 0 \end{align} \]

– 這裡要注意一點，若拉格朗日乘數等於0，在計算b時不能使用，而且也完全不影響到W向量的結果。

findCoefs <- function(a, y, X){
  nonzero <-  abs(a) > 5e-6
  W <- rowSums(sapply(which(nonzero), function(i) a[i]*y[i]*X[i,]))
  b <- mean(sapply(which(nonzero), function(i) y[i]-X[i,]%*%W))
  result <- c(b, W)
  names(result) = c("w0", "w1", "w2")
  return(result)
}

coefs = findCoefs(qpsol, y, X)
print(coefs)

##         w0         w1         w2 
##  0.9999975 -0.9999963  0.9999963

A = -coefs[1]/coefs[3]
B = -coefs[2]/coefs[3]

plot(x1, x2, col =  y + 3, pch = 19)
abline(a = A, b = B, lwd = 2, lty = 1)

第二節：對偶問題(7)

在轉換為對偶問題後，我們試著用原來無法分割的資料來進行實驗，發現居然現在能找到一條線了。

library(quadprog)

data(iris)
sub.iris <- iris[51:150,]
x1 <- sub.iris[,1]
x2 <- sub.iris[,2]
y <- as.integer(sub.iris[,5]) * 2 - 5

n.sample = 100
small.value = 5e-6

X = cbind(x1, x2)

D.matrix = (y%*%t(y))*(X%*%t(X))
D.matrix = D.matrix + small.value*diag(n.sample)
A.matrix = t(rbind(matrix(y, ncol = n.sample), diag(n.sample)))

fit = solve.QP(Dmat = D.matrix, dvec = rep(1, n.sample), Amat = A.matrix, bvec = rep(0, n.sample + 1), meq = 1, factorized = FALSE)
qpsol <- fit$solution

coefs = findCoefs(qpsol, y, X)

A = -coefs[1]/coefs[3]
B = -coefs[2]/coefs[3]

plot(x1, x2, col =  y + 3, pch = 19)
abline(a = A, b = B, lwd = 2, lty = 1)

但這個過程其實是有問題的，畢竟剛剛數學推導時是不允許這種事情發生的。但結果是好的，又該如何解釋呢?

– 先記住這個結果，我們下節課會介紹軟邊界的SVM，就會從數學上介紹這到底是怎麼回事。

第三節：利用套件做出SVM(1)

到了這一步，我們終於可以運用套件來做SVM了。

– 自己寫函數雖然對學習它的原理很有幫助，但畢竟不是好的解決之道(為何要重新發明輪子?)，套件「e1071」是在R裡面做SVM計算最常用的套件，讓我們用它來對iris分類一下

– 需要注意的是這裡有個參數【cost】，這是「違反\(y(w^Tx_i + b) \geq 1\)所要付出的代價」，我們目前完全不允許他違反，所以先設一個非常大的代價，下節課我們會介紹如何將這個參數帶入SVM：

library(e1071)

x1 = c(0, 2, 2, 3)
x2 = c(0, 2, 0, 0)
y = c(1, 1, -1, -1)

svm.model = svm(y ~ x1 + x2, kernel = "linear", scale = FALSE, type = "C-classification", cost = 1e5)

讓我們來解讀一下結果的數字(請參考函數說明)

– 這就是「支持向量」(所有\(\alpha_i > 0\)的點)

svm.model$SV

##   x1 x2
## 1  0  0
## 2  2  2
## 3  2  0

– 這是每個「支持向量」的拉格朗日乘數乘以標籤\(\alpha_i y_i\)

svm.model$coefs

##      [,1]
## [1,]  0.5
## [2,]  0.5
## [3,] -1.0

– 這是負數截距\(-b\)

svm.model$rho

## [1] -1

我們可以透過支持向量求得w1及w2的數字

W.vector = rowSums(sapply(1:length(svm.model$coefs), function(i) svm.model$coefs[i]*svm.model$SV[i,]))
w0 = -svm.model$rho
w1 = W.vector[1]
w2 = W.vector[2]

看來我們又能利用w0、w1、w2畫出分割線了

A0 = -w0/w2
A1 = A0 + 1/w2
A2 = A0 - 1/w2
B = -w1/w2

plot(x1, x2, col =  y + 3, pch = 19, cex = 3)
abline(a = A0, b = B, lwd = 2, lty = 1)
abline(a = A1, b = B, lwd = 2, lty = 2)
abline(a = A2, b = B, lwd = 2, lty = 2)

這是與分割線的「垂直距離」(非最短距離)，可以用來畫ROC曲線

svm.model$decision.values

##   1/-1
## 1    1
## 2    1
## 3   -1
## 4   -2

第三節：利用套件做出SVM(2)

讓我們試著用點不同的資料吧！

data(iris)
sub.iris = iris[1:100,c(1, 2, 5)]
sub.iris[,3] = as.factor(as.character(sub.iris[,3]))
X1 = sub.iris[,1]
X2 = sub.iris[,2]
Y = as.integer(sub.iris[,3])

svm.model = svm(Species ~ Sepal.Length + Sepal.Width, data = sub.iris, kernel = "linear", scale = FALSE, type = "C-classification", cost = 1e5)

W.vector = rowSums(sapply(1:length(svm.model$coefs), function(i) svm.model$coefs[i]*svm.model$SV[i,]))
w0 = -svm.model$rho
w1 = W.vector[1]
w2 = W.vector[2]

A0 = -w0/w2
A1 = A0 + 1/w2
A2 = A0 - 1/w2
B = -w1/w2

plot(X1, X2, col =  Y * 2 + 2, pch = 19)
abline(a = A0, b = B, lwd = 2, lty = 1)
abline(a = A1, b = B, lwd = 2, lty = 2)
abline(a = A2, b = B, lwd = 2, lty = 2)

利用函數「predict」可以用SVM MODEL來預測新資料

predict(svm.model, data.frame(Sepal.Length = 5, Sepal.Width = 4))

##      1 
## setosa 
## Levels: setosa versicolor

predict(svm.model, data.frame(Sepal.Length = 6, Sepal.Width = 3))

##          1 
## versicolor 
## Levels: setosa versicolor

用內建的函數來畫畫吧

plot(svm.model, sub.iris)

練習2：重現套件預測機制

你已經大約知道線性SVM的原理了，現在我希望你能重現函數「predict」的預測結果，但不使用該函數

– 關鍵是你如何求得「w0、w1、w2」

– 提示：先想想「decision.values」是怎樣算出來的

data(iris)
sub.iris = iris[1:100,c(1, 2, 5)]
sub.iris[,3] = as.factor(as.character(sub.iris[,3]))
X1 = sub.iris[,1]
X2 = sub.iris[,2]
Y = as.integer(sub.iris[,3])

svm.model = svm(Species ~ Sepal.Length + Sepal.Width, data = sub.iris, kernel = "linear", scale = FALSE, type = "C-classification", cost = 1e5)

W.vector = rowSums(sapply(1:length(svm.model$coefs), function(i) svm.model$coefs[i]*svm.model$SV[i,]))
w0 = -svm.model$rho
w1 = W.vector[1]
w2 = W.vector[2]

predict(svm.model, data.frame(Sepal.Length = 5, Sepal.Width = 4), decision.values = TRUE)

##      1 
## setosa 
## attr(,"decision.values")
##   setosa/versicolor
## 1          6.791608
## Levels: setosa versicolor

練習2答案

答案就在剛剛的推導以及程式碼，我們需要把「w0、w1、w2」都求出來。

– 使用Sepal.Length = 5以及Sepal.Width = 4的點。

W.vector = rowSums(sapply(1:length(svm.model$coefs), function(i) svm.model$coefs[i]*svm.model$SV[i,]))
w0 = -svm.model$rho
w1 = W.vector[1]
w2 = W.vector[2]

w1 * 5 + w2 * 4 + w0

## Sepal.Length 
##     6.791608

你是不是更了解SVM的運作原理了？

課程小結

支持向量機是一個強大的分類模型，其背後的數學原理非常複雜，而目前為止的支持向量機還非常不好用。

– 目前的支持向量機要求所有的點都必須滿足\(y(w^Tx_i + b) \geq 1\)，也就是我們只能求有完美線性分割解的問題，而這個邏輯斯回歸也能解出一樣的結果，並且邏輯斯回歸還能面對不能完美分割的情形。

– 我們還是沒有解決高維空間投影的部分，現在如果想要多幾個變項還是需要手動進行特徵工程，有點不方便。

雖然他還存在上述的問題，但透過對偶問題的求解我們發現現在我們求解的關鍵已經變成在求拉格朗日乘數(\(\alpha_i\))了，並且注意我們在求二次規劃時所用的關鍵矩陣\(D\)的維度已經從「變項數」變成了「樣本數」，而這是很重要的關鍵。

– 透過這個關鍵的轉換，我們才有可能面對「無限多維」的高維空間，我們下節課會繼續介紹支持向量機的擴展，讓你感受這個強大的分類器為什麼會在深度學習突破到來之前統治世界。

– 至今為止，支持向量機可能還是在數值分析中最好用的工具之一。

機器學習及演算法