程序示例–基于 SMO 的 SVM 模型

在这里，我们会实现一个基于 SMO 的 SVM 模型，在其中，提供了简化版 SMO 和 完整版 SMO 的实现。

简化版 SMO：不使用启发式方法选择 (α(i),α(j))(α^{(i)},α^{(j)})(α(i),α(j)) ，始终在整个样本集上选择违反了 KKT 条件的 α(i)α^{(i)}α(i) ， α(j)α^{(j)}α(j) 则随机选择。
完整版 SMO：使用启发式方法选择 (α(i),α(j))(α^{(i)},α^{(j)})(α(i),α(j)) 。在整个训练集或者非边界样本中选择违反了 KKT 条件的 α(i)α^{(i)}α(i) 。并选择使得 ∣E(i)−E(j)∣|E^{(i)}−E^{(j)}|∣E(i)−E(j)∣ 达到最大的 alpha(j)alpha^{(j)}alpha(j) 。

核函数

模型支持了核函数，并提供了线性核函数 和 高斯核（又称之为径向基核函数 RBF）函数 使用：

# coding: utf8
# svm/smo.pydef linearKernel():"""线性核函数"""def calc(X, A):return X * A.Treturn calcdef rbfKernel(delta):"""rbf核函数"""gamma = 1.0 / (2 * delta**2)def calc(X, A):return np.mat(rbf_kernel(X, A, gamma=gamma))return calc

RBF 核的实现我们使用了性能更高的 sklearn.metrics.pairwise 提供的的 rbf_kernel

训练初始化

def getSmo(X, y, C, tol, maxIter, kernel=linearKernel()):"""SMOArgs:X 训练样本y 标签集C 正规化参数tol 容忍值maxIter 最大迭代次数K 所用核函数Returns:trainSimple 简化版训练算法train 完整版训练算法predict 预测函数"""m, n = X.shape# 存放核函数的转化结果K = kernel(X, X)# Cache存放预测误差，用以加快计算速度ECache = np.zeros((m,2))# ...

权值向量 www

我们知道，模型 f(x)f(x)f(x) 仅仅与支持向量有关，所以，权值 www 可以定义为：
w=∑i=1ma(i)y(i)(x(i))Tw=\sum_{i=1}^ma^{(i)}y^{(i)}(x^{(i)})^Tw=i=1∑ma(i)y(i)(x(i))T=∑i∈1a(i)y(i)(x(i))T=\sum_{i∈1}a^{(i)}y^{(i)}(x^{(i)})^T=i∈1∑a(i)y(i)(x(i))T

其中， SSS 表示了属于支持向量的样本坐标集合。

# coding: utf8
# svm/smo.pydef getSmo(X, y, C, tol, maxIter, kernel=linearKernel()):# ...def w(alphas, b, supportVectorsIndex, supportVectors):return (np.multiply(alphas[supportVectorsIndex], y[supportVectorsIndex]).T * supportVectors).T

预测误差 E(i)E^{(i)}E(i)

E(i)=f(x(i))−y(i)E^{(i)}=f(x^{(i)})-y^{(i)}E(i)=f(x(i))−y(i)=∑j=1ma(j)y(j)k(x(j),x(i))+b−y(i)=\sum_{j=1}^ma^{(j)}y^{(j)}k(x^{(j)},x^{(i)})+b-y^{(i)}=j=1∑ma(j)y(j)k(x(j),x(i))+b−y(i)

# coding: utf8
# svm/smo.pydef getSmo(X, y, C, tol, maxIter, kernel=linearKernel()):# ...def E(i, alphas, b):"""计算预测误差Args:i ialphas alphasb bReturns:E_i 第i个样本的预测误差"""FXi = float(np.multiply(alphas, y).T * K[:, i]) + bE = FXi - float(y[i])return E

预测函数

f(x)=∑i∈Sa(i)y(i)(x(i))T+bf(x)=\sum_{i∈S}a^{(i)}y^{(i)}(x^{(i)})^T+bf(x)=i∈S∑a(i)y(i)(x(i))T+b

