机器学习(4)-支持向量机的理解与代码实现（下）

这一篇接着上一篇机器学习(4)-支持向量机的理解与代码实现（上）继续介绍SMO算法以及对应的代码实现

序列最小最优化算法

上一篇我们讲到将原问题的优化问题转化为求解对偶问题的最优解，最终转化为求解如下约束规划问题：
$$\begin{gathered} \underset{\alpha }{\min }L(w,b,\alpha )=\frac{1}{2}\sum _{i=1}^N\sum _{j=1}^N{\alpha }_i{\alpha }_j{y}_i{y}_{j}K(x_i,x_j)-\sum _{i=1}^N{\alpha }_i \\ s.t.\sum _{i=1}^N{\alpha }_i{y}_i=0 \\ 0\le {\alpha }_i\le C,i=1,2,...,N \end{gathered}$$
上一篇是一个最大优化问题，这里我们对优化函数取一个负号，转化为最小优化问题。在实际运用中，随着训练数据的不断增加，上述优化问题需要求解的参数也等比例的增加，同时求解多个参数的算法将变得异常的非常和抵消，为了解决这个问题，Platt在1998年提出了SMO算法。基本思想如下：每次从要求解的参数中随机选取两个，固定其他参数，将原问题转为二元优化问题，通过零求导的梯度等于零，不断的更新参数直至所有的参数都满足KKT条件从而求得最优解。
下面详细讲解一下优化过程，不失一般性，假设我们选择的两个变量是${\alpha }_1$和${\alpha }_2$，去掉优化问题中与${\alpha }_1$和${\alpha }_2$无关的常数项(其他参数现在是固定的相当于已知量，在求导的过程中这些项的导数都为零)得到SMO需要优化的问题：
$$\begin{gathered} \underset{{\alpha }_1,{\alpha }_2}{\min }W({\alpha }_1,{\alpha }_2)=\frac{1}{2}{K}_{11}{\alpha }_1^2+\frac{1}{2}{K}_{22}{\alpha }_2^2+y_1{\alpha }_1\sum _{i=3}^N{y}_i{\alpha }_i{K}_{i1}+y_2{\alpha }_2\sum _{i=3}^N{y}_i{\alpha }_i{K}_{i2}+y_1{\alpha }_1{y}_2{\alpha }_2{K}_{12}-{\alpha }_1-{\alpha }_2 \\ s.t.{\alpha }_1{y}_1+{\alpha }_2{y}_2=-\sum _{i=3}^N{\alpha }_i{y}_i=\varsigma \\ 0\le {\alpha }_i\le C,i=1,2,...,N \end{gathered}$$
其中$K_{ij}=K(x_i,x_j),$对于等式${\alpha }_1{y}_1+{\alpha }_2{y}_2=-{\sum }_{i=3}^N{\alpha }_i{y}_i=\varsigma$左右同乘$y_1$可得${\alpha }_1=y_1(\varsigma -{\alpha }_2{y}_2)$,同时我们令$v_j={\sum }_{i=3}^N{y}_i{\alpha }_i{K}_{ij}$，将$v_1,v_2,{\alpha }_1$代入优化目标函数可以的到一个只和${\alpha }_2$有关的函数：
$$W({\alpha }_2)=\frac{1}{2}{K}_{11}(\varsigma -{\alpha }_2{y}_2{)}^2+\frac{1}{2}{K}_{22}{\alpha }_2^2+(\varsigma -{\alpha }_2{y}_2)v_1+y_2{\alpha }_2{v}_2+(\varsigma -{\alpha }_2{y}_2)y_2{\alpha }_2{K}_{12}-y_1(\varsigma -{\alpha }_2{y}_2)-{\alpha }_2$$
为了求得极小值我们对目标函数进行求导：
$$\frac{\partial W({\alpha }_2)}{\partial {\alpha }_2}={\alpha }_2(K_{11}+K_{22}-2K_{12})+\varsigma y_2(K_{12}-K_{11})+y_2(v_2-v_1)+y_1{y}_2-1$$
对于输入$x$模型的复测输出$f(x)={\sum }_{i=1}^N{y}_i{\alpha }_{i}K(x_i,x)+b$，我们用$f(x)$来表示$v_1$和$v_2$:
$$\begin{gathered} v_1=f(x_1)-b-y_1{\alpha }_1{K}_{11}-y_2{\alpha }_2{K}_{12} \\ {v}_2=f(x_2)-b-y_1{\alpha }_1{K}_{12}-y_2{\alpha }_2{K}_{22} \end{gathered}$$
前面我们已经知道${\alpha }_1=y_1(\varsigma -{\alpha }_2{y}_2)$,代入上式可以得到：
$$v_1-v_2=f(x_1)-f(x_2)+y_2{\alpha }_2(K_{11}+K_{22}-2K_{12})+\varsigma (K_{12}-K_{11})$$
因为预测值$f(x)$是根据更新前的参数得到的，我们令更新前的参数为${\alpha }^{old}$和$b^{old}$,带求解的更新参数为${\alpha }^{new}$和$b^{new}$,则$v_1-v_2$可以用旧的参数表示，同时我们令预测值和真实值的差值为$E_i=f(x_i)-y_i$将$v_1-v_2$代入梯度可得：
$$\begin{array}{rl}\frac{\partial W({\alpha }_2)}{\partial {\alpha }_2}& ={\alpha }_2(K_{11}+K_{22}-2K_{12})-{\alpha }_2^{old}(K_{11}+K_{22}-2K_{12})+y_2(y_1-f(x_1))-y_2(y_2-f(x_2))\\ & ={\alpha }_2(K_{11}+K_{22}-2K_{12})-{\alpha }_2^{old}(K_{11}+K_{22}-2K_{12})-y_2(E_1-E_2)\end{array}$$
令梯度等于零，可得${\alpha }_2$的更新值为：
$${\alpha }_2^{new}={\alpha }_2^{old}+\frac{y_2(E_1-E_2)}{\eta },\eta =K_{11}+K_{22}-2K_{12}$$
我们现在求出了${\alpha }_2$的更新值，但是这个是没有考虑参数的约束条件$0\le {\alpha }_i\le C$的情况下，下面我们根据约束条件对参数进行裁剪,我们重新命名${\alpha }_2^{new}$为${\alpha }_2^{new,unc}$,进行裁剪后的参数为${\alpha }_2^{new}$,根据约束条件和$y$的取值${\alpha }_2^{new}$的取值范围可以用如下两幅图表示：
左图是$y_1\ne y_2$的情况这个时候是一条斜率为一的线段，可以得到最小值和最大值分别为：
$$L=\max (0,{\alpha }_2^{old}-{\alpha }_1^{old}),H=min(C,C+{\alpha }_2^{old}-{\alpha }_1^{old})$$
当$y_1=y_2$时，如右图所示，此时最小和最大值分别为：
$$L=\max (0,{\alpha }_2^{old}+{\alpha }_1^{old}-C),H=\min (C,{\alpha }_2^{old}+{\alpha }_1^{old})$$
根据约束条件可以得到最终${\alpha }_2^{new}$的表达式如下：
$${\alpha }_2^{new}=\{\begin{array}{lr}L,& {\alpha }_2^{new,unc}<L\\ {\alpha }_2^{new,unc},& L\le {\alpha }_2^{new,unc}\le H\\ H,& {\alpha }_2^{new,unc}>H\end{array}$$
可以根据${\alpha }_2^{new}$可求得${\alpha }_1^{new}$：
$${\alpha }_1^{new}={\alpha }_1^{old}+y_1{y}_2({\alpha }_2^{old}-{\alpha }_2^{new})$$
得到了更新后的$\alpha$,就可以由此计算$b^{new}$了，当$0<{\alpha }_1^{new}<C$时，由KKT条件可知：
$$\sum _{i=1}^N{\alpha }_i{y}_i{K}_{i1}+b=y_1$$
其中除了${\alpha }_1$和${\alpha }_2$其他的都没有改变：
$$b_1^{new}=y_1-\sum _{i=3}^N{\alpha }_i{y}_i{K}_{i1}-{\alpha }_1^{new}{y}_1{K}_{11}-{\alpha }_2^{new}{y}_2{K}_{21}$$
右边等式的前两项可以用$E_1$表示为：
$$y_1-\sum _{i=3}^N{\alpha }_i{y}_i{K}_{i1}=-E_1+{\alpha }_1^{old}{y}_1{K}_{11}+{\alpha }_2^{old}{y}_2{K}_{21}+b^{old}$$
代入可得：
$$b_1^{new}=-E_1-y_1{K}_{11}({\alpha }_1^{new}-{\alpha }_1^{old})-y_2{K}_{21}({\alpha }_2^{new}-{\alpha }_2^{old})+b^{old}$$
同理当$0<{\alpha }_2^{new}<C$时：
$$b_2^{new}=-E2-y_1{K}_{12}({\alpha }_1^{new}-{\alpha }_1^{old})-y_2{K}_{22}({\alpha }_2^{new}-{\alpha }_2^{old})+b^{old}$$
当$0<{\alpha }_i^{new}<C,i=1,2$同时满足时，此时$b_1^{new}$和$b_2^{new}$是相等的，如果${\alpha }_1^{new}$和${\alpha }_2^{new}$等于C或者零时，$b_1^{new}$和$b_2^{new}$以及他们之间的任何值都是符合KKT条件的，这时我们可以取他们的平均值

代码实现

下面就是用代码实现SMO算法训练SVM的过程，我们可以把SVM抽象为一个类，然后主要分析下在求解${\alpha }_1^{new}$和${\alpha }_2^{new}$需要用到哪些值以及需要哪些功能，把这些功能分离出来用单独的函数实现，实际使用的时候再调用，就是把功能模块化。求解过程大致可以分为一下几个模块（1）任意两个样本之间的核函数值，因为会多次用到所以可以先一次性计算出来保存好（2）预测值和真实值之间的误差，而且在训练过程中还需要不断更新，我们可以用一个向量存储所有样本的误差，及时更新（3）如何挑选${\alpha }_1$和${\alpha }_2$,就是挑选哪两个参数用来更新，具体的可以看下李航老师书上介绍的方法（4）训练SVM模型（5）对测试集进行预测，大致可以分为上面五个模块，下面是具体的代码实现：

class SVM(Model):
	"""docstring for SVM"""
	def __init__(self, C=200, trainfile='../mnist_train.csv', sigma=10):
		super(SVM, self).__init__(trainfile)
		self.C = C
		self.sigma = sigma
		self.n, self.m = self.data.shape #(samplenum, featurenum)
		self.alpha = np.zeros(self.n)
		# self.alpha = np.random.random(self.n)
		self.b = 0
		self.E = np.zeros(self.n)
		self.g = np.zeros(self.n)
		self.data, self.label = self.preprocess(self.data, self.label)
		self.kernel = self.calGaussianKernel(self.data, self.data)
		# print(self.kernel)
		self.cal_E()

	def calGaussianKernel(self, data1, data2):
		n1 = data1.shape[0]
		n2 = data2.shape[0]
		data1_1 = np.sum(data1 * data1, axis=-1, keepdims=True) # (data1·data2) inner product
		data2_2 = np.sum(data2 * data2, axis=-1, keepdims=True) # (data2·data2)
		data1_2 = np.dot(data1, data2.T)
		data = np.tile(data1_1,(1, n2)) + np.tile(data2_2.T, (n1, 1)) - 2*data1_2
		data = np.exp(-data / (2 * self.sigma**2))
		return data

	def cal_E(self):
		self.g = np.dot(self.alpha * self.label, self.kernel) + self.b
		self.E = self.g - self.label

	def IsSatisfyKKT(self, i):
		if self.alpha[i] == 0 and self.label[i]*self.g[i] >= 1:
			return True
		if 0 < self.alpha[i] < self.C and self.label[i]*self.g[i] == 1:
			return True
		if self.alpha[i] == self.C and self.label[i]*self.g[i] <=1:
			return True
		return False

	def getAlpha2(self, i):
		maxdiff = 0
		index = i
		# print(self.E)
		for j in range(self.n):
			if abs(self.E[i] - self.E[j]) > maxdiff:
				maxdiff = abs(self.E[i] - self.E[j])
				index = j
		while(index == i):
			index = random.randint(0, self.n-1)
		return index

	@calTime
	def train(self, maxiter=50):
		it = 0
		stopiter = False
		while(it < maxiter and not stopiter):
			stopiter = True
			it += 1
			# print(self.alpha)
			for i in range(self.n):
				if not self.IsSatisfyKKT(i):
					j = self.getAlpha2(i)
					# print('i is {}, j is {}'.format(i, j))
					k11 = self.kernel[i][i]
					k22 = self.kernel[j][j]
					k12 = self.kernel[i][j]
					E1 = self.E[i]
					E2 = self.E[j]
					y1 = self.label[i]
					y2 = self.label[j]
					alpha1_old = self.alpha[i]
					alpha2_old = self.alpha[j]
					b_old = self.b
					# print('alpha2_old is {}, alpha1_old is {}'.format(alpha2_old, alpha1_old))
					eta = k11 + k22 - 2 * k12
					# print('k11 k12 k22',k11,k12,k22)
					alpha2_new = alpha2_old + y2 * (E1 - E2) / eta
					# print(y2 * (E1 - E2) / eta)
					# print('y1 is {}, y2 is {}'.format(y1, y2))
					if y1 == y2:
						L = max(0, alpha2_old + alpha1_old - self.C)
						H = min(self.C, alpha2_old + alpha1_old)
					else:
						L = max(0, alpha2_old - alpha1_old)
						H = min(self.C, self.C + alpha2_old - alpha1_old)
					if L == H:
						continue
					if alpha2_new < L:
						alpha2_new = L
					elif alpha2_new > H:
						alpha2_new = H
					# print('alpha2_new after cutting is {}'.format(alpha2_new - alpha2_old))
					alpha1_new = alpha1_old + y1*y2*(alpha2_old - alpha2_new)

					b1_new = -E1 - y1*k11*(alpha1_new - alpha1_old) - y2*k12*(alpha2_new - alpha2_old) + b_old
					b2_new = -E2 - y2*k22*(alpha2_new - alpha2_old) - y1*k12*(alpha1_new - alpha1_old) + b_old

					if 0 < alpha1_new < self.C:
						b_new = b1_new
					elif 0 < alpha2_new < self.C:
						b_new = b2_new
					else:
						b_new = (b1_new + b2_new) / 2
					# print('alpha2_new is {}, alpha1_new is {}, bnew is {}'.format(alpha2_new, alpha1_new, b_new))
					self.b = b_new
					self.alpha[i] = alpha1_new
					self.alpha[j] = alpha2_new
					self.cal_E()
					if abs(alpha2_new - alpha2_old) >= 0.00001:
						stopiter = False
						# print(self.alpha)
			print('iterate {}times'.format(it))



	def preprocess(self, data, label):
		label[label>0] = 1
		label[label==0] = -1
		
		
		#不能直接 a /= 255,会报错 TypeError: No loop matching the specified signature and casting was found for ufunc true_divide
		#或者可以先将data转化为float
		data = data / 255 
		return data, label

	@calTime
	def predict(self, data):
		kernel = self.calGaussianKernel(self.data, data)
		g = np.dot(self.alpha * self.label, kernel)
		# print(g)
		return np.sign(g)

---

顺便记录一下自己踩过的坑，一开始模型预测准确率一直为零，无论我怎么对数据做预处理调整模型的参数（其实这个时候我就应该反思了，我怎么调整参数输出都没有改变，至少说明除了模型代码肯定还有其他地方有问题），以及检查这个算法的代码都没有任何改变。最后打印我预测值才发现输出是None，但是我的预测函数是由返回值的，然后我就意识到是装饰器的原因了，我一开始写的装饰器没有返回值，所以导致我的predict一直是None，这是我第一次自己写代码用到装饰器，果然你以为你会了和你实际会了还有很远的距离，完整的代码后期我会整理放到github上,有错误的地方欢迎大家指正~