《Machine Learning in Action》—— 剖析支援向量機，優化SMO

玩世不恭的Coder發表於2020-11-16
Mac優化

《Machine Learning in Action》—— 剖析支援向量機，優化SMO

薄霧濃雲愁永晝，瑞腦銷金獸。

愁的很，上次不是更新了一篇關於支援向量機的文章嘛，《Machine Learning in Action》—— 剖析支援向量機，單手狂撕線性SVM。雖然效果還算不錯，資料集基本都能夠分類正確，模型訓練效率的話也還說的過去，但這是基於我們訓練樣本資料集比較少、迭代次數比較少的前提下。

假如說我們資料集比較大，而且還需要迭代不少次數的話，上一篇文章中使用到的SMO演算法的效率可就不敢恭維了，訓練的速度可堪比龜龜。月黑風高夜，殺人放火天。不對不對，月黑風高夜，瘋狂肝文時。既然一般的SMO演算法的效率低下，那怎麼說也得進一步優化才行吶。

就在前幾天還聽見收音機上說，國內外有許多如雷貫耳的大佬都在不斷研究新演算法來進一步提高SMO演算法的訓練效率。聞此一言，Taoye心中大喜："如果我能蹦躂出一個新的優化演算法，哥哥我聲名遠揚的大好機會就來了啊，雄霸天下的抱負就指日可待了啊！哈哈哈哈！！！"想法雖好，可是該怎麼優化呢？在這薄霧濃雲、月黑風高之夜，Taoye的思緒漫天飛，愁的很。

《Machine Learning in Action》—— 剖析支援向量機，優化SMO

有心栽花花不開，無心插柳柳成蔭。

明明已經知道SMO演算法待優化的地方太多了，可是就是不知如何下手，想了老半天，腦闊疼的厲害。此刻實驗室空無一人，Taoye轉著座椅，雙目望向窗外，皎潔的月光總是給人無限遐想。

罷了罷了，與其木訥在這有心栽花花不開，倒不如出去轉悠轉悠，說不定能捕獲個意想不到的收穫，給我來了無心插柳柳成蔭呢？說走我們就走哇（調調有點不對勁~~），關上空調，披件外套，鎖上室門，雙手插袋朝外走去。

或許是空調吹太久了，亦或實驗室呆太久了，出來的瞬間一股神清氣爽的感覺湧上心頭，五音不全的我此時還真想高歌一曲。頃之，微涼，好在出門前披了件外套。活動活動筋骨，朝湖邊走去。走著走著，不知不覺來到了步行橋，風平浪靜的湖面，沒有一絲波紋蕩起。左轉，低頭看著湖面中胡子拉碴且憔悴的自己，此時我的眼角又再一次。。。┭┮﹏┭┮ 。。抬頭看著湖邊零星幾對小情侶，或有說有笑、或呢喃竊語、或打情罵俏，滋滋滋，有點兒意思，只有我只身一人還在想著如何優化SMO演算法。

等會兒。。。小情侶？？？優化SMO？？？

我記得在前面那篇文章中寫到的SMO演算法的核心思想裡，正是不斷迭代成雙成對的\(\alpha_i\)和\(\alpha_j\)，只不過那個時候的這對“小情侶”是隨意配對，所以導致的排列組合的可能性太多，從而拉低了整個模型訓練的效率。假如說，我那個時候選取這對“小情侶”的時候並不是無意的，而是有意識、有條件的去選取，不就可以避開大量沒必要的可能性計算麼，這樣一來不就可以大大提高模型的訓練效率了麼？？？

握了一棵草，乖乖，我真是個天才，這都被我想到了？？？心頭靈光一閃，猶如飢渴春筍聽到的第一聲驚雷，屁顛屁顛的朝實驗室跑去，而一對對小情侶異樣且詫異的眼神朝我看來。。。

上述內容部分虛構，僅做引出下文之用。

在這之前，我們先靜下心來分析一下上篇文章中SMO的核心演算法：

上圖是我們在講解手撕線性SVM中所寫到的linear_smo方法，詳細請看：《Machine Learning in Action》—— 剖析支援向量機，單手狂撕線性SVM。其中Taoye圈出了三個程式碼塊來給大家介紹下：