# coding: utf8
# svm/smo.pydef getSmo(X, y, C, tol, maxIter, kernel=linearKernel()):# ...def predict(X, alphas, b, supportVectorsIndex, supportVectors):"""计算权值向量Args:X 预测矩阵alphas alphasb bsupportVectorsIndex 支持向量坐标supportVectors 支持向量Returns:predicts 预测结果"""Ks = kernel(supportVectors, X)predicts = (np.multiply(alphas[supportVectorsIndex], y[supportVectorsIndex]).T * Ks + b).Tpredicts = np.sign(predicts)return predicts

选择 α(i)α^{(i)}α(i)

# coding: utf8
# svm/smo.pydef getSmo(X, y, C, tol, maxIter, kernel=linearKernel()):# ...def select(i, alphas, b):"""alpha对选择"""Ei = E(i, alphas, b)# 选择违背KKT条件的，作为alpha2Ri = y[i] * Eiif (Ri < -tol and alphas[i] < C) or \(Ri > tol and alphas[i] > 0):# 选择第二个参数j = selectJRand(i)Ej = E(j, alphas, b)# j, Ej = selectJ(i, Ei, alphas, b)# get boundsif y[i] != y[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:return 0, alphas, bKii = K[i, i]Kjj = K[j, j]Kij = K[i, j]eta = 2.0 * Kij - Kii - Kjjif eta >= 0:return 0, alphas, biOld = alphas[i].copy()jOld = alphas[j].copy()alphas[j] = jOld - y[j] * (Ei - Ej) / etaif alphas[j] > H:alphas[j] = Helif alphas[j] < L:alphas[j] = Lif abs(alphas[j] - jOld) < tol:alphas[j] = jOldreturn 0, alphas, balphas[i] = iOld + y[i] * y[j] * (jOld - alphas[j])# update ECacheupdateE(i, alphas, b)updateE(j, alphas, b)# update bbINew = b - Ei - y[i] * (alphas[i] - iOld) * Kii - y[j] * \(alphas[j] - jOld) * KijbJNew = b - Ej - y[i] * (alphas[i] - iOld) * Kij - y[j] * \(alphas[j] - jOld) * Kjjif alphas[i] > 0 and alphas[i] < C:bNew = bINewelif alphas[j] > 0 and alphas[j] < C:bNew = bJNewelse:bNew = (bINew + bJNew) / 2return 1, alphas, belse:return 0, alphas, b

选择 α(j)α^{(j)}α(j)

# coding: utf8
# svm/smo.pydef getSmo(X, y, C, tol, maxIter, kernel=linearKernel()):# ...def selectJRand(i):""""""j = iwhile j == i:j = int(np.random.uniform(0, m))return jdef selectJ(i, Ei, alphas, b):"""选择权值 {% math %}\alpha^{(i)}{% endmath %}"""maxJ = 0; maxDist=0; Ej = 0ECache[i] = [1, Ei]validCaches = np.nonzero(ECache[:, 0])[0]if len(validCaches) > 1:for k in validCaches:if k==i: continueEk = E(k, alphas, b)dist = np.abs(abs(Ei-Ek))if maxDist < dist:Ej = EkmaxJ = kmaxDist = distreturn maxJ, Ejelse:### 随机选择j = selectJRand(i)Ej = E(j, alphas, b)return j, Ej

完整版 SMO 训练方法：

# coding: utf8
# svm/smo.pydef getSmo(X, y, C, tol, maxIter, kernel=linearKernel()):# ...def train():"""完整版训练算法Returns:alphas alphasw wb bsupportVectorsIndex 支持向量的坐标集supportVectors 支持向量iterCount 迭代次数"""numChanged = 0examineAll = TrueiterCount = 0alphas = np.mat(np.zeros((m, 1)))b = 0# 如果所有alpha都遵从 KKT 条件，则在整个训练集上迭代# 否则在处于边界内 (0, C) 的 alpha 中迭代while (numChanged > 0 or examineAll) and (iterCount < maxIter):numChanged = 0if examineAll:for i in range(m):changed, alphas, b = select(i, alphas, b)numChanged += changedelse:nonBoundIds = np.nonzero((alphas.A > 0) * (alphas.A < C))[0]for i in nonBoundIds:changed, alphas, b = select(i, alphas, b)numChanged += changediterCount += 1if examineAll:examineAll = Falseelif numChanged == 0:examineAll = TruesupportVectorsIndex = np.nonzero(alphas.A > 0)[0]supportVectors = np.mat(X[supportVectorsIndex])return alphas, w(alphas, b, supportVectorsIndex, supportVectors), b, \supportVectorsIndex, supportVectors, iterCount

