目錄
前言
一、引論
二、理論鋪墊
線性可分性（linear separability）
超平面
決策邊界
支持向量（support vector）
損失函數（loss function）
經驗風險（empirical risk）與結構風險（structural risk）
核方法
常見的核函數
三、算法流程
SMO序列最小優化算法
Python sklearn代碼實現：
Python源代碼實現+手寫字識別分類：
點關注，防走丟，如有纰漏之處，請留言指教，非常感謝
參閱：
博主大大小小參與過數十場數學建模比賽，SVM經常在各種建模比賽的優秀論文上見到該模型，一般直接使用SVM算法是比較少的，現在都是在此基礎理論之上提出優化算法。但是SVM的基礎理論是十分重要的思想，放眼整個分類算法中，SVM是最好的現成的分類器。這裡說的‘現成’指的是分類器不加修改即可直接使用。在神經網絡沒有出現之前，SVM的優化模型可以算得上是預測分類神器了，在機器學習中SVM仍舊是最為出名的算法之一了，本篇博客將致力於將SVM算法以及原理每一個知識點都講明白，希望沒有講明白的點大家可以在評論區指出。
我們使用SVM支持向量機一般用於分類，得到低錯誤率的結果。SVM能夠對訓練集意外的數據點做出很好的分類決策。那麼首先我們應該從數據層面上去看SVM到底是如何做決策的，這裡來看這樣一串數據集集合在二維平面坐標系上描繪的圖：
現在我們需要考慮，是否能夠畫出一條直線將圓形點和星星點分開。像first第一張圖片來看，圓點和星點就分的很開，很容易就可以在圖中畫出一條直線將兩組數據分開。而看第二張圖片，圓點和星點幾乎都聚合在一起，要區分的話十分困難。
我們要劃線將他們區分開來的話，有有無數條可以畫，但是我們難以找到一條最好區分度最高的線條將它們幾乎完全區分。那麼在此我們需要了解兩個關於數據集的基本概念：
第一張圖片的數據就是線性可分的，肉眼可見的線性可分。而第二張圖片的數據一看就是線性不可分的。當然我們這樣直率的觀測也沒有錯，但是這僅是建立在二維平面數據可視化的基礎之上，若是特征維度更高，例如：
超出我們可視化的技術之外我們如何判斷他們是否可分呢？
所有要判斷線性可分性需要對線性可分性下個定義：
在分類問題中給定輸入數據和學習目標：
,
,其中輸入數據的每個樣本都包含多個特征並由此構成特征空間（feature space）：
,而學習目標為二元變量
表示負類（negative class）和正類（positive class）。
若輸入數據所在的特征空間存在作為決策邊界（decision boundary）的超平面將學習目標按正類和負類分開，並使任意樣本的點到平面距離(point to plane distance)大於等於1：
decision boundary：
point to plane distance:
則稱該分類問題具有線性可分性，參數
,分別為超平面的法向量和截距。
到這裡肯定有很多小伙伴會問啥是決策邊界什麼是超平面啊，莫慌馬上給出定義：
而對機器學習來說，涉及的多是高維空間（多維度）的數據分類，高維空間的SVM，即為超平面。機器學習的最終目的就是要找到最合適的（也即最優的）一個分類超平面（Hyper plane），從而應用這個最優分類超平面將特征數據很好地區分為兩類。
大家看圖就懂了：
在數學中，超平面（Hyperplane）是n維歐式空間中余維度等於1的線性子空間。這是平面中的直線、空間中的平面之推廣。設F為域（可考慮
）。
n維空間
中的超平面是由方程：
定義的子集，其中
是不全為零的常數。
在線性代數的脈絡下，F-矢量空間V中的超平面是指形如：
的子空間，其中
是任一非零的線性映射。
在射影幾何中，同樣可定義射影空間
中的超平面。在齊次坐標（
)下超平面可由以下方程定義：
，其中
是不全為零的常數。
總而言之：超平面H是從n維空間到n-1維空間的一個映射子空間，它有一個n維向量和一個實數定義。設d是n維歐式空間R中的一個非零向量，a是實數，則R中滿足條件dX=a的點X所組成的集合稱為R中的一張超平面。
SVM是一種優化的分類算法，其動機是尋找一個最佳的決策邊界，使得從決策邊界與各組數據之間存在margin，並且需要使各側的margin最大化。那麼這個決策邊界就是不同類之間的界限。
總而言之：在具有兩個類的統計分類問題中，決策邊界或決策表面是超平面，其將基礎向量空間劃分為兩個集合，一個集合。 分類器將決策邊界一側的所有點分類為屬於一個類，而將另一側的所有點分類為屬於另一個類。
在了解了超平面和決策邊界我們發現SVM的核心任務是找到一個超平面作為決策邊界。那麼滿足該條件的決策邊界實際上構造了2個平行的超平面作為間隔邊界以判別樣本的分類：
所有在上間隔邊界上方的樣本屬於正類，在下間隔邊界下方的樣本屬於負類。兩個間隔邊界的距離
被定義為邊距（margin），位於間隔邊界上的正類和負類樣本為支持向量（support vector）。
拿我們的first的數據集來看：
無論我們怎麼畫直線總有損失的點沒有正確分類到 。在一個分類問題不具有線性可分性時，使用超平面作為決策邊界會帶來分類損失，即部分支持向量不再位於間隔邊界上，而是進入了間隔邊界內部，或落入決策邊界的錯誤一側。
損失函數可以對分類損失進行量化，其按數學意義可以得到的形式是0-1損失函數：
0-1損失函數是一個不連續的分段函數，不利於求解其最小化問題，因此在應用可構造其代理損失（surrogate loss）。代理損失是與原損失函數具有相合性（consistency）的損失函數，最小化代理損失所得的模型參數也是最小化原損失函數的解。也就是說我們需要損失函數計算最小的決策邊界。這裡給出二元分類（binary classification）中0-1損失函數的代理損失：
名稱表達式鉸鏈損失函數（hinge loss function）
交叉熵損失函數（cross-entropy loss function）
指數損失函數（exponential loss function）
SVM使用的是鉸鏈損失函數:
.
按統計學習理論，分類器在經過學習並應用於新數據時會產生風險，風險的類型可分為經驗風險和結構風險。
empirical risk：
structural risk：
式中
表示分類器，經驗風險由損失函數定義，描述了分類器所給出的分類結果的准確程度；結構風險由分類器參數矩陣的范數定義，描述了分類器自身的復雜程度以及穩定程度，復雜的分類器容易產生過擬合，因此是不穩定的。若一個分類器通過最小化經驗風險和結構風險的線性組合以確定其模型參數：
則對該分類器的求解是一個正則化問題，常數
是正則化系數。當
=2時，該式被稱為
正則化或Tikhonov正則化（Tikhonov regularization）。SVM的結構風險
表示，在線性可分問題下硬邊界SVM的經驗風險可以歸0，因此其是一個完全最小化結構風險的分類器；在線性不可分問題中，軟邊界SVM的經驗風險不可歸0，因此其是一個正則化分類器，最小化結構風險和經驗風險的線性組合。
一些線性不可分的問題可能是非線性可分的，即特征空間存在超平面將正類和負類分開。使用非線性函數可以將非線性可分問題從原始的特征空間映射至更高維的空間（希爾伯特空間（Hilbert space）
）從而轉化為線性可分問題，此時作為決策邊界的超平面表示如下：
式中
為映射函數。由於映射函數具有復雜的形式，難以計算其內積，因此可使用核方法（kernel method），即定義映射函數的內積為核函數：
以回避內積的顯式計算。
徑向基函數核（RBF kernel）
拉普拉斯核（Laplacian kernel）
Sigmoid核（Sigmoid kernel）當多項式核的階為1時，其被稱為線性核，對應的非線性分類器退化為線性分類器。RBF核也被稱為高斯核（Gaussian kernel），其對應的映射函數將樣本空間映射至無限維空間。
在了解了上述算法原理以後，我們便對SVM算法有了個大致清晰的認知，那麼SVM是如何挖掘到損失最小的超平面的呢？
SVM的求解可以使用二次凸優化問題的數值方法，例如內點法（Interior Point Method, IPM）和序列最小優化算法（Sequential Minimal Optimization, SMO），在擁有充足學習樣本時也可使用隨機梯度下降。
這裡我們采用SMO序列最小優化算法計算：
SMO是一種坐標下降法（coordinate descent）以迭代方式求解SVM的對偶問題，其設計是在每個迭代步選擇拉格朗日乘子中的兩個變量
並固定其他參數，將原優化問題簡化至一維子可行域，此時約束條件有如下等價形式 ：
將上式右側帶入SVM的對偶問題並消去求和項中的
可以得到僅關於
的二次規劃問題，該優化問題有閉式解可以快速計算。在此基礎上，SMO有如下計算框架:
可以證明，在二次凸優化問題中，SMO的每步迭代都嚴格地優化了SVM的對偶問題，且迭代會在有限步後收斂於全局極大值 。SMO算法的迭代速度與所選取乘子對KKT條件的偏離程度有關，因此SMO通常采用啟發式方法選取拉格朗日乘子。
為了方便展示效果，我們仍舊用first的數據集進行：
sklearn.svm.SVC語法格式為：
class sklearn.svm.SVC( *,
C=1.0,
kernel='rbf',
degree=3,
gamma='scale',
coef0=0.0,
shrinking=True,
probability=False,
tol=0.001,
cache_size=200,
class_weight=None,
verbose=False,
max_iter=- 1,
decision_function_shape='ovr',
break_ties=False,
random_state=None)該實現基於 libsvm，更多詳細信息可以去官網上查閱：sklearn.svm.SVC.
# 導入模塊
import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm, datasets
# 鸢尾花數據
iris = datasets.load_iris()
X = iris.data[:, :2] # 為便於繪圖僅選擇2個特征
y = iris.target
# 測試樣本（繪制分類區域）
xlist1 = np.linspace(X[:, 0].min(), X[:, 0].max(), 200)
xlist2 = np.linspace(X[:, 1].min(), X[:, 1].max(), 200)
XGrid1, XGrid2 = np.meshgrid(xlist1, xlist2)
# 非線性SVM：RBF核，超參數為0.5，正則化系數為1，SMO迭代精度1e-5, 內存占用1000MB
svc = svm.SVC(kernel='rbf', C=1, gamma=0.5, tol=1e-5, cache_size=1000).fit(X, y)
# 預測並繪制結果
Z = svc.predict(np.vstack([XGrid1.ravel(), XGrid2.ravel()]).T)
Z = Z.reshape(XGrid1.shape)
plt.contourf(XGrid1, XGrid2, Z, cmap=plt.cm.hsv)
plt.contour(XGrid1, XGrid2, Z, colors=('k',))
plt.scatter(X[:, 0], X[:, 1], c=y, edgecolors='k', linewidth=1.5, cmap=plt.cm.hsv)
plt.show() 
from numpy import *
def selectJrand(i, m): # 在某個區間范圍內隨機選擇一個整數
# i為第一個alhpa的下表，m是所有alpha的數目
j = i # we want to select any J not equal to i
while (j == i):
j = int(random.uniform(0, m))
return j
def clipAlpha(aj, H, L): # 在數值太大的時候對其進行調整
# aj是H是下限，是L的上限
if aj > H:
aj = H
if L > aj:
aj = L
return aj
def kernelTrans(X, A, kTup): # calc the kernel or transform data to a higher dimensional space
# X是數據，A是
m, n = shape(X)
K = mat(zeros((m, 1)))
# 這裡為了簡單，我們只內置了兩種核函數，但是原理是一樣的，需要的話再寫其他類型就是了
# 線性核函數：k(x,x_i)=x*x_i,它不需要再傳入參數。
if kTup[0] == 'lin':
K = X * A.T # linear kernel，線性核函數
elif kTup[0] == 'rbf': # 高斯核函數，傳入的參數作為detla
for j in range(m):
deltaRow = X[j, :] - A
K[j] = deltaRow * deltaRow.T # 2范數
K = exp(K / (-1 * kTup[1] ** 2)) # divide in NumPy is element-wise not matrix like Matlab
# numpy中，/表示對矩陣元素進行計算而不是計算逆（MATLAB）
else:
raise NameError('Houston We Have a Problem -- \
That Kernel is not recognized')
return K
# 定義了一個類來進行SMO算法
class optStruct:
def __init__(self, dataMatIn, classLabels, C, toler, kTup): # Initialize the structure with the parameters
# kTup儲存核函數信息，它是一個元組，元組第一個元素是一個描述核函數類型的字符串，其他兩個元素是核函數可能需要的參數
self.X = dataMatIn
self.labelMat = classLabels
self.C = C
self.tol = toler
self.m = shape(dataMatIn)[0] # m是樣本個數，也是a的個數
self.alphas = mat(zeros((self.m, 1))) # 初始化a序列，都設置為0
self.b = 0
# 第一列給出的是eCache是否有效的標志位，而第二位是實際的E值
self.eCache = mat(zeros((self.m, 2))) # first column is valid flag
# 如果第一位是1，說明現在的這個Ek是有效的
self.K = mat(zeros((self.m, self.m))) # 使用核函數計算後的數據，就是內積矩陣，方便直接調用內積結果
for i in range(self.m):
self.K[:, i] = kernelTrans(self.X, self.X[i, :], kTup)
def calcEk(oS, k): # 計算第k個樣本的Ek
fXk = float(multiply(oS.alphas, oS.labelMat).T * oS.K[:, k] + oS.b)
Ek = fXk - float(oS.labelMat[k])
return Ek
def selectJ(i, oS, Ei): # this is the second choice -heurstic, and calcs Ej
# 選定了a_1之後選擇a_2
# 選擇a_2
maxK = -1;
maxDeltaE = 0;
Ej = 0 # 選擇|E1-E2|最大的E2並返回E2和j
# 先將E1存進去，以便於之後的統一化進行
oS.eCache[i] = [1, Ei] # set valid #choose the alpha that gives the maximum delta E
'''
numpy函數返回非零元素的目錄。
返回值為元組， 兩個值分別為兩個維度， 包含了相應維度上非零元素的目錄值。
可以通過a[nonzero(a)]來獲得所有非零值。
.A的意思是：getArray(),也就是將矩陣轉換為數組
'''
# 獲取哪些樣本的Ek是有效的，ValidEcacheList裡面存的是所有有效的樣本行Index
validEcacheList = nonzero(oS.eCache[:, 0].A)[0]
# 對每一個有效的Ecache都比較一遍
if (len(validEcacheList)) > 1:
for k in validEcacheList: # loop through valid Ecache values and find the one that maximizes delta E
if k == i: continue # don't calc for i, waste of time
# 如果選到了和a1一樣的，就繼續，因為a1和a2必須選不一樣的樣本
Ek = calcEk(oS, k)
deltaE = abs(Ei - Ek)
if (deltaE > maxDeltaE):
maxK = k;
maxDeltaE = deltaE;
Ej = Ek
return maxK, Ej
else: # in this case (first time around) we don't have any valid eCache values
j = selectJrand(i, oS.m)
Ej = calcEk(oS, j)
return j, Ej
def updateEk(oS, k): # after any alpha has changed update the new value in the cache
Ek = calcEk(oS, k)
oS.eCache[k] = [1, Ek]
# 計算Ek並保存在類oS中
# 內循環尋找合適的a_2
def innerL(i, oS):
Ei = calcEk(oS, i) # 為什麼這裡要重新算呢？因為a_1剛剛更新了，和存儲的不一樣
if ((oS.labelMat[i] * Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or (
(oS.labelMat[i] * Ei > oS.tol) and (oS.alphas[i] > 0)):
# 如果a1選的合適的話，不合適就直接結束了
# 剩下的邏輯都一樣，只不過不是使用x_ix_j，而是使用核函數
j, Ej = selectJ(i, oS, Ei) # this has been changed from selectJrand
alphaIold = oS.alphas[i].copy();
alphaJold = oS.alphas[j].copy();
if (oS.labelMat[i] != oS.labelMat[j]):
L = max(0, oS.alphas[j] - oS.alphas[i])
H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
else:
L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
H = min(oS.C, oS.alphas[j] + oS.alphas[i])
if L == H: print("L==H"); return 0
eta = 2.0 * oS.K[i, j] - oS.K[i, i] - oS.K[j, j] # changed for kernel
if eta >= 0: print("eta>=0"); return 0
oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej) / eta
oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
updateEk(oS, j) # added this for the Ecache
if (abs(oS.alphas[j] - alphaJold) < 0.00001): print("j not moving enough"); return 0
oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alphaJold - oS.alphas[j]) # update i by the same amount as j
updateEk(oS, i) # added this for the Ecache #the update is in the oppostie direction
b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.K[i, i] - oS.labelMat[j] * (
oS.alphas[j] - alphaJold) * oS.K[i, j]
b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.K[i, j] - oS.labelMat[j] * (
oS.alphas[j] - alphaJold) * oS.K[j, j]
if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]):
oS.b = b1
elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]):
oS.b = b2
else:
oS.b = (b1 + b2) / 2.0
return 1
else:
return 0
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)): # 默認核函數是線性，參數為0（那就是它本身了）
# 這個kTup先不管，之後用
oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler, kTup) # 初始化這一結構
iter = 0
# entireSet是控制開關，一次循環對所有樣本點都遍歷，第二次就只遍歷非邊界點
entireSet = True;
alphaPairsChanged = 0
while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
alphaPairsChanged: int = 0
if entireSet: # go over all
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: # go over non-bound (railed) alphas
# 把大於0且小於C的a_i挑出來，這些是非邊界點，只從這些點上遍歷
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): # 如果alpha沒有更新,計算全樣本遍歷
entireSet = True
print("iteration number: %d" % iter)
return oS.b, oS.alphas
# 利用SVM進行分類，返回的是函數間隔，大於0屬於1類，小於0屬於2類。
def calcWs(alphas, dataArr, classLabels):
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 img2vector(filename):
# 一共有1024個特征
returnVect = zeros((1, 1024))
fr = open(filename)
for i in range(32):
lineStr = fr.readline()
for j in range(32):
returnVect[0, 32 * i + j] = int(lineStr[j])
return returnVect
def loadImages(dirName):
from os import listdir
hwLabels = []
# 利用listdir讀文件名，這裡的label寫在了文件名裡面
trainingFileList = listdir(dirName) # load the training set
m = len(trainingFileList)
trainingMat = zeros((m, 1024))
for i in range(m):
fileNameStr = trainingFileList[i]
fileStr = fileNameStr.split('.')[0] # take off .txt
classNumStr = int(fileStr.split('_')[0])
if classNumStr == 9:
hwLabels.append(-1)
else:
hwLabels.append(1)
trainingMat[i, :] = img2vector('%s/%s' % (dirName, fileNameStr))
return trainingMat, hwLabels
# 手寫識別問題
def testDigits(kTup=('rbf', 10)):
dataArr, labelArr = loadImages('trainingDigits')
b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, kTup)
datMat = mat(dataArr);
labelMat = mat(labelArr).transpose()
svInd = nonzero(alphas.A > 0)[0]
sVs = datMat[svInd]
labelSV = labelMat[svInd];
print("there are %d Support Vectors" % shape(sVs)[0])
m, n = shape(datMat)
errorCount = 0
for i in range(m):
kernelEval = kernelTrans(sVs, datMat[i, :], kTup)
predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
if sign(predict) != sign(labelArr[i]): errorCount += 1
print("the training error rate is: %f" % (float(errorCount) / m))
dataArr, labelArr = loadImages('testDigits')
errorCount = 0
datMat = mat(dataArr);
labelMat = mat(labelArr).transpose()
m, n = shape(datMat)
for i in range(m):
kernelEval = kernelTrans(sVs, datMat[i, :], kTup)
predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
if sign(predict) != sign(labelArr[i]): errorCount += 1
print("the test error rate is: %f" % (float(errorCount) / m))
testDigits(kTup=('rbf', 20))
以上就是本期全部內容。我是fanstuck ，有問題大家隨時留言討論 ，我們下期見。
SVM通俗詳解
SVM百度百科
損失函數百度百科
超平面百度百科
Computer graduation design Python+djang library book borrowing and returning management system (source code + system + mysql database + Lw document)
項目介紹論文闡述了圖書管理系統,並對該系統的需求分析及系統需
