SVM算法代码优化建议：提高效率和性能

"""#################################################\n# SVM: support vector machine\n# Author : yuehonggong\n# Date : 2023-05-06\n# Email : 105533927@qq.com\n#################################################"""\nfrom numpy import \nimport time\nimport matplotlib.pyplot as plt\nfrom six.moves import xrange\n\n\n# 计算内核值\ndef calcKernelValue(matrix_x, sample_x, kernelOption):\n kernelType = kernelOption[0]\n numSamples = matrix_x.shape[0]\n kernelValue = mat(zeros((numSamples, 1)))\n\n if kernelType == 'linear':\n kernelValue = matrix_x * sample_x.T\n elif kernelType == 'rbf':\n sigma = kernelOption[1]\n if sigma == 0:\n sigma = 1.0\n for i in xrange(numSamples):\n diff = matrix_x[i, :] - sample_x\n kernelValue[i] = exp(diff * diff.T / (-2.0 * sigma ** 2))\n else:\n raise NameError('Not support kernel type! You can use linear or rbf!')\n return kernelValue\n\n\n# 计算给定的训练集和核类型的核矩阵\ndef calcKernelMatrix(train_x, kernelOption):\n numSamples = train_x.shape[0]\n kernelMatrix = mat(zeros((numSamples, numSamples)))\n for i in xrange(numSamples):\n kernelMatrix[:, i] = calcKernelValue(train_x, train_x[i, :], kernelOption)\n return kernelMatrix\n\n\n# 定义仅用于存储变量和数据的结构\nclass SVMStruct:\n def init(self, dataSet, labels, C, toler, kernelOption):\n self.train_x = dataSet # each row stands for a sample\n self.train_y = labels # corresponding label\n self.C = C # slack variable\n self.toler = toler # termination condition for iteration\n self.numSamples = dataSet.shape[0] # number of samples\n self.alphas = mat(zeros((self.numSamples, 1))) # Lagrange factors for all samples\n self.b = 0\n self.errorCache = mat(zeros((self.numSamples, 2)))\n self.kernelOpt = kernelOption\n self.kernelMat = calcKernelMatrix(self.train_x, self.kernelOpt)\n\n\n# 计算阿尔法 k 的误差\ndef calcError(svm, alpha_k):\n output_k = float(multiply(svm.alphas, svm.train_y).T * svm.kernelMat[:, alpha_k] + svm.b)\n error_k = output_k - float(svm.train_y[alpha_k])\n return error_k\n\n\n# 在优化 Alpha K 后更新 Alpha K 的错误缓存\ndef updateError(svm, alpha_k):\n error = calcError(svm, alpha_k)\n svm.errorCache[alpha_k] = [1, error]\n\n\n# 选择步长最大的阿尔法J\ndef selectAlpha_j(svm, alpha_i, error_i):\n svm.errorCache[alpha_i] = [1, error_i] # mark as valid(has been optimized)\n candidateAlphaList = nonzero(svm.errorCache[:, 0].A)[0] # mat.A return array\n maxStep = 0; \n alpha_j = 0; \n error_j = 0\n\n # 找到具有最大迭代步长的 alpha\n if len(candidateAlphaList) > 1:\n for alpha_k in candidateAlphaList:\n if alpha_k == alpha_i:\n continue\n error_k = calcError(svm, alpha_k)\n if abs(error_k - error_i) > maxStep:\n maxStep = abs(error_k - error_i)\n alpha_j = alpha_k\n error_j = error_k\n # 如果第一次进入这个循环，我们随机选择 Alpha J\n else:\n alpha_j = alpha_i\n while alpha_j == alpha_i:\n alpha_j = int(random.uniform(0, svm.numSamples))\n error_j = calcError(svm, alpha_j)\n\n return alpha_j, error_j\n\n\n# 用于优化 alpha i 和 alpha j 的内部循环\ndef innerLoop(svm, alpha_i):\n error_i = calcError(svm, alpha_i)\n\n ### check and pick up the alpha who violates the KKT condition\n ## satisfy KKT condition\n # 1) yif(i) >= 1 and alpha == 0 (outside the boundary)\n # 2) yif(i) == 1 and 0<alpha< C (on the boundary)\n # 3) yif(i) <= 1 and alpha == C (between the boundary)\n ## violate KKT condition\n # because y[i]*E_i = y[i]*f(i) - y[i]^2 = y[i]*f(i) - 1, so\n # 1) if y[i]E_i < 0, so yif(i) < 1, if alpha < C, violate!(alpha = C will be correct)\n # 2) if y[i]E_i > 0, so yif(i) > 1, if alpha > 0, violate!(alpha = 0 will be correct)\n # 3) if y[i]E_i = 0, so yif(i) = 1, it is on the boundary, needless optimized\n if (svm.train_y[alpha_i] * error_i < -svm.toler) and (svm.alphas[alpha_i] < svm.C) or \n (svm.train_y[alpha_i] * error_i > svm.toler) and (svm.alphas[alpha_i] > 0):\n\n # step 1: 查找 alpha j\n alpha_j, error_j = selectAlpha_j(svm, alpha_i, error_i)\n alpha_i_old = svm.alphas[alpha_i].copy()\n alpha_j_old = svm.alphas[alpha_j].copy()\n\n # step 2: 计算 alpha j 的边界 L 和 H\n if svm.train_y[alpha_i] != svm.train_y[alpha_j]:\n L = max(0, svm.alphas[alpha_j] - svm.alphas[alpha_i])\n H = min(svm.C, svm.C + svm.alphas[alpha_j] - svm.alphas[alpha_i])\n else:\n L = max(0, svm.alphas[alpha_j] + svm.alphas[alpha_i] - svm.C)\n H = min(svm.C, svm.alphas[alpha_j] + svm.alphas[alpha_i])\n if L == H:\n return 0\n\n # step 3: 计算 ETA（样本 i 和 j 的相似性）\n eta = 2.0 * svm.kernelMat[alpha_i, alpha_j] - svm.kernelMat[alpha_i, alpha_i] \n - svm.kernelMat[alpha_j, alpha_j]\n if eta >= 0:\n return 0\n\n # step 4: 更新 alpha j\n svm.alphas[alpha_j] -= svm.train_y[alpha_j] * (error_i - error_j) / eta\n\n # step 5: 剪辑alpha j\n if svm.alphas[alpha_j] > H:\n svm.alphas[alpha_j] = H\n if svm.alphas[alpha_j] < L:\n svm.alphas[alpha_j] = L\n\n # step 6: 如果 Alpha J 移动不够，只需返回\n if abs(alpha_j_old - svm.alphas[alpha_j]) < 0.00001:\n updateError(svm, alpha_j)\n return 0\n\n # step 7: 优化后更新 alpha I alpha J\n svm.alphas[alpha_i] += svm.train_y[alpha_i] * svm.train_y[alpha_j] \n * (alpha_j_old - svm.alphas[alpha_j])\n\n # step 8: 更新阈值 B\n b1 = svm.b - error_i - svm.train_y[alpha_i] * (svm.alphas[alpha_i] - alpha_i_old) \n * svm.kernelMat[alpha_i, alpha_i] \n - svm.train_y[alpha_j] * (svm.alphas[alpha_j] - alpha_j_old) \n * svm.kernelMat[alpha_i, alpha_j]\n b2 = svm.b - error_j - svm.train_y[alpha_i] * (svm.alphas[alpha_i] - alpha_i_old) \n * svm.kernelMat[alpha_i, alpha_j] \n - svm.train_y[alpha_j] * (svm.alphas[alpha_j] - alpha_j_old) \n * svm.kernelMat[alpha_j, alpha_j]\n if (0 < svm.alphas[alpha_i]) and (svm.alphas[alpha_i] < svm.C):\n svm.b = b1\n elif (0 < svm.alphas[alpha_j]) and (svm.alphas[alpha_j] < svm.C):\n svm.b = b2\n else:\n svm.b = (b1 + b2) / 2.0\n\n # step 9: 在优化 Alpha I、J 和 B 后更新 alpha i、j 的错误缓存\n updateError(svm, alpha_j)\n updateError(svm, alpha_i)\n\n return 1\n else:\n return 0\n\n\n# 主要培训程序\ndef trainSVM(train_x, train_y, C, toler, maxIter, kernelOption=('rbf', 1.0)):\n # 计算训练时间\n startTime = time.time()\n\n # init data struct for svm\n svm = SVMStruct(mat(train_x), mat(train_y), C, toler, kernelOption)\n\n # start training\n entireSet = True\n alphaPairsChanged = 0\n iterCount = 0\n # Iteration termination condition:\n # Condition 1: reach max iteration\n # Condition 2: no alpha changed after going through all samples,\n # in other words, all alpha (samples) fit KKT condition\n while (iterCount < maxIter) and ((alphaPairsChanged > 0) or entireSet):\n alphaPairsChanged = 0\n\n # 更新所有训练示例的 Alpha\n if entireSet:\n for i in xrange(svm.numSamples):\n alphaPairsChanged += innerLoop(svm, i)\n print('---iter:%d entire set, alpha pairs changed:%d' % (iterCount, alphaPairsChanged))\n iterCount += 1\n # 在 alpha 不是 0 和不是 C（不在边界上）的示例中更新 alpha\n else:\n nonBoundAlphasList = nonzero((svm.alphas.A > 0) * (svm.alphas.A < svm.C))[0]\n for i in nonBoundAlphasList:\n alphaPairsChanged += innerLoop(svm, i)\n print('---iter:%d non boundary, alpha pairs changed:%d' % (iterCount, alphaPairsChanged))\n iterCount += 1\n\n # 所有示例和非边界示例的交替循环\n if entireSet:\n entireSet = False\n elif alphaPairsChanged == 0:\n entireSet = True\n\n print('Congratulations, training complete! Took %fs!' % (time.time() - startTime))\n return svm\n\n\n# testing your trained svm model given test set\ndef testSVM(svm, test_x, test_y):\n test_x = mat(test_x)\n test_y = mat(test_y)\n numTestSamples = test_x.shape[0]\n supportVectorsIndex = nonzero(svm.alphas.A > 0)[0]\n supportVectors = svm.train_x[supportVectorsIndex]\n supportVectorLabels = svm.train_y[supportVectorsIndex]\n supportVectorAlphas = svm.alphas[supportVectorsIndex]\n matchCount = 0\n for i in xrange(numTestSamples):\n kernelValue = calcKernelValue(supportVectors, test_x[i, :], svm.kernelOpt)\n predict = kernelValue.T * multiply(supportVectorLabels, supportVectorAlphas) + svm.b\n if sign(predict) == sign(test_y[i]):\n matchCount += 1\n accuracy = float(matchCount) / numTestSamples\n return accuracy\n\n\n# 显示仅可用于二维数据的已训练 SVM 模型\ndef showSVM(svm):\n if svm.train_x.shape[1] != 2:\n print("Sorry! I can not draw because the dimension of your data is not 2!")\n return 1\n\n # draw all samples\n for i in xrange(svm.numSamples):\n if svm.train_y[i] == -1:\n plt.plot(svm.train_x[i, 0], svm.train_x[i, 1], 'or')\n elif svm.train_y[i] == 1:\n plt.plot(svm.train_x[i, 0], svm.train_x[i, 1], 'ob')\n\n # 标记支持向量\n supportVectorsIndex = nonzero(svm.alphas.A > 0)[0]\n for i in supportVectorsIndex:\n plt.plot(svm.train_x[i, 0], svm.train_x[i, 1], 'oy')\n\n # 绘制分类线\n w = zeros((2, 1))\n for i in supportVectorsIndex:\n w += multiply(svm.alphas[i] * svm.train_y[i], svm.train_x[i, :].T)\n min_x = min(svm.train_x[:, 0])[0, 0]\n max_x = max(svm.train_x[:, 0])[0, 0]\n y_min_x = float(-svm.b - w[0] * min_x) / w[1]\n y_max_x = float(-svm.b - w[0] * max_x) / w[1]\n plt.plot([min_x, max_x], [y_min_x, y_max_x], '-g')\n plt.rcParams['font.sans-serif'] = ['KaiTi']\n plt.title("201530508 李涵博")\n plt.rcParams['axes.unicode_minus'] = False\n plt.show()\n\n\n以下是对代码进行优化的几点建议：\n\n1. 减少函数调用次数：代码中有许多函数调用，可以尝试减少函数的调用次数，以提高程序的执行效率。例如，在计算内核矩阵时，可以将计算内核值的函数直接嵌入到计算内核矩阵的函数中，避免多次调用。\n\n2. 使用向量化操作：在计算内核值、更新变量等涉及到矩阵运算的地方，可以利用numpy的向量化操作来提高运算速度。例如，可以使用numpy的multiply函数来替代循环中的乘法运算。\n\n3. 优化循环结构：代码中有多个循环结构，可以尝试优化循环的逻辑，减少循环次数。例如，在训练模型的过程中，可以根据alpha的变化情况来判断是否需要继续迭代，避免无效的循环。\n\n4. 使用更高效的数据结构：代码中使用了mat对象来表示矩阵，可以尝试使用numpy的ndarray来替代mat对象，以提高运算速度。\n\n5. 并行计算：如果计算量较大，可以考虑使用并行计算来加速程序的执行。可以使用多线程或多进程来并行计算，以提高程序的效率。\n\n6. 减少不必要的操作：可以检查代码中是否存在不必要的操作，如多余的条件判断、重复的计算等，可以尝试去除这些不必要的操作，以提高程序的执行效率。\n\n7. 使用更高效的算法：如果可能的话，可以考虑使用更高效的算法来解决问题。可以参考相关文献或其他资源，寻找更优的算法来替代当前的实现。\n\n以上是对代码进行优化的一些建议，具体可以根据实际情况进行调整和优化。\n