简化版 SMO 训练方法：

# coding: utf8
# svm/smo.pydef getSmo(X, y, C, tol, maxIter, kernel=linearKernel()):# ...def trainSimple():"""简化版训练算法Returns:alphas alphasw wb bsupportVectorsIndex 支持向量的坐标集supportVectors 支持向量iterCount 迭代次数"""numChanged = 0iterCount = 0alphas = np.mat(np.zeros((m, 1)))b = 0L = 0H = 0while iterCount < maxIter:numChanged = 0for i in range(m):Ei = E(i, alphas, b)Ri = y[i] * Ei# 选择违背KKT条件的，作为alpha2if (Ri < -tol and alphas[i] < C) or \(Ri > tol and alphas[i] > 0):# 选择第二个参数j = selectJRand(i)Ej = E(j, alphas, b)# get boundsif y[i] != y[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:continueKii = K[i, i]Kjj = K[j, j]Kij = K[i, j]eta = 2.0 * Kij - Kii - Kjjif eta >= 0:continueiOld = alphas[i].copy();jOld = alphas[j].copy()alphas[j] = jOld - y[j] * (Ei - Ej) / etaif alphas[j] > H:alphas[j] = Helif alphas[j] < L:alphas[j] = Lif abs(alphas[j] - jOld) < tol:alphas[j] = jOldcontinuealphas[i] = iOld + y[i] * y[j] * (jOld - alphas[j])# update bbINew = b - Ei - y[i] * (alphas[i] - iOld) * Kii - y[j] * \(alphas[j] - jOld) * KijbJNew = b - Ej - y[i] * (alphas[i] - iOld) * Kij - y[j] * \(alphas[j] - jOld) * Kjjif alphas[i] > 0 and alphas[i] < C:b = bINewelif alphas[j] > 0 and alphas[j] < C:b = bJNewelse:b = (bINew + bJNew) / 2.0numChanged += 1if numChanged == 0:iterCount += 1else:iterCount = 0supportVectorsIndex = np.nonzero(alphas.A > 0)[0]supportVectors = np.mat(X[supportVectorsIndex])return alphas, w(alphas, b, supportVectorsIndex, supportVectors), b, \supportVectorsIndex, supportVectors, iterCount

完整代码

# coding: utf8
# svm/smo.pyimport numpy as npfrom sklearn.metrics.pairwise import rbf_kernel"""
svm模型
"""def linearKernel():"""线性核函数"""def calc(X, A):return X * A.Treturn calcdef rbfKernel(delta):"""rbf核函数"""gamma = 1.0 / (2 * delta**2)def calc(X, A):return np.mat(rbf_kernel(X, A, gamma=gamma))return calcdef getSmo(X, y, C, tol, maxIter, kernel=linearKernel()):"""SMOArgs:X 训练样本y 标签集C 正规化参数tol 容忍值maxIter 最大迭代次数K 所用核函数Returns:trainSimple 简化版训练算法train 完整版训练算法predict 预测函数"""# 存放核函数的转化结果K = kernel(X, X)# Cache存放预测误差，用以加快计算速度ECache = np.zeros((m,2))def predict(X, alphas, b, supportVectorsIndex, supportVectors):"""计算权值向量Args:X 预测矩阵alphas alphasb bsupportVectorsIndex 支持向量坐标集supportVectors 支持向量Returns:predicts 预测结果"""Ks = kernel(supportVectors, X)predicts = (np.multiply(alphas[supportVectorsIndex], y[supportVectorsIndex]).T * Ks + b).Tpredicts = np.sign(predicts)return predictsdef w(alphas, b, supportVectorsIndex, supportVectors):"""计算权值Args:alphas alphasb bsupportVectorsIndex 支持向量坐标supportVectors 支持向量Returns:w 权值向量"""return (np.multiply(alphas[supportVectorsIndex], y[supportVectorsIndex]).T * supportVectors).Tdef E(i, alphas, b):"""计算预测误差Args:i ialphas alphasb bReturns:E_i 第i个样本的预测误差"""FXi = float(np.multiply(alphas, y).T * K[:, i]) + bE = FXi - float(y[i])return Edef updateE(i, alphas, b):ECache[i] = [1, E(i, alphas, b)]def selectJRand(i):""""""j = iwhile j == i:j = int(np.random.uniform(0, m))return jdef selectJ(i, Ei, alphas, b):"""选择权值 {% math %}\alpha^{(i)}{% endmath %}"""maxJ = 0; maxDist=0; Ej = 0ECache[i] = [1, Ei]validCaches = np.nonzero(ECache[:, 0])[0]if len(validCaches) > 1:for k in validCaches:if k==i: continueEk = E(k, alphas, b)dist = np.abs(abs(Ei-Ek))if maxDist < dist:Ej = EkmaxJ = kmaxDist = distreturn maxJ, Ejelse:### 随机选择j = selectJRand(i)Ej = E(j, alphas, b)return j, Ejdef select(i, alphas, b):"""alpha对选择"""Ei = E(i, alphas, b)# 选择违背KKT条件的，作为alpha2Ri = y[i] * Eiif (Ri < -tol and alphas[i] < C) or \(Ri > tol and alphas[i] > 0):# 选择第二个参数j = selectJRand(i)Ej = E(j, alphas, b)# j, Ej = selectJ(i, Ei, alphas, b)# get boundsif y[i] != y[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:return 0, alphas, bKii = K[i, i]Kjj = K[j, j]Kij = K[i, j]eta = 2.0 * Kij - Kii - Kjjif eta >= 0:return 0, alphas, biOld = alphas[i].copy()jOld = alphas[j].copy()alphas[j] = jOld - y[j] * (Ei - Ej) / etaif alphas[j] > H:alphas[j] = Helif alphas[j] < L:alphas[j] = Lif abs(alphas[j] - jOld) < tol:alphas[j] = jOldreturn 0, alphas, balphas[i] = iOld + y[i] * y[j] * (jOld - alphas[j])# update ECacheupdateE(i, alphas, b)updateE(j, alphas, b)# update bbINew = b - Ei - y[i] * (alphas[i] - iOld) * Kii - y[j] * \(alphas[j] - jOld) * KijbJNew = b - Ej - y[i] * (alphas[i] - iOld) * Kij - y[j] * \(alphas[j] - jOld) * Kjjif alphas[i] > 0 and alphas[i] < C:bNew = bINewelif alphas[j] > 0 and alphas[j] < C:bNew = bJNewelse:bNew = (bINew + bJNew) / 2return 1, alphas, belse:return 0, alphas, bdef train():"""完整版训练算法Returns:alphas alphasw wb bsupportVectorsIndex 支持向量的坐标集supportVectors 支持向量iterCount 迭代次数"""numChanged = 0examineAll = TrueiterCount = 0alphas = np.mat(np.zeros((m, 1)))b = 0# 如果所有alpha都遵从 KKT 条件，则在整个训练集上迭代# 否则在处于边界内 (0, C) 的 alpha 中迭代while (numChanged > 0 or examineAll) and (iterCount < maxIter):numChanged = 0if examineAll:for i in range(m):changed, alphas, b = select(i, alphas, b)numChanged += changedelse:nonBoundIds = np.nonzero((alphas.A > 0) * (alphas.A < C))[0]for i in nonBoundIds:changed, alphas, b = select(i, alphas, b)numChanged += changediterCount += 1if examineAll:examineAll = Falseelif numChanged == 0:examineAll = TruesupportVectorsIndex = np.nonzero(alphas.A > 0)[0]supportVectors = np.mat(X[supportVectorsIndex])return alphas, w(alphas, b, supportVectorsIndex, supportVectors), b, \supportVectorsIndex, supportVectors, iterCountdef trainSimple():"""简化版训练算法Returns:alphas alphasw wb bsupportVectorsIndex 支持向量的坐标集supportVectors 支持向量iterCount 迭代次数"""numChanged = 0iterCount = 0alphas = np.mat(np.zeros((m, 1)))b = 0L = 0H = 0while iterCount < maxIter:numChanged = 0for i in range(m):Ei = E(i, alphas, b)Ri = y[i] * Ei# 选择违背KKT条件的，作为alpha2if (Ri < -tol and alphas[i] < C) or \(Ri > tol and alphas[i] > 0):# 选择第二个参数j = selectJRand(i)Ej = E(j, alphas, b)# get boundsif y[i] != y[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:continueKii = K[i, i]Kjj = K[j, j]Kij = K[i, j]eta = 2.0 * Kij - Kii - Kjjif eta >= 0:continueiOld = alphas[i].copy();jOld = alphas[j].copy()alphas[j] = jOld - y[j] * (Ei - Ej) / etaif alphas[j] > H:alphas[j] = Helif alphas[j] < L:alphas[j] = Lif abs(alphas[j] - jOld) < tol:alphas[j] = jOldcontinuealphas[i] = iOld + y[i] * y[j] * (jOld - alphas[j])# update bbINew = b - Ei - y[i] * (alphas[i] - iOld) * Kii - y[j] * \(alphas[j] - jOld) * KijbJNew = b - Ej - y[i] * (alphas[i] - iOld) * Kij - y[j] * \(alphas[j] - jOld) * Kjjif alphas[i] > 0 and alphas[i] < C:b = bINewelif alphas[j] > 0 and alphas[j] < C:b = bJNewelse:b = (bINew + bJNew) / 2.0numChanged += 1if numChanged == 0:iterCount += 1else:iterCount = 0supportVectorsIndex = np.nonzero(alphas.A > 0)[0]supportVectors = np.mat(X[supportVectorsIndex])return alphas, w(alphas, b, supportVectorsIndex, supportVectors), b, \supportVectorsIndex, supportVectors, iterCountreturn trainSimple, train, predict

参考资料

《机器学习实战》
《机器学习》

5.7 程序示例--基于 SMO 的 SVM 模型-机器学习笔记-斯坦福吴恩达教授相关推荐

7.3 程序示例--PCA 模型-机器学习笔记-斯坦福吴恩达教授
程序示例–PCA 模型 # coding: utf8 # pca/pca.pyimport numpy as npdef normalize(X):"""数据标准化处理A ...
5.5 SVM补充-机器学习笔记-斯坦福吴恩达教授
SVM补充决策边界 Coursera 上 ML 的课程对 SVM 介绍有限,参看了周志华教授的<机器学习>一书后,补充了当中对于 SVM 的介绍. 首先,我们考虑用更传统的权值定义式来描 ...
8.7 程序示例--异常检测-机器学习笔记-斯坦福吴恩达教授
程序示例–异常检测异常检测模型提供了一般高斯分布模型和多元高斯分布模型.其中,多元高斯分布模型被限制到了同轴分布: # coding: utf8 # anomaly_detection/anoma ...
5.11 程序示例--垃圾邮件检测-机器学习笔记-斯坦福吴恩达教授
程序示例–垃圾邮件检测邮件内容的预处理下面展示了一封常见的 email,邮件内容包含了一个 URL (http://www.rackspace.com/),一个邮箱地址(groupname-uns ...
5.10 程序示例--模型选择-机器学习笔记-斯坦福吴恩达教授
程序示例–模型选择在新的一组样本中,我们将通过交叉验证集选择模型,参数 CCC 和高斯核的参数 δδδ 我们都将在以下 8 个值中选取测试,则总共构成了 8×8=648×8=648×8=64 个模 ...
5.9 程序示例--非线性分类-机器学习笔记-斯坦福吴恩达教授
程序示例–非线性分类接下来,我们采用高斯核函数来解决非线性可分问题,由于数据集较大,我们使用性能更好的完整版 SMO 算法进行训练: # coding: utf8 # svm/test_non_li ...
5.8 程序示例--线性分类-机器学习笔记-斯坦福吴恩达教授
程序示例–线性分类首先,我们使用线性核函数来训练线性可分问题,这里,我们使用的是简化版 SMO 算法: # coding: utf8 # svm/test_linear import smo imp ...
9.5 程序示例--推荐系统-机器学习笔记-斯坦福吴恩达教授
程序示例–推荐系统推荐模型在推荐模型中,我们将暴露: 训练接口 train() 预测接口 predict(Theta, X) 获得推荐接口 getTopRecommends(Theta, X, i ...
7.5 程序示例--PCA for 数据可视化-机器学习笔记-斯坦福吴恩达教授
程序示例–PCA for 数据可视化我们有一张小鸟的图片,这是一个三通道彩色图像: 我们将图片的像素按颜色进行聚类,并在三维空间观察聚类成果: 似乎在三维空间可视化不是那么直观,借助于PCA,我们将 ...

5.7 程序示例--基于 SMO 的 SVM 模型-机器学习笔记-斯坦福吴恩达教授