虽然SVM本身算法理论，水比较深，很难懂
但是基本原理却非常直观易懂，就是找到与训练集中支持向量有最大间隔的超平面
形式化的描述：
其中需要满足m个约束条件，m为数据集大小，即数据集中的每个数据点function margin都是>=1，因为之前假设所有支持向量，即离超平面最近的点，的function margin为1
对于这种有约束条件的最优化问题，用拉格朗日定理，于是得到如下的形式，
现在我们的目的就是求出最优化的m个拉格朗日算子，因为通过他们我们可以间接的算出w和b，从而得到最优超平面
考虑到某些离群值，加入松弛变量（slack variables），软间隔的版本为，
可以看到只是拉格朗日算子α的取值范围发生了变化，并且需要满足如下的KKT条件，
对于每个训练数据点都有一个对应的α
α对于绝大部分非支持向量而言，都是等于0的，这个由之前KKT条件可以推出
α等于C表示为离群值，个别function margin小于1的异常值
支持向量的α应该在0到C之间
所以这组KKT条件很容易理解，前面假设支持向量的function margin，即wx+b，为1
下面就是如何解这个问题？
之前大家都是用二次规划求解工具来求解，总之，这个计算量是非常大的
后来Platt提出SMO算法来优化这个求解过程，看过Andrew讲义就知道，这个算法是坐标上升算法的变种
下面就具体看看怎样实现SMO算法，关于算法的实现Andrew的讲义已经不够用，需要看Platt的论文
或参考这个blog，写的比较清楚，
SMO的简化版本
因为对于SMO算法，为了保持KKT条件，每次需要选取一对，即两个，α来进行优化
出于简单考虑，先随机的来选取α，后面再考虑启发式的来选取
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
dataMatrix = mat(dataMatIn); labelMat = mat(classLabels).transpose()
b = 0; m,n = shape(dataMatrix) #初始化b为0
alphas = mat(zeros((m,1))) #初始化所有α为0
iter = 0
while (iter < maxIter):
alphaPairsChanged = 0
for i in range(m): #遍历所有的训练集
fXi = float(multiply(alphas,labelMat).T * (dataMatrix*dataMatrix[i,:].T)) + b #1.计算wx+b
Ei = fXi - float(labelMat[i]) #和真实值比，计算误差
if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
j = selectJrand(i,m) #随机选择第二个α
fXj = float(multiply(alphas,labelMat).T * (dataMatrix*dataMatrix[j,:].T)) + b
Ej = fXj - float(labelMat[j])
alphaIold = alphas[i].copy(); #存下变化前两个α的值
alphaJold = alphas[j].copy();
if (labelMat[i] != labelMat[j]): #2.算出H和L，即α的边界值
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 = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - \ #3.计算新的α
dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
if eta >= 0: print "eta>=0"; continue
alphas[j] -= labelMat[j]*(Ei - Ej)/eta #更新αj
alphas[j] = clipAlpha(alphas[j],H,L) #和边界比较，不能超过边界
if (abs(alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; continue
alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j]) #根据αj的更新，反过来更新αi
b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*\ #4.更新b
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)
if (alphaPairsChanged == 0): iter += 1 #记录α没有变化的次数，如果多次α没有变化，说明已经收敛
else: iter = 0
print "iteration number: %d" % iter
return b,alphas
1. 计算wx+b
由于w是可以用α间接算出的，公式如下，
2. 计算α2的取值范围，即边界值，H和L，参考上面的blog
由约束条件，可以得到
等式右边是常数，而y不是1，就是-1，所以α1和α2是上图的线性关系（y1=1，y2=-1的case）
图上画出线不同的两种边界值的case，可以比较直观的看出，
对于y1=1，y2=1的case同理差不多
3. 计算新的α
从讲义中，只知道算出下面式子的极值，就可以得到新的α
但没有给出具体解的过程，其实还是比较复杂的，参考上面的blog
得到的解是，
其中，
K表示内积
看到这个公式就不难理解上面的代码的计算过程了
SMO的完整版本
关键的变化就是α的选择，采用启发式的来选择，从而大大优化算法速度
首先如何选取第一个α？
因为我们是需要通过最优化W来更新α的，所以对于边界值，即α=0或C的case，α是不会改变的，因为边界值的改变一定不会让分类效果变得更好
所以我们希望可以不断的优化支持向量的α，因为这个才是真正有效的优化，而不需要在边界α上浪费时间，这样可以大大提高算法效率
但问题是，所有α的初始值都是0，所以直接找非0的α是不可行的，所以这里的方法是，
先遍历所有α进行更新，这样会产生一些非0的α，然后对于这些α进行不断优化，直到收敛
然后再遍历所有α进行更新，也许会得到一些新的非0的α。。。。。
代码如下，
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=(‘lin‘, 0)):
oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler)
iter = 0
entireSet = True #第一遍是全遍历
alphaPairsChanged = 0
while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
alphaPairsChanged = 0
if entireSet:
for i in range(oS.m):
alphaPairsChanged += innerL(i,oS) #调用内循环去更新α
print "fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)
iter += 1
else:
nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
for i in nonBoundIs:
alphaPairsChanged += innerL(i,oS)
print "non-bound, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)
iter += 1
if entireSet: entireSet = False #一次全遍历后，切换成选择遍历
elif (alphaPairsChanged == 0): entireSet = True #当前α不改变了，再全遍历
return oS.b,oS.alphas
上面的blog中，全遍历应该选取那些违反KKT条件的α进行优化，这样应该更合理一些，那些符合KTT条件的没必要做优化
有空可以把这个实现改一下，看看会不会更快一些
如何选择第二个α？
相对于之前的random选取，改进为，选择最大误差步长的，即Ei-Ej最大的
即当
为正时选择负的绝对值最大的
，反之，选择正值最大的
这个直观上看，是有道理的，
因为另外一个α，需要根据当前α的改动，做相反的改动，以保证所有α的和不变
具体的代码参考原书，因为和简单版的相似
只是要缓存和更新所有的E值，并在选择第二个α时，做下计算选出Ei-Ej最大的
对非线性的复杂数据集使用kernels核函数
基础的SVM只能用于线性可分的数据集，如果不可分怎么办？
上面就是典型的例子，通常解决的方法，就是将低维数据转换为高维数据
比如上面的数据在二维看是线性不可分的，但是如果到了3维，可能蓝色的点和红色的点在高度上有些差异，这样就可以用一个平面，线性可分了
所以当数据非线性可分，当把数据的维度升高时，就有可能是线性可分的
对于SVM的对偶形式而言，只会计算内积形式，参考讲义
而我们如果把内积替换成核函数，就可以达成将低维数据到高维数据的转换，为什么？
因为所有的核函数，一定是等价于两个高维向量的内积，你不一定可以找出这个高维向量是什么，因为可能是无限维的向量，你也不需要care
所以用核函数替换当前低维向量的内积，就相当于用高维向量的内积替换低维向量的内积
下面看看具体的代码实现
这里实现的是最为流行的radial bias function的高斯版本，高斯核，对应于无限维的向量内积
对于m个数据点的训练集，我们要用高斯核函数，算出每对数据点间的核函数值，即m×m的结果矩阵
这样做训练或预测时，当需要计算内积时，直接查表即可
self.K = mat(zeros((self.m,self.m))) #初始化m*m矩阵
for i in range(self.m): #对于一个数据点，算出和其他所有数据点之间的核函数值，填入表内
self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)
下面看下kernelTrans
def kernelTrans(X, A, kTup):
m,n = shape(X)
K = mat(zeros((m,1)))
if kTup[0]==‘lin‘: K = X * A.T #线性核
elif kTup[0]==‘rbf‘: #高斯核
for j in range(m):
deltaRow = X[j,:] - A
K[j] = deltaRow*deltaRow.T
K = exp(K /(-1*kTup[1]**2))
else: raise NameError(‘Houston We Have a Problem -- That Kernel is not recognized‘)
return K
其中kTup[0]表示核类型，这里给出两种核的实现
kTup[1]是计算核需要的参数，对于高斯核就是a
剩下使用时，就把所有需要计算内积的地方，都替换成查表即可
Machine Learning in Action -- Support Vector Machines