diff --git a/docs/6.支持向量机.md b/docs/6.支持向量机.md index 96c4e29b..11da7481 100644 --- a/docs/6.支持向量机.md +++ b/docs/6.支持向量机.md @@ -144,3 +144,11 @@ SMO伪代码大致如下： 如果两个向量都不能被优化，退出内循环 如果所有向量都没被优化，增加迭代数目，继续下一次循环 ``` + +> SVM简化版：应用简化版SMO算法处理小规模数据集 + +代码可参考 svm-simple.py + +> SVM完整版：使用完整 Platt SMO算法加速优化 + +* 优化点：选择alpha的方式不同。 \ No newline at end of file diff --git a/src/python/6.SVM/svmMLiA.py b/src/python/6.SVM/svm-complete.py similarity index 91% rename from src/python/6.SVM/svmMLiA.py rename to src/python/6.SVM/svm-complete.py index d90b4a01..75252d5a 100644 --- a/src/python/6.SVM/svmMLiA.py +++ b/src/python/6.SVM/svm-complete.py @@ -8,8 +8,7 @@ Chapter 5 source file for Machine Learing in Action @author: Peter/geekidentity/片刻 """ from numpy import * -import pylab -from time import sleep +import matplotlib.pyplot as plt def loadDataSet(fileName): @@ -95,9 +94,9 @@ def smoSimple(dataMatIn, classLabels, C, toler, maxIter): # 记录alpha是否已经进行优化，每次循环时设为0，然后再对整个集合顺序遍历 alphaPairsChanged = 0 for i in range(m): - print 'alphas=', alphas - print 'labelMat=', labelMat - print 'multiply(alphas, labelMat)=', multiply(alphas, labelMat) + # print 'alphas=', alphas + # print 'labelMat=', labelMat + # print 'multiply(alphas, labelMat)=', multiply(alphas, labelMat) # 我们预测的类别 y = w^Tx[i]+b; 其中因为 w = Σ(1~n) a[n]lable[n]x[n] fXi = float(multiply(alphas, labelMat).T*(dataMatrix*dataMatrix[i, :].T)) + b # 预测结果与真实结果比对，计算误差Ei @@ -139,7 +138,7 @@ def smoSimple(dataMatIn, classLabels, C, toler, maxIter): alphas[j] -= labelMat[j]*(Ei - Ej)/eta # 并使用辅助函数，以及L和H对其进行调整 alphas[j] = clipAlpha(alphas[j], H, L) - # 检查alpha[j]是否有轻微的改变，如果是的话，就退出for循环。 + # 检查alpha[j]是否只是轻微的改变，如果是的话，就退出for循环。 if (abs(alphas[j] - alphaJold) < 0.00001): print("j not moving enough") continue @@ -165,11 +164,82 @@ def smoSimple(dataMatIn, classLabels, C, toler, maxIter): return b, alphas +def calcWs(alphas, dataArr, classLabels): + """ + 基于alpha计算w值 + Args: + alphas 拉格朗日乘子 + dataArr feature数据集 + classLabels 目标变量数据集 + + Returns: + wc 回归系数 + """ + X = mat(dataArr) + labelMat = mat(classLabels).transpose() + m, n = shape(X) + w = zeros((n, 1)) + for i in range(m): + w += multiply(alphas[i] * labelMat[i], X[i, :].T) + return w + + +def plotfig_SVM(xMat, yMat, ws, b, alphas): + """ + 参考地址： + http://blog.csdn.net/maoersong/article/details/24315633 + http://www.cnblogs.com/JustForCS/p/5283489.html + http://blog.csdn.net/kkxgx/article/details/6951959 + """ + + xMat = mat(xMat) + yMat = mat(yMat) + + # b原来是矩阵，先转为数组类型后其数组大小为（1,1），所以后面加[0]，变为(1,) + b = array(b)[0] + fig = plt.figure() + ax = fig.add_subplot(111) + + # 注意flatten的用法 + ax.scatter(xMat[:, 0].flatten().A[0], xMat[:, 1].flatten().A[0]) + + # x最大值，最小值根据原数据集dataArr[:, 0]的大小而定 + x = arange(-1.0, 10.0, 0.1) + + # 根据x.w + b = 0 得到，其式子展开为w0.x1 + w1.x2 + b = 0, x2就是y值 + y = (-b-ws[0, 0]*x)/ws[1, 0] + ax.plot(x, y) + for i in range(len(yMat)): + if yMat[i, 0] > 0: + ax.plot(xMat[i, 0], xMat[i, 1], 'cx') + else: + ax.plot(xMat[i, 0], xMat[i, 1], 'kp') + + # 找到支持向量，并在图中标红 + for i in range(100): + if alphas[i] > 0.0: + ax.plot(xMat[i, 0], xMat[i, 1], 'ro') + plt.show() + + if __name__ == "__main__": # 获取特征和目标变量 dataArr, labelArr = loadDataSet('input/6.SVM/testSet.txt') # print labelArr - smoSimple(dataArr, labelArr, 0.6, 0.001, 40) + + # b是常量值， alphas是拉格朗日乘子 + b, alphas = smoSimple(dataArr, labelArr, 0.6, 0.001, 40) + print '/n/n/n' + print 'b=', b + print 'alphas[alphas>0]=', alphas[alphas > 0] + print 'shape(alphas[alphas > 0])=', shape(alphas[alphas > 0]) + for i in range(100): + if alphas[i] > 0: + print dataArr[i], labelArr[i] + # 画图 + ws = calcWs(alphas, dataArr, labelArr) + plotfig_SVM(dataArr, labelArr, ws, b, alphas) + def kernelTrans(X, A, kTup): # calc the kernel or transform data to a higher dimensional space @@ -383,26 +453,6 @@ def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)): return oS.b, oS.alphas -def calcWs(alphas, dataArr, classLabels): - """ - 基于alpha计算w值 - Args: - alphas: - dataArr: - classLabels: - - Returns: - - """ - X = mat(dataArr) - labelMat = mat(classLabels).transpose() - m, n = shape(X) - w = zeros((n, 1)) - for i in range(m): - w += multiply(alphas[i] * labelMat[i], X[i, :].T) - return w - - def testRbf(k1=1.3): dataArr, labelArr = loadDataSet('testSetRBF.txt') b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) # C=200 important diff --git a/src/python/6.SVM/svm-simple.py b/src/python/6.SVM/svm-simple.py new file mode 100644 index 00000000..d59647c9 --- /dev/null +++ b/src/python/6.SVM/svm-simple.py @@ -0,0 +1,241 @@ +#!/usr/bin/python +# coding:utf8 + +""" +Created on Nov 4, 2010 +Update on 2017-03-21 +Chapter 5 source file for Machine Learing in Action +@author: Peter/geekidentity/片刻 +""" +from numpy import * +import matplotlib.pyplot as plt + + +def loadDataSet(fileName): + """ + 对文件进行逐行解析，从而得到第行的类标签和整个数据矩阵 + Args: + fileName 文件名 + Returns: + dataMat 数据矩阵 + labelMat 类标签 + """ + dataMat = [] + labelMat = [] + fr = open(fileName) + for line in fr.readlines(): + lineArr = line.strip().split('\t') + dataMat.append([float(lineArr[0]), float(lineArr[1])]) + labelMat.append(float(lineArr[2])) + return dataMat, labelMat + + +def selectJrand(i, m): + """ + 随机选择一个整数 + Args: + i 第一个alpha的下标 + m 所有alpha的数目 + Returns: + j 返回一个不为i的随机数，在0~m之间的整数值 + """ + j = i + while j == i: + j = int(random.uniform(0, m)) + return j + + +def clipAlpha(aj, H, L): + """clipAlpha(调整aj的值，使aj处于 L<=aj<=H) + Args: + aj 目标值 + H 最大值 + L 最小值 + Returns: + aj 目标值 + """ + if aj > H: + aj = H + if L > aj: + aj = L + return aj + + +def smoSimple(dataMatIn, classLabels, C, toler, maxIter): + """smoSimple + + Args: + dataMatIn 数据集 + classLabels 类别标签 + C 松弛变量(常量值)，允许有些数据点可以处于分隔面的错误一侧。 + 控制最大化间隔和保证大部分的函数间隔小于1.0这两个目标的权重。 + 可以通过调节该参数达到不同的结果。 + toler 容错率 + maxIter 退出前最大的循环次数 + Returns: + b 模型的常量值 + alphas 拉格朗日乘子 + """ + dataMatrix = mat(dataMatIn) + # 矩阵转制 和 .T 一样的功能 + labelMat = mat(classLabels).transpose() + m, n = shape(dataMatrix) + + # 初始化 b和alphas(alpha有点类似权重值。) + b = 0 + alphas = mat(zeros((m, 1))) + + # 没有任何alpha改变的情况下遍历数据的次数 + iter = 0 + while (iter < maxIter): + # w = calcWs(alphas, dataMatIn, classLabels) + # print("w:", w) + + # 记录alpha是否已经进行优化，每次循环时设为0，然后再对整个集合顺序遍历 + alphaPairsChanged = 0 + for i in range(m): + # print 'alphas=', alphas + # print 'labelMat=', labelMat + # print 'multiply(alphas, labelMat)=', multiply(alphas, labelMat) + # 我们预测的类别 y = w^Tx[i]+b; 其中因为 w = Σ(1~n) a[n]lable[n]x[n] + fXi = float(multiply(alphas, labelMat).T*(dataMatrix*dataMatrix[i, :].T)) + b + # 预测结果与真实结果比对，计算误差Ei + Ei = fXi - float(labelMat[i]) + + # 约束条件 (KKT条件是解决最优化问题的时用到的一种方法。我们这里提到的最优化问题通常是指对于给定的某一函数，求其在指定作用域上的全局最小值。) + # 0<=alphas[i]<=C，但由于0和C是边界值，我们无法进行优化，因为需要增加一个alphas和降低一个alphas。 + # 表示发生错误的概率：labelMat[i]*Ei 如果超出了 toler， 才需要优化。至于正负号，我们考虑绝对值就对了。 + if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)): + + # 如果满足优化的条件，我们就随机选取非i的一个点，进行优化比较 + j = selectJrand(i, m) + # 预测j的结果 + fXj = float(multiply(alphas, labelMat).T*(dataMatrix*dataMatrix[j, :].T)) + b + Ej = fXj - float(labelMat[j]) + alphaIold = alphas[i].copy() + alphaJold = alphas[j].copy() + + # L和H用于将alphas[j]调整到0-C之间。如果L==H，就不做任何改变，直接执行continue语句 + if (labelMat[i] != labelMat[j]): + L = max(0, alphas[j] - alphas[i]) + H = min(C, C + alphas[j] - alphas[i]) + else: + L = max(0, alphas[j] + alphas[i] - C) + H = min(C, alphas[j] + alphas[i]) + if L == H: + print("L==H") + continue + + # eta是alphas[j]的最优修改量，如果eta==0，需要退出for循环的当前迭代过程 + # 如果ETA为0，那么计算新的alphas[j]就比较麻烦了, 为什么呢？ 因为2个值一样。 + # 2ab <= a^2 + b^2 + eta = 2.0 * dataMatrix[i, :]*dataMatrix[j, :].T - dataMatrix[i, :]*dataMatrix[i,:].T - dataMatrix[j, :]*dataMatrix[j, :].T + if eta >= 0: + print("eta>=0") + continue + + # 计算出一个新的alphas[j]值 + alphas[j] -= labelMat[j]*(Ei - Ej)/eta + # 并使用辅助函数，以及L和H对其进行调整 + alphas[j] = clipAlpha(alphas[j], H, L) + # 检查alpha[j]是否只是轻微的改变，如果是的话，就退出for循环。 + if (abs(alphas[j] - alphaJold) < 0.00001): + print("j not moving enough") + continue + # 然后alphas[i]和alphas[j]同样进行改变，虽然改变的大小一样，但是改变的方向正好相反 + alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j]) + # 在对alpha[i], alpha[j] 进行优化之后，给这两个alpha值设置一个常数b。 + b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T + b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T + if (0 < alphas[i]) and (C > alphas[i]): + b = b1 + elif (0 < alphas[j]) and (C > alphas[j]): + b = b2 + else: + b = (b1 + b2)/2.0 + alphaPairsChanged += 1 + print("iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged)) + # 在for循环外，检查alpha值是否做了更新，如果在更新则将iter设为0后继续运行程序 + if (alphaPairsChanged == 0): + iter += 1 + else: + iter = 0 + print("iteration number: %d" % iter) + return b, alphas + + +def calcWs(alphas, dataArr, classLabels): + """ + 基于alpha计算w值 + Args: + alphas 拉格朗日乘子 + dataArr feature数据集 + classLabels 目标变量数据集 + + Returns: + wc 回归系数 + """ + X = mat(dataArr) + labelMat = mat(classLabels).transpose() + m, n = shape(X) + w = zeros((n, 1)) + for i in range(m): + w += multiply(alphas[i] * labelMat[i], X[i, :].T) + return w + + +def plotfig_SVM(xMat, yMat, ws, b, alphas): + """ + 参考地址： + http://blog.csdn.net/maoersong/article/details/24315633 + http://www.cnblogs.com/JustForCS/p/5283489.html + http://blog.csdn.net/kkxgx/article/details/6951959 + """ + + xMat = mat(xMat) + yMat = mat(yMat) + + # b原来是矩阵，先转为数组类型后其数组大小为（1,1），所以后面加[0]，变为(1,) + b = array(b)[0] + fig = plt.figure() + ax = fig.add_subplot(111) + + # 注意flatten的用法 + ax.scatter(xMat[:, 0].flatten().A[0], xMat[:, 1].flatten().A[0]) + + # x最大值，最小值根据原数据集dataArr[:, 0]的大小而定 + x = arange(-1.0, 10.0, 0.1) + + # 根据x.w + b = 0 得到，其式子展开为w0.x1 + w1.x2 + b = 0, x2就是y值 + y = (-b-ws[0, 0]*x)/ws[1, 0] + ax.plot(x, y) + for i in range(len(yMat)): + if yMat[i, 0] > 0: + ax.plot(xMat[i, 0], xMat[i, 1], 'cx') + else: + ax.plot(xMat[i, 0], xMat[i, 1], 'kp') + + # 找到支持向量，并在图中标红 + for i in range(100): + if alphas[i] > 0.0: + ax.plot(xMat[i, 0], xMat[i, 1], 'ro') + plt.show() + + +if __name__ == "__main__": + # 获取特征和目标变量 + dataArr, labelArr = loadDataSet('input/6.SVM/testSet.txt') + # print labelArr + + # b是常量值， alphas是拉格朗日乘子 + b, alphas = smoSimple(dataArr, labelArr, 0.6, 0.001, 40) + print '/n/n/n' + print 'b=', b + print 'alphas[alphas>0]=', alphas[alphas > 0] + print 'shape(alphas[alphas > 0])=', shape(alphas[alphas > 0]) + for i in range(100): + if alphas[i] > 0: + print dataArr[i], labelArr[i] + # 画图 + ws = calcWs(alphas, dataArr, labelArr) + plotfig_SVM(dataArr, labelArr, ws, b, alphas)