第一個程式碼塊，我們可以發現程式碼行為for i in range(m):，想必大家都知道這是一個迴圈語句，在這個方法中它具體表達的意思是根據樣本數量一次選取索引\(i\)，然後通過這個索引來確定\(\alpha_i\)的選擇，所以它最終會把所有的樣本都“走一遍”。

第二個程式碼塊，是根據\(i\)的值重新在\(m\)中選取一個不與\(i\)相同的\(j\)，然後根據這個\(j\)來修改對應的\(\alpha_j\)

第三個就是我們大量的矩陣、向量進行計算的程式碼塊了，我們可以發現，無論前兩個\(i、j\)的選取是怎樣的，第三個程式碼塊都會去執行、計算，然而有些計算完全是沒必要的，這樣就大大拉跨了整個SMO演算法的效率，這可不是我們想要的。

綜上，我們需要在\(\alpha\)的選取上做點文章，使其在一定的跳過第三個程式碼塊的計算。

  • 第一個\(\alpha_i\)的選取

在上一篇文章中，我們也有提到，大多數樣本對於決策面的確定都是無用的，只有少數部分的樣本點才能確定具體的決策面。而\(\alpha\)與樣本之間滿足如下關係：

\[ \left\{ \begin{array}{c} \alpha_i=0 \quad <=> \quad y_i(w^T+b)\geq1\\ \alpha_i=C \quad <=> \quad y_i(w^T+b)\leq1\\ 0\leq\alpha_i\leq C \quad <=> \quad y_i(w^T+b)=1\\ \end{array} \right. \]

對於第一個\(\alpha_i\)的選取，初始化的\(\alpha\)向量為0，所以第一次迭代是對所有的樣本點進行。待第一次迭代完成之後，此時的\(\alpha\)向量已經更新完成了，在之後的迭代過程中，我們就不需要對所有的樣本進行遍歷了，而是選取在\((0, C)\)區間的\(\alpha_i\)值即可，因為其他的\(\alpha\)值對於最終決策面的確定沒有什麼影響。

對於第一個\(\alpha_i\)的選取，核心程式碼思想如下所示，也就是我們的外層迴圈。其中data_struct是一個資料結構，其內部儲存了一些公有地的屬性，這個我們後面會講：

"""
    Author：Taoye
    微信公眾號：玩世不恭的Coder
    Explain：外層迴圈，選取第一個合適alpha
    Parameters：
        x_data: 樣本屬性特徵矩陣
        y_label: 屬性特徵對應的標籤
        C：懲罰引數
        toler：容錯率
        max_iter：迭代次數
    Return:
        b: 決策面的引數b
        alphas：獲取決策面引數w所需要的alphas
"""
def outer_smo(self, data_struct, x_data, y_label, C, toler, max_iter):
	iter_num, ergodic_flag, alpha_optimization_num = 0, True, 0
	while (iter_num < max_iter) and ((alpha_optimization_num > 0) or (ergodic_flag)):
		alpha_optimization_num = 0
		if ergodic_flag:
			for i in range(data_struct.m):
				alpha_optimization_num += self.inner_smo(i, data_struct)
				print("遍歷所有樣本資料：第%d次迭代，樣本為：%d，alpha優化的次數：%d" % (iter_num, i, alpha_optimization_num))
			iter_num += 1
		else:
			no_zero_index = np.nonzero((data_struct.alphas.A > 0) * (data_struct.alphas.A < C))[0]
			for i in no_zero_index:
				alpha_optimization_num += self.inner_smo(i, data_struct)
				print("非邊界遍歷樣本資料：第%d次迭代，樣本為：%d，alpha優化的次數：%d" % (iter_num, i, alpha_optimization_num))
			iter_num += 1
		if ergodic_flag: ergodic_flag = False
		elif alpha_optimization_num == 0: ergodic_flag = True
		print("迭代次數：%d" % iter_num)
	return data_struct.b, data_struct.alphas
  • 第二個\(\alpha_j\)的選取

對於第二個\(\alpha_j\)，我們不妨先分析一下前篇文章中SMO演算法最終優化之後的\(\alpha_2^{new}\)：

\[\alpha_2^{new}=\frac{y_2(E_1-E_2)}{\eta}+\alpha_2^{old} \]

我們可以發現\(\alpha_2\)主要是靠更新迭代來進行優化，而\(\alpha_2^{old}\)是已知的，我們沒有選擇的權利，但是\(\frac{y_2(E_1-E_2)}{\eta}\)這一部分的值我們是可控的，也就是說我們要選擇合適的\(\alpha_j\)來使得後一部分的值儘可能的大，從而達到快速修改\(\alpha\)向量的目的，才能更快速的實現訓練飽和。

簡單來說就是\(\alpha_2^{new}\)的變化依賴於\(|E_1-E_2|\)，當該絕對值越大，\(\alpha_2\)的變化也就越大。也就是說，\(E_1\)為正的時候，我們的要選擇儘可能小的\(E_2\)，\(E_1\)為負的時候，我們要選擇儘可能大的\(E_2\)。根據這個思想我們就能讓這對“小情侶”完美的配對了，美滋滋~~

對於第一個\(\alpha_i\)的選取，核心程式碼思想如下所示，也就是我們的內層迴圈的所需要選取\(\alpha_j\)內容：

def select_appropriate_j(self, i, data_struct, E_i):
	max_k, max_delta_E, E_j = -1, 0, 0
	data_struct.E_cache[i] = [1, E_i]
	valid_E_cache_list = np.nonzero(data_struct.E_cache[:, 0].A)[0]
	if (len(valid_E_cache_list) > 1):
		for k in valid_E_cache_list:
			if k == i: continue
			E_k = self.calc_E(data_struct.alphas, data_struct.y_label, data_struct.x_data, data_struct.b, k)
			delta_E = abs(E_i - E_k)
			if (delta_E > max_delta_E): max_k, max_delta_E, E_j = k, delta_E, E_k
		return max_k, E_j
	else: 
		j = self.random_select_alpha_j(i, data_struct.m)
		E_j = self.calc_E(data_struct.alphas, data_struct.y_label, data_struct.x_data, data_struct.b, j)
	return j, E_j

為了方便我們使用有關資料集和模型的一些公共資源，以及方便對它們進行操作，我們需要單獨封裝一個資料結構（當然了，不封裝也沒什麼問題），該資料結構有關屬性解釋如下：

  1. x_data：etablish_data隨機生成資料集中的屬性矩陣
  2. y_label：etablish_data隨機生成資料集中的標籤
  3. C：懲罰引數
  4. toler：容錯率
  5. m：資料樣本數
  6. alphas：SMO演算法所需要訓練的\(\alpha\)向量
  7. b：SMO演算法所需要訓練的\(b\)引數
  8. E_cache：用於儲存誤差，第一列為有效標誌位，第二列為樣本索引對應的誤差E

此外，為了提高程式碼的擴充套件性和靈活性，還單獨抽離了一個方法update_E_k，主要用於更新data_struct物件中的E_cache屬性：

def update_E_k(self, data_struct, k):
	E_k = self.calc_E(data_struct.alphas, data_struct.y_label, data_struct.x_data, data_struct.b, k)                                        #計算Ek
	data_struct.E_cache[k] = [1,E_k]

完整程式碼：

import numpy as np
import pylab as pl
from matplotlib import pyplot as plt

class DataStruct:
	def __init__(self, x_data, y_label, C, toler):
		self.x_data = x_data
		self.y_label = y_label
		self.C = C
		self.toler = toler
		self.m = x_data.shape[0]
		self.alphas = np.mat(np.zeros((self.m, 1)))
		self.b = 0
		self.E_cache = np.mat(np.zeros((self.m, 2)))

class OptimizeLinearSVM:
	def __init__(self):
		pass

	"""
	    Author: Taoye
	    微信公眾號: 玩世不恭的Coder
	    Explain: 用於生成訓練資料集
	    Parameters:
	        data_number: 樣本資料數目
	    Return:
	        x_data: 資料樣本的屬性矩陣
	        y_label: 樣本屬性所對應的標籤
	"""
	def etablish_data(self, data_number):
		np.random.seed(38)
		x_data = np.concatenate((np.add(np.random.randn(data_number, 2), [3, 3]),       
	                             np.subtract(np.random.randn(data_number, 2), [3, 3])),
	                             axis = 0)      # random隨機生成資料，+ -3達到不同類別資料分隔的目的 
		temp_data = np.zeros([data_number])
		temp_data.fill(-1)
		y_label = np.concatenate((temp_data, np.ones([data_number])), axis = 0)
		return x_data, y_label

	"""
	    Author: Taoye
	    微信公眾號: 玩世不恭的Coder
	    Explain: 隨機選取alpha_j
	    Parameters:
	        alpha_i_index: 第一個alpha的索引
	        alpha_number: alpha總數目
	    Return:
	        alpha_j_index: 第二個alpha的索引
	"""
	def random_select_alpha_j(self, alpha_i_index, alpha_number):
		alpha_j_index = alpha_i_index
		while alpha_j_index == alpha_i_index:
			alpha_j_index = np.random.randint(0, alpha_number)
		return alpha_j_index

	"""
	    Author: Taoye
	    微信公眾號: 玩世不恭的Coder
	    Explain: 使得alpha_j在[L, R]區間之內
	    Parameters:
	        alpha_j: 原始alpha_j
	        L: 左邊界值
	        R: 右邊界值
	    Return:
	        L,R,alpha_j: 修改之後的alpha_j
	"""
	def modify_alpha(self, alpha_j, L, R):
		if alpha_j < L: return L
		if alpha_j > R: return R
		return alpha_j

	"""
	    Author: Taoye
	    微信公眾號: 玩世不恭的Coder
	    Explain: 計算誤差並返回
	"""
	def calc_E(self, alphas, y_label, x_data, b, i):
		f_x_i = float(np.dot(np.multiply(alphas, y_label).T, x_data * x_data[i, :].T)) + b
		return f_x_i - float(y_label[i])

	"""
	    Author: Taoye
	    微信公眾號: 玩世不恭的Coder
	    Explain: 計算eta並返回
	"""
	def calc_eta(self, x_data, i, j):
		eta = 2.0 * x_data[i, :] * x_data[j, :].T \
	            - x_data[i, :] * x_data[i, :].T \
	            - x_data[j, :] * x_data[j,:].T
		return eta

	"""
	    Author: Taoye
	    微信公眾號: 玩世不恭的Coder
	    Explain: 計算b1, b2並返回
	"""
	def calc_b(self, b, x_data, y_label, alphas, alpha_i_old, alpha_j_old, E_i, E_j, i, j):
		b1 = b - E_i \
	         - y_label[i] * (alphas[i] - alpha_i_old) * x_data[i, :] * x_data[i, :].T \
	         - y_label[j] * (alphas[j] - alpha_j_old) * x_data[i, :] * x_data[j, :].T
		b2 = b - E_j \
	         - y_label[i] * (alphas[i] - alpha_i_old) * x_data[i, :] * x_data[j, :].T \
	         - y_label[j] * (alphas[j] - alpha_j_old) * x_data[j, :] * x_data[j, :].T
		return b1, b2

	def select_appropriate_j(self, i, data_struct, E_i):
		max_k, max_delta_E, E_j = -1, 0, 0
		data_struct.E_cache[i] = [1, E_i]
		valid_E_cache_list = np.nonzero(data_struct.E_cache[:, 0].A)[0]
		if (len(valid_E_cache_list) > 1):
			for k in valid_E_cache_list:
				if k == i: continue
				E_k = self.calc_E(data_struct.alphas, data_struct.y_label, data_struct.x_data, data_struct.b, k)
				delta_E = abs(E_i - E_k)
				if (delta_E > max_delta_E): max_k, max_delta_E, E_j = k, delta_E, E_k
			return max_k, E_j
		else: 
			j = self.random_select_alpha_j(i, data_struct.m)
			E_j = self.calc_E(data_struct.alphas, data_struct.y_label, data_struct.x_data, data_struct.b, j)
		return j, E_j

	def update_E_k(self, data_struct, k):
		E_k = self.calc_E(data_struct.alphas, data_struct.y_label, data_struct.x_data, data_struct.b, k)                                        #計算Ek
		data_struct.E_cache[k] = [1,E_k]
	
	"""
	    Author: Taoye
	    微信公眾號: 玩世不恭的Coder
	    Explain: smo內層
	"""
	def inner_smo(self, i, data_strcut):
		E_i = self.calc_E(data_strcut.alphas, data_struct.y_label, data_struct.x_data, data_struct.b, i)      # 呼叫calc_E方法計算樣本i的誤差
		if ((data_struct.y_label[i] * E_i < -data_struct.toler) and (data_struct.alphas[i] < data_struct.C)) or ((data_struct.y_label[i] * E_i > data_struct.toler) and (data_struct.alphas[i] > 0)):
			j, E_j = self.select_appropriate_j(i, data_strcut, E_i)              # 選取一個恰當的j
			alpha_i_old, alpha_j_old = data_struct.alphas[i].copy(), data_struct.alphas[j].copy()
			if (data_struct.y_label[i] != data_struct.y_label[j]):               # 確保alphas在[L, R]區間內
				L, R = max(0, data_struct.alphas[j] - data_struct.alphas[i]), min(data_struct.C, data_struct.C + data_struct.alphas[j] - data_struct.alphas[i])
			else:
				L, R = max(0, data_struct.alphas[j] + data_struct.alphas[i] - data_struct.C), min(data_struct.C, data_struct.alphas[j] + data_struct.alphas[i])
			if L == R: print("L==R"); return 0          # L==R時選取下一個樣本
			eta = self.calc_eta(data_struct.x_data, i, j)                  # 計算eta值
			if eta >= 0: print("eta>=0"); return 0
			data_struct.alphas[j] -= data_struct.y_label[j] * (E_i - E_j) / eta
			data_struct.alphas[j] = self.modify_alpha(data_struct.alphas[j], L, R)     # 修改alpha[j]
			self.update_E_k(data_strcut, j)
			if (abs(data_strcut.alphas[j] - alpha_j_old) < 0.000001): print("alpha_j修改太小了"); return 0
			data_struct.alphas[i] += data_strcut.y_label[j] * data_strcut.y_label[i] * (alpha_j_old - data_strcut.alphas[j])
			self.update_E_k(data_strcut, i)
			b1, b2= self.calc_b(data_struct.b, data_struct.x_data, data_struct.y_label, data_struct.alphas, alpha_i_old, alpha_j_old, E_i, E_j, i, j)    # 計算b值
			if (0 < data_struct.alphas[i]) and (data_struct.C > data_struct.alphas[i]): data_struct.b = b1
			elif (0 < data_struct.alphas[j]) and (data_struct.C > data_struct.alphas[j]): data_struct.b = b2
			else: data_struct.b = (b1 + b2)/2.0
			return 1
		else: return 0

	"""
	    Author：Taoye
	    微信公眾號：玩世不恭的Coder
	    Explain：外層迴圈，選取第一個合適alpha
	    Parameters：
	        x_data: 樣本屬性特徵矩陣
	        y_label: 屬性特徵對應的標籤
	        C：懲罰引數
	        toler：容錯率
	        max_iter：迭代次數
	    Return:
	        b: 決策面的引數b
	        alphas：獲取決策面引數w所需要的alphas
	"""
	def outer_smo(self, data_struct, x_data, y_label, C, toler, max_iter):
		iter_num, ergodic_flag, alpha_optimization_num = 0, True, 0
		while (iter_num < max_iter) and ((alpha_optimization_num > 0) or (ergodic_flag)):
			alpha_optimization_num = 0
			if ergodic_flag:
				for i in range(data_struct.m):
					alpha_optimization_num += self.inner_smo(i, data_struct)
					print("遍歷所有樣本資料：第%d次迭代，樣本為：%d，alpha優化的次數：%d" % (iter_num, i, alpha_optimization_num))
				iter_num += 1
			else:
				no_zero_index = np.nonzero((data_struct.alphas.A > 0) * (data_struct.alphas.A < C))[0]
				for i in no_zero_index:
					alpha_optimization_num += self.inner_smo(i, data_struct)
					print("非邊界遍歷樣本資料：第%d次迭代，樣本為：%d，alpha優化的次數：%d" % (iter_num, i, alpha_optimization_num))
				iter_num += 1
			if ergodic_flag: ergodic_flag = False
			elif alpha_optimization_num == 0: ergodic_flag = True
			print("迭代次數：%d" % iter_num)
		return data_struct.b, data_struct.alphas

	"""
	    Author: Taoye
	    微信公眾號: 玩世不恭的Coder
	    Explain: 根據公式計算出w權值向量
	    Parameters:
	        x_data: 樣本屬性特徵矩陣
	        y_label: 屬性特徵對應的標籤
	        alphas：linear_smo方法所返回的alphas向量
	    Return:
	        w: 決策面的引數w
	"""
	def calc_w(self, x_data, y_label, alphas):
		x_data, y_label, alphas = np.array(x_data), np.array(y_label), np.array(alphas)
		return np.dot((np.tile(y_label.reshape(1, -1).T, (1, 2)) * x_data).T, alphas).tolist()

	"""
	    Author: Taoye
	    微信公眾號: 玩世不恭的Coder
	    Explain: 繪製出分類結果
	    Parameters:
	        x_data: 樣本屬性特徵矩陣
	        y_label: 屬性特徵對應的標籤
	        w：決策面的w引數
	        b：決策面的引數b
	"""
	def plot_result(self, x_data, y_label, w, b):
		data_number, _ = x_data.shape; middle = int(data_number / 2)
		plt.scatter(x_data[:, 0], x_data[:, 1], c = y_label, cmap = pl.cm.Paired)
		x1, x2 = np.max(x_data), np.min(x_data)
		w1, w2 = w[0][0], w[1][0]
		y1, y2 = (-b - w1 * x1) / w2, (-b - w1 * x2) / w2
		plt.plot([float(x1), float(x2)], [float(y1), float(y2)])    # 繪製決策面
		for index, alpha in enumerate(alphas):
			if alpha > 0:
				b_temp = - w1 * x_data[index][0] - w2 * x_data[index][1]
				y1_temp, y2_temp = (-b_temp - w1 * x1) / w2, (-b_temp - w1 * x2) / w2
				plt.plot([float(x1), float(x2)], [float(y1_temp), float(y2_temp)], "k--")    # 繪製支援向量
				plt.scatter(x_data[index][0], x_data[index][1], s=150, c='none', alpha=0.7, linewidth=2, edgecolor='red')   # 圈出支援向量
		plt.show()

if __name__ == '__main__':
	optimize_linear_svm = OptimizeLinearSVM()
	x_data, y_label = optimize_linear_svm.etablish_data(50)
	data_struct = DataStruct(np.mat(x_data), np.mat(y_label).T, 0.8, 0.00001)
	b, alphas = optimize_linear_svm.outer_smo(data_struct, x_data, y_label, data_struct.C, data_struct.toler, 10)
	w = optimize_linear_svm.calc_w(x_data, y_label, alphas)
	optimize_linear_svm.plot_result(x_data, y_label, w, b)

優化後的分類結果：

這期的內容沒那麼多，主要是優化了一下上期內容中的SMO演算法，從而在一定程度上提高模型的訓練效率，由於是手動實現該線性SVM演算法，所以模型可能達不到那些框架內建的效能，有興趣的讀者可自行慢慢優化。

關於線性SVM應該就暫時寫到這裡了，後期的話應該會更新非線性相關的內容。

我是Taoye，愛專研，愛分享，熱衷於各種技術，學習之餘喜歡下象棋、聽音樂、聊動漫，希望藉此一畝三分地記錄自己的成長過程以及生活點滴，也希望能結實更多志同道合的圈內朋友，更多內容歡迎來訪微信公主號：玩世不恭的Coder

參考資料：

[1] 《機器學習實戰》：Peter Harrington 人民郵電出版社
[2] 《統計學習方法》：李航 第二版 清華大學出版社

推薦閱讀

《Machine Learning in Action》—— 剖析支援向量機，單手狂撕線性SVM
print( "Hello，NumPy！" )
幹啥啥不行，吃飯第一名
Taoye滲透到一家黑平臺總部，背後的真相細思極恐
《大話資料庫》-SQL語句執行時，底層究竟做了什麼小動作？
那些年，我們玩過的Git，真香
基於Ubuntu+Python+Tensorflow+Jupyter notebook搭建深度學習環境
網路爬蟲之頁面花式解析
手握手帶你瞭解Docker容器技術
一文詳解Hexo+Github小白建站
​開啟ElasticSearch、kibana、logstash的正確方式

相關文章