理论部分见:http://www.cnblogs.com/biyeymyhjob/archive/2012/07/17/2591592.html
http://blog.youkuaiyun.com/sunmenggmail/article/details/7445035
//这个版本的代码仅可以算是初步实现,一些代码优化还有停止条件的选择还有待改进。
#include <iostream>
#include <cmath>
using namespace std;
#define N 12 //样本空间大小
#define C 4 //惩罚参数
#define M 2 //样本维数
double b = 4;//原始问题常数参数b
double W[M]; //原始问题参数
double A[N] = {0}; //初始化对偶问题参数
//训练样本
double X[N][M] = {
{1, 2},
{1, 1},
{-2, 2},
{2, 1},
{-2, 1},
{-1, 1},
{1, 4},
{2, 3},
{3, 2},
{-1, 4},
{-2, 3},
{-3, 2}
};
double Y[] = {-1, -1, -1,-1,-1,-1, 1, 1, 1, 1, 1, 1};
//kernel tricks
double Kernel (double *x, double *y, int length)
{
double sum = 0;
for (int k = 0; k < M; k++)
sum += pow(x[k] * y[k] ,2);//将x=(x1, x2)映射到z = (z1, z2) = (x1*x1, x2*x2)上,将x空间的椭圆变换成新空间里的直线,当然,这只是实验,核函数还有很多,具体选什么样的核函数还要视情况而定。
//return pow(sum, 2);
return sum;
}
//kernel tricks for samples
double K (int i, int j)
{
return Kernel(X[i], X[j], M);
}
//tool function
double Convert(double *x)
{
int j = 0;
double sum = 0;
for (j = 0; j < N; j++)
{
sum += A[j] * Y[j] * Kernel(X[j], x, M);
}
return sum;
}
//求Ei
double E(int i)
{
return Convert(X[i]) + b - Y[i];
}
//求函数间隔,与Ei仅差一个Yi
double g(double *x)
{
return Convert(x) + b;
}
//选择两个变量进行坐标下降
bool choose(int &a1, int &a2)
{
int i;
for (i = 0; i < N; ++i)
{
//选择违反KKT的a作为第一个的变量
if ((A[i] == 0 && Y[i] * g(X[i]) < 1) || (A[i] < C && A[i] > 0 && Y[i] * g(X[i]) != 1) || (A[i] == C && Y[i]*g(X[i]) > 1))
{
a1 = i;
break;
}
else if( i == N - 1)
{
return true; //当所有参数都满足KKT时停止,无法满足时?TODO:将可行间隙比率 增加为停止条件
}
}
double e1 = E(a1);
double e, e2, max = 0;
for (i = 0; i < N; i++)
{
if (i != a1)
{
e = E(i);
if (abs(e1-e) > max)
{
max = abs(e1 - e);
a2 = i;
e2 = e;
}
}
}
double L, H;
if (Y[a1] != Y[a2])
{
L = A[a2] - A[a1] > 0 ? A[a2] - A[a1] : 0;
H = C + A[a2] -A[a1] < C ? C + A[a2] -A[a1] :C;
}
else
{
L = A[a2] + A[a1] - C > 0 ? A[a2] + A[a1] - C : 0;
H = A[a2] + A[a1] < C ? A[a2] + A[a1] :C;
}
double temp = A[a2] + Y[a2] * (e1 - e2) / (K(a1, a1) + K(a2, a2) -2 * K (a1, a2));
double old = A[a2];
A[a2] = temp < L ? L : (temp > H ? H : temp); //更新a2
A[a1] = A[a1] + Y[a1]* Y[a2] * (old - A[a2]);//更新a1
return false;
}
double compute_b(int num, int a1, int a2)
{
double sum = 0;
int t = num == 1 ? a1: a2;
for (int i = 0; i < N; i++)
{
if (i != a1 && i!= a2)
sum += A[i] * Y[i] * K(i, t);
}
return Y[t] - sum - A[a1] * Y[a1] * K(a1, t) - A[a2] * Y[a2] * K(a2, t);
}
//SMO 算法解最优化对偶问题
void SMO()
{
int a1, a2;
double b1, b2;
int i =0;
while(true)
{
i++;
if (i > 100)
break;//当停止条件改进后,就不再需要
if (choose(a1, a2))
{
//满足KKT就停止
return;
}
else//更新b
{
b1 = compute_b(1, a1, a2);
b2 = compute_b(2, a1, a2);
if ((A[a1] < C && A[a1] > 0) && (A[a2] < C && A[a2]) > 0)
{
b = b1 ;
}
if ((A[a1] == 0 || A[a1] == C) && (A[a2] == 0 || A[a2] == C) )
{
b = (b1 + b2) / 2;
}
}
}
}
int main()
{
SMO();
double test[] = {0, -5};
cout << g(test) << endl;
return 0;
}
实验结果为:2.94999 ,说明(0,-5)点是正的,其他结果同样可以验证。
SVM分类结果图: