目前玩机器学习的小伙伴，上来就是使用现有的sklearn机器学习包，写两行代码，调调参数就能跑起来，看似方便，实则有时不利于个人能力发展，要知道现在公司需要的算法工程师，不仅仅只是会调参（这种工作，入门几个月的人就可以干了），而是要深入底层，能优化代码，能自己搭。
本文章适合以下几类人：
1）初学者，了解机器学习的实现过程
2）想提升自己的代码能力
第一步：原理
什么是SVM？
在机器学习中，支持向量机 (英语: support vector machine，常简称为SVM，又名支持向量网络) 是在分类与回归分析中分析数据的监督式学习模型与相关的学习算法。给定一组训练实例，每个训练实例被标记为属于两个类别中的一个或另一个，SVM训练算法创 建一个将新的实例分配给两个类别之一的模型，使其成为非概率二元线性分类器。SVM模型是将实例表示为空间中的点，这样映射就使得单独类别的实例被尽可能宽的明显的间隔分开。然后，将新的实例映射到同一空间，并基于它们落在间隔的哪一侧来预测所属类别。
除了进行线性分类之外，SVM还可以使用所谓的核技巧有效地进行非线性分类，将其输入隐式映射到高维特征空间中。
当数据末被标记时，不能进行监督式学习，需要用非监督式学习，它会尝试找出数据到簇的自然聚类，并将新数据映射到这些已形成的簇。将支持向量机改进的聚类算法被称为支持向量聚类，当数据末被标记或者仅一些数据被标记时，支持向量聚类经常在工业应用中用作分类步骤的预处理。
SVM的推导，百度一下很多，这里就做一下搬运工了，可参考：【机器学习实战】线性支持向量机Python实现-CSDN博客
以下为线性核SVM
以下为rbf核SVM
第二步：代码实现（有简化版SMO，原始SMO，基于核的SMO）
#include<math.h>
#include<iostream>
#include<string.h>
#include<vector>
#include<map>
#include<time.h>
#include <stdlib.h>
#include <string>
#include <fstream>
#include <sstream>
#include <random>
#include<time.h>
#include "matrix.h"
using namespace std;
double max(double a, double b) {
return (a > b ? a : b);
}
double min(double a, double b) {
return (a < b ? a : b);
}
/***没有优化的版本***/
//生成随机数函数
inline int selectJrand(const int i, const int m)
{
int randVal = i;
while(randVal == i)
{
randVal = rand() % (m);
}
return randVal;
}
//阈值函数
inline double clipAlpha(double &aj, const double &H, const double &L)
{
if (aj > H)
{
aj = H;
}
if (aj < L)
{
aj = L;
}
return aj;
}
//data数据里面第i行乘以第j行的转置，返回一个double型数值
double dataIJ(Matrix data, int i, int j)
{
double dataVal = 0.0;
Matrix dataI;
Matrix dataJ;
Matrix dataJT;
dataI.initMatrix(&dataI, 1, data.row);
dataJ.initMatrix(&dataJ, 1, data.row);
dataJT.initMatrix(&dataJT, data.row, 1);
dataI = dataI.getOneCol(data, i);
dataJ = dataJ.getOneCol(data, j);
dataJT.transposematrix(dataJ, &dataJT);
for (int k = 0; k < data.row; k++)
{
dataVal += dataI.mat[0][k] * dataJT.mat[k][0];
}
return dataVal;
}
//简化版SMO算法 返回b，alpha则用引用
double smoSimple(Matrix dataMatIn, Matrix classLabels, Matrix &alphas, const double C, const double toler, const int maxIter)
{
}
/***没有优化的版本***/
/***优化的版本***/
//参数保持在数据结构体中
struct OS{
Matrix X;
Matrix labelMat;
double C;
double tol;
int m;
Matrix alpahs;
double b;
Matrix eCache;
};
class SMOP{
public:
OS os;
public:
int initOs(Matrix dataMatIn, Matrix classLabels, double C, double toler)
{
os.X.initMatrix(&os.X, dataMatIn.col, dataMatIn.row);
os.X.copy(dataMatIn, &os.X);
os.labelMat.initMatrix(&os.labelMat, classLabels.col, classLabels.row);
os.labelMat.copy(classLabels, &os.labelMat);
os.C = C;
os.tol = toler;
os.m = dataMatIn.col;
os.alpahs.initMatrix(&os.alpahs, os.m, 1, 0.0);
os.b = 0.0;
os.eCache.initMatrix(&os.eCache, os.m, 2, 0.0);
return 0;
}
double calcEk(int k)
{
double Ek = 0.0;
double fXk = 0.0;
Matrix alphasMultsLabel;
alphasMultsLabel.initMatrix(&alphasMultsLabel, os.m, 1);
alphasMultsLabel.multsmatrixDot(&alphasMultsLabel, os.alpahs, os.labelMat);
Matrix alphasMultsLabelT;
alphasMultsLabelT.initMatrix(&alphasMultsLabelT, 1, os.m);
alphasMultsLabelT.transposematrix(alphasMultsLabel, &alphasMultsLabelT);
Matrix osXk;
osXk.initMatrix(&osXk, 1, os.X.row);
osXk = osXk.getOneCol(os.X, k);
Matrix osXkT;
osXkT.initMatrix(&osXkT, os.X.row, 1);
osXkT.transposematrix(osXk, &osXkT);
Matrix osXMultisOsXkT;
osXMultisOsXkT.initMatrix(&osXMultisOsXkT, os.X.col, 1);
osXMultisOsXkT.multsmatrix(&osXMultisOsXkT, os.X, osXkT);
for (int i = 0; i < os.m; i++)
{
fXk += alphasMultsLabelT.mat[0][i] * osXMultisOsXkT.mat[i][0];
}
fXk += os.b;
Ek = fXk - os.labelMat.mat[k][0];
return Ek;
}
//如果处理返回 maxK和Ej 否则返回j和Ej
double selectJ(int i, double Ei , double &Ej)
{
int maxK = os.m - 1;
double maxDeltaE = 0.0;
Ej = 0.0;
os.eCache.mat[i][0] = 1;
os.eCache.mat[i][1] = Ei;
vector<int>validEcacheList;
//统计第一列不为0的序号
for (int k = 0; k < os.m; k++)
{
if (os.eCache.mat[k][0] != 0)
{
validEcacheList.push_back(k);
}
}
if (validEcacheList.size() > 1)
{
for (int k = 0; k < validEcacheList.size(); k++)
{
if (validEcacheList[k] == i)
{
continue;
}
double Ek = calcEk(validEcacheList[k]);
double deltaE = abs(Ei - Ek);
if (deltaE > maxDeltaE)
{
maxK = validEcacheList[k];
maxDeltaE = deltaE;
Ej = Ek;
}
}
return maxK;
}
else
{
int j = selectJrand(i, os.m);
Ej = calcEk(j);
return j;
}
}
void updateEk(int k)
{
double Ek = calcEk(k);
os.eCache.mat[k][0] = 1;
os.eCache.mat[k][1] = Ek;
}
int innerL(int i)
{
double Ei = calcEk(i);
if (((os.labelMat.mat[i][0] * Ei < -os.tol) && (os.alpahs.mat[i][0] < os.C)) ||
((os.labelMat.mat[i][0] * Ei > os.tol) && (os.alpahs.mat[i][0] > 0)))
{
double Ej = 0.0;
int j = selectJ(i, Ei, Ej);
double alphaIold = os.alpahs.mat[i][0];
double alphaJold = os.alpahs.mat[j][0];
double L = 0.0;
double H = 0.0;
if (os.labelMat.mat[i][0] != os.labelMat.mat[j][0])
{
L = max(0.0, os.alpahs.mat[j][0] - os.alpahs.mat[i][0]);
H = min(os.C, os.C + os.alpahs.mat[j][0] - os.alpahs.mat[i][0]);
}
else
{
L = max(0.0, os.alpahs.mat[j][0] + os.alpahs.mat[i][0] - os.C);
H = min(os.C, os.alpahs.mat[j][0] + os.alpahs.mat[i][0]);
}
if (abs(L - H) < 0.00001)
{
cout<<"L == H"<<endl;
return 0;
}
double eta = 0.0;
eta += 2.0 * dataIJ(os.X, i, j);
eta -= dataIJ(os.X, i, i);
eta -= dataIJ(os.X, j, j);
if (eta >= 0)
{
cout<<"eta >= 0"<<endl;
return 0;
}
os.alpahs.mat[j][0] -= os.labelMat.mat[j][0] * (Ei - Ej) / eta;
os.alpahs.mat[j][0] = clipAlpha(os.alpahs.mat[j][0], H, L);
updateEk(j);
if (abs(os.alpahs.mat[j][0] - alphaJold) < 0.00001)
{
cout<<"j not moving enought"<<endl;
return 0;
}
os.alpahs.mat[i][0] += os.labelMat.mat[j][0] * os.labelMat.mat[i][0] *
(alphaJold - os.alpahs.mat[j][0]);
updateEk(i);
double b1 = os.b - Ei;
b1 -= os.labelMat.mat[i][0] * (os.alpahs.mat[i][0] - alphaIold) * dataIJ(os.X, i, i);
b1 -= os.labelMat.mat[j][0] * (os.alpahs.mat[j][0] - alphaJold) * dataIJ(os.X, i, j);
double b2 = os.b - Ej;
b2 -= os.labelMat.mat[i][0] * (os.alpahs.mat[i][0] - alphaIold) * dataIJ(os.X, i, j);
b2 -= os.labelMat.mat[j][0] * (os.alpahs.mat[j][0] - alphaJold) * dataIJ(os.X, j, j);
if ((0 < os.alpahs.mat[i][0]) && (os.C > os.alpahs.mat[i][0]))
{
os.b = b1;
}
else if ((0 < os.alpahs.mat[j][0]) && (os.C > os.alpahs.mat[j][0]))
{
os.b = b2;
}
else
os.b = (b1 + b2) / 2.0;
return 1;
}
else
return 0;
}
double smoP(Matrix dataMatIn, Matrix classLabels, Matrix &alphas, double C, double toler,int maxIter)
{
}
};
/***优化的版本+核函数***/
//参数保持在数据结构体中
struct OSK{
Matrix X;
Matrix labelMat;
double C;
double tol;
int m;
Matrix alpahs;
double b;
Matrix eCache;
Matrix k;
};
struct KTup{
string stringType;
double val;
};
class SMOPK{
public:
OSK os;
public:
int initOs(Matrix dataMatIn, Matrix classLabels, double C, double toler, KTup kTup)
{
os.X.initMatrix(&os.X, dataMatIn.col, dataMatIn.row);
os.X.copy(dataMatIn, &os.X);
os.labelMat.initMatrix(&os.labelMat, classLabels.col, classLabels.row);
os.labelMat.copy(classLabels, &os.labelMat);
os.C = C;
os.tol = toler;
os.m = dataMatIn.col;
os.alpahs.initMatrix(&os.alpahs, os.m, 1, 0.0);
os.b = 0.0;
os.eCache.initMatrix(&os.eCache, os.m, 2, 0.0);
os.k.initMatrix(&os.k, os.m, os.m);
Matrix Ktemp;
Ktemp.initMatrix(&Ktemp, os.m, 1);
Matrix A;
A.initMatrix(&A, 1, os.X.row);
for (int i = 0; i < os.m; i++)
{
A = A.getOneCol(os.X, i);
Ktemp = kernelTrans(os.X, A, kTup);
cout<<endl;
for (int j = 0; j < os.m; j++)
{
os.k.mat[j][i] = Ktemp.mat[j][0];
}
}
os.k.print(os.k);
return 0;
}
Matrix kernelTrans(Matrix X, Matrix &A, KTup kTup)
{
int m = X.col;
int n = X.row;
Matrix K;
K.initMatrix(&K,m,1);
Matrix AT;
AT.initMatrix(&AT, A.row, 1);
AT.transposematrix(A, &AT);
if (kTup.stringType == "lin")
{
K.multsmatrix(&K, X, AT);
}
else if (kTup.stringType == "rbf")
{
for (int j = 0; j < m; j++)
{
for (int i = 0; i < n; i++)
{
K.mat[j][0] += (X.mat[j][i] - A.mat[0][i]) * (X.mat[j][i] - A.mat[0][i]);
}
K.mat[j][0] = exp(K.mat[j][0] / (-1 * pow(kTup.val,2)));
}
}
else
{
cout<<"That Kernel is not recognized"<<endl;
exit(-1);
}
return K;
}
double calcEk(int k)
{
double Ek = 0.0;
double fXk = 0.0;
Matrix alphasMultsLabel;
alphasMultsLabel.initMatrix(&alphasMultsLabel, os.m, 1);
alphasMultsLabel.multsmatrixDot(&alphasMultsLabel, os.alpahs, os.labelMat);
Matrix alphasMultsLabelT;
alphasMultsLabelT.initMatrix(&alphasMultsLabelT, 1, os.m);
alphasMultsLabelT.transposematrix(alphasMultsLabel, &alphasMultsLabelT);
Matrix osXk;
osXk.initMatrix(&osXk, os.X.col, 1);
osXk = osXk.getOneRow( os.k, k); //获得第k列
for (int i = 0; i < os.m; i++)
{
fXk += alphasMultsLabelT.mat[0][i] * osXk.mat[i][0];
}
fXk += os.b;
Ek = fXk - os.labelMat.mat[k][0];
return Ek;
}
//如果处理返回 maxK和Ej 否则返回j和Ej
double selectJ(int i, double Ei , double &Ej)
{
int maxK = os.m - 1;
double maxDeltaE = 0.0;
Ej = 0.0;
os.eCache.mat[i][0] = 1;
os.eCache.mat[i][1] = Ei;
vector<int>validEcacheList;
//统计第一列不为0的序号
for (int k = 0; k < os.m; k++)
{
if (os.eCache.mat[k][0] != 0)
{
validEcacheList.push_back(k);
}
}
if (validEcacheList.size() > 1)
{
for (int k = 0; k < validEcacheList.size(); k++)
{
if (validEcacheList[k] == i)
{
continue;
}
double Ek = calcEk(validEcacheList[k]);
double deltaE = abs(Ei - Ek);
if (deltaE > maxDeltaE)
{
maxK = validEcacheList[k];
maxDeltaE = deltaE;
Ej = Ek;
}
}
return maxK;
}
else
{
int j = selectJrand(i, os.m);
Ej = calcEk(j);
return j;
}
}
void updateEk(int k)
{
double Ek = calcEk(k);
os.eCache.mat[k][0] = 1;
os.eCache.mat[k][1] = Ek;
}
int innerL(int i)
{
double Ei = calcEk(i);
if (((os.labelMat.mat[i][0] * Ei < -os.tol) && (os.alpahs.mat[i][0] < os.C)) ||
((os.labelMat.mat[i][0] * Ei > os.tol) && (os.alpahs.mat[i][0] > 0)))
{
double Ej = 0.0;
int j = selectJ(i, Ei, Ej);
double alphaIold = os.alpahs.mat[i][0];
double alphaJold = os.alpahs.mat[j][0];
double L = 0.0;
double H = 0.0;
if (os.labelMat.mat[i][0] != os.labelMat.mat[j][0])
{
L = max(0.0, os.alpahs.mat[j][0] - os.alpahs.mat[i][0]);
H = min(os.C, os.C + os.alpahs.mat[j][0] - os.alpahs.mat[i][0]);
}
else
{
L = max(0.0, os.alpahs.mat[j][0] + os.alpahs.mat[i][0] - os.C);
H = min(os.C, os.alpahs.mat[j][0] + os.alpahs.mat[i][0]);
}
if (abs(L - H) < 0.00001)
{
cout<<"L == H"<<endl;
return 0;
}
double eta = 0.0;
eta += 2.0 * os.k.mat[i][j];
eta -= os.k.mat[i][i];
eta -= os.k.mat[j][j];
if (eta >= 0)
{
cout<<"eta >= 0"<<endl;
return 0;
}
os.alpahs.mat[j][0] -= os.labelMat.mat[j][0] * (Ei - Ej) / eta;
os.alpahs.mat[j][0] = clipAlpha(os.alpahs.mat[j][0], H, L);
updateEk(j);
if (abs(os.alpahs.mat[j][0] - alphaJold) < 0.00001)
{
cout<<"j not moving enought"<<endl;
return 0;
}
os.alpahs.mat[i][0] += os.labelMat.mat[j][0] * os.labelMat.mat[i][0] *
(alphaJold - os.alpahs.mat[j][0]);
updateEk(i);
double b1 = os.b - Ei;
b1 -= os.labelMat.mat[i][0] * (os.alpahs.mat[i][0] - alphaIold) * os.k.mat[i][i];
b1 -= os.labelMat.mat[j][0] * (os.alpahs.mat[j][0] - alphaJold) * os.k.mat[i][j];
double b2 = os.b - Ej;
b2 -= os.labelMat.mat[i][0] * (os.alpahs.mat[i][0] - alphaIold) * os.k.mat[i][j];
b2 -= os.labelMat.mat[j][0] * (os.alpahs.mat[j][0] - alphaJold) * os.k.mat[j][j];
if ((0 < os.alpahs.mat[i][0]) && (os.C > os.alpahs.mat[i][0]))
{
os.b = b1;
}
else if ((0 < os.alpahs.mat[j][0]) && (os.C > os.alpahs.mat[j][0]))
{
os.b = b2;
}
else
os.b = (b1 + b2) / 2.0;
return 1;
}
else
return 0;
}
double smoP(Matrix dataMatIn, Matrix classLabels, Matrix &alphas, double C, double toler,int maxIter, KTup kTup)
{
}
};
int main()
{
srand((unsigned)time(NULL));
dataToMatrix dtm;
cout<<"load data"<<endl;
cout<<"------------"<<endl;
char file[20] = "G:/data/svm.txt";
dtm.loadData(&dtm,file);
Matrix dataMatrix;
dataMatrix.loadMatrix(&dataMatrix,dtm);
Matrix labelMatrix;
labelMatrix.initMatrix(&labelMatrix, dataMatrix.col,1);
labelMatrix = labelMatrix.getOneRow(dataMatrix,dataMatrix.row - 1);
dataMatrix.deleteOneRow(&dataMatrix, dataMatrix.row);
Matrix alphas;
alphas.initMatrix(&alphas, dataMatrix.col, 1, 0);
// cout<<smoSimple(dataMatrix, labelMatrix, alphas, 0.6, 0.001, 40)<<endl;
// SMOP smop;
// cout<<smop.smoP(dataMatrix, labelMatrix, alphas, 0.6, 0.001, 40)<<endl;
SMOPK smopk;
KTup kTup;
kTup.stringType = "rbf";
kTup.val = 1.3;
cout<<smopk.smoP(dataMatrix, labelMatrix, alphas, 0.6, 0.001,40, kTup)<<endl;
for (int i = 0; i < 100; i++)
{
if (alphas.mat[i][0] > 0.0)
{
for (int j = 0; j < dataMatrix.row; j++)
{
cout<<dataMatrix.mat[i][j]<<" ";
}
cout<<labelMatrix.mat[i][0]<<endl;
}
}
return 0;
}
第三步：运行过程
运行结果
用到的软件是vs2010以上的版本都可以，不用额外配置什么，没调包，会用这个软件进行c++开发，就会使用这个软件
此程序由于不调用任何外源库，所以读者可以看清楚每一个算法的原理，要想学好机器学习算法，必须打好基础，不要好高骛远，另外，程序都是有备注，应该很好理解的，实在不懂，可以来问店主
代码的下载路径（新窗口打开链接）：机器学习算法支持向量机SVM之c++实现（不调用外源库）
有问题可以私信或者留言，有问必答
54
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
