支持向量机（SVM）入门C#代码实现smo（四，算法编译调试成功，针对开头举例，数学推导与程序运行结果一致）

最新推荐文章于 2023-12-27 08:16:53 发布

工业机器视觉设计和实现

最新推荐文章于 2023-12-27 08:16:53 发布

阅读量537

点赞数

分类专栏：人工智能浅尝支持向量机svm

本文链接：https://blog.csdn.net/ganggangwawa/article/details/108510083

版权

人工智能浅尝同时被 2 个专栏收录

75 篇文章 7 订阅

订阅专栏

支持向量机svm

6 篇文章 1 订阅

订阅专栏

注释：前三篇参考豆子的博客和支持向量机三重境界

int RandomSelect(int i)
{
int j = i;
while (i == j)
{
Random rm = new Random();
j = rm.Next(0, 2);
}
return j;
}

int kk(int i, int j)
{
int k11 = 18; int k12 = 21; int k21 = k12; int k13 = 6; int k31 = 6; int k23 = 7; int k32 = 7; int k22 = 25; int k33 = 2;
int temp = 0;
switch (i)
{
case 0://1
{
if (j == 0) {temp=k11; } else if (j == 1) {temp=k12; } else { temp=k13;}
}
break;
case 1: { if (j == 0) {temp=k21; } else if (j == 1) {temp=k22; } else {temp=k23; } }//2
break;
case 2: { if (j == 0) { temp = k31; } else if (j ==1) { temp = k32; } else { temp = k33; } }//3
break;
}
return temp;
}
void smo()
{
List<Point> xx=new List<Point>();
Point temppt=new Point(3,3);
xx.Add(temppt);
temppt=new Point(4,3);
xx.Add(temppt);
temppt=new Point(1,1);
xx.Add(temppt);

List<int> yi= new List<int>();
int tempyi = 1;
yi.Add(tempyi);
tempyi = 1;
yi.Add(tempyi);
tempyi = -1;
yi.Add(tempyi);
/*
       * 默认输入参数值
       * C: regularization parameter
       * tol: numerical tolerance
       * max passes
       */
       double C = 1; //对不在界内的惩罚因子
       double tol = 0.01;//容忍极限值
       int maxPasses = 5; //表示没有改变拉格朗日乘子的最多迭代次数
List<int> boundAlpha = new List<int>();
       /*
       * 初始化a[], b, passes
       */

double[] a = new double[xx.Count];//拉格朗日乘子
       //this.a = a;

       //将乘子初始化为0
for (int i = 0; i < xx.Count; i++)
{
           a[i] = 0;
       }
       int passes = 0;
double b = 0;

       while (passes < maxPasses) {
           //表示改变乘子的次数（基本上是成对改变的）
           int num_changed_alphas = 0;
for (int i = 0; i < xx.Count; i++)
{
               //表示特定阶段由a和b所决定的输出与真实yi的误差
               //参照公式(7)
               //double Ei = getE(i);
int k11 = 18; int k12 = 21; int k21 = k12; int k13 = 6; int k31 = 6; int k23 = 7; int k32 = 7; int k22 = 25; int k33 = 2;
double fxi = 0;
switch (i)
{
case 0://1
{ fxi = a[0] * yi[0] * k11 + a[1] * yi[1] * k21 + a[2] * yi[2] * k31 + b; }
break;
case 1://2
{ fxi = a[0] * yi[0] * k12 + a[1] * yi[1] * k22 + a[2] * yi[2] * k32 + b; }
break;
case 2://3
{ fxi = a[0] * yi[0] * k13 + a[1] * yi[1] * k23 + a[2] * yi[2] * k33 + b; }
break;
}
double Ei =fxi-yi[i];
               /*;
               * 把违背KKT条件的ai作为第一个
               * 满足KKT条件的情况是：
               * yi*f(i) >= 1 and alpha == 0 (正确分类)
               * yi*f(i) == 1 and 0<alpha < C (在边界上的支持向量)
               * yi*f(i) <= 1 and alpha == C (在边界之间)
               *
               *
               *
               * ri = y[i] * Ei = y[i] * f(i) - y[i]^2 >= 0
               * 如果ri < 0并且alpha < C 则违反了KKT条件
               * 因为原本ri < 0 应该对应的是alpha = C
               * 同理，ri > 0并且alpha > 0则违反了KKT条件
               * 因为原本ri > 0对应的应该是alpha =0
               */
               if ((yi[i] * Ei < -tol && a[i] < C) ||
                   (yi[i] * Ei > tol && a[i] > 0))
               {
                   /*
                   * ui*yi=1边界上的点 0 < a[i] < C
                   * 找MAX|E1 - E2|
                   */
                   int j;
                   /*
                   * boundAlpha表示x点处于边界上所对应的
                   * 拉格朗日乘子a的集合
                   */
if (boundAlpha.Count > 0)
{
//参照公式(5)
int maxK = -1; //用于保存临时最大索引
double maxDeltaE = 0; //用于保存临时最大差值--->|Ei-Ej|
//j = findMax(Ei, boundAlpha);
for (int ii = 0; ii < boundAlpha.Count; ii++)
{
int k = boundAlpha[ii];
if (k == i) continue;
double fxk = 0;
switch (k)
{
case 0://1
{ fxk = a[0] * yi[0] * k11 + a[1] * yi[1] * k21 + a[2] * yi[2] * k31 + b; }
break;
case 1://2
{ fxk = a[0] * yi[0] * k12 + a[1] * yi[1] * k22 + a[2] * yi[2] * k32 + b; }
break;
case 2://3
{ fxk = a[0] * yi[0] * k13 + a[1] * yi[1] * k23 + a[2] * yi[2] * k33 + b; }
break;
}
double Ek = fxk - yi[k];
double deltaE = Math.Abs(Ei - Ek);
if (deltaE > maxDeltaE)
{
maxK = k;
maxDeltaE = deltaE;
}
}
j = maxK;//0--2
}
else
{
//如果边界上没有，就随便选一个j != i的aj
j = RandomSelect(i);//1--3,与0---2矛盾，这个函数已更正未0--2,已经全部更改为0---2
}

                   //double Ej = getE(j);
double fxj = 0;
switch (j)
{
case 0://1
{ fxj = a[0] * yi[0] * k11 + a[1] * yi[1] * k21 + a[2] * yi[2] * k31 + b; }
break;
case 1://2
{ fxj = a[0] * yi[0] * k12 + a[1] * yi[1] * k22 + a[2] * yi[2] * k32 + b; }
break;
case 2://3
{ fxj = a[0] * yi[0] * k13 + a[1] * yi[1] * k23 + a[2] * yi[2] * k33 + b; }
break;
}
double Ej = fxj - yi[j];
                   //保存当前的ai和aj
                   double oldAi = a[i];
                   double oldAj = a[j];

                   /*
                   * 计算乘子的范围U, V
                   * 参考公式(4)
                   */
                   double L, H;
                   if (yi[i] != yi[j]) {
                       L = Math.Max(0, a[j] - a[i]);
                       H = Math.Min(C, C - a[i] + a[j]);
                   } else {
                       L = Math.Max(0, a[i] + a[j] - C);
                       H = Math.Min(0, a[i] + a[j]);
                   }


                   /*
                   * 如果eta等于0或者大于0 则表明a最优值应该在L或者U上
                   */
                   double eta = 2 * kk(i, j) - kk(i, i) - kk(j, j);//公式(3)

                   if (eta >= 0)
                       continue;

                   a[j] = a[j] - yi[j] * (Ei - Ej)/ eta;//公式(2)
                   if (0 < a[j] && a[j] < C)
                       boundAlpha.Add(j);

                   if (a[j] < L)
                       a[j] = L;
                   else if (a[j] > H)
                       a[j] = H;
if (L == H) continue;

                   if (Math.Abs(a[j] - oldAj) < 1e-5)
                       continue;
                   a[i] = a[i] + yi[i] * yi[j] * (oldAj - a[j]);
                   if (0 < a[i] && a[i] < C)
                       boundAlpha.Add(i);

                   /*
                   * 计算b1， b2
                   * 参照公式(6)
                   */
                   double b1 = b - Ei - yi[i] * (a[i] - oldAi) * kk(i, i) - yi[j] * (a[j] - oldAj) * kk(i, j);
                   double b2 = b - Ej - yi[i] * (a[i] - oldAi) * kk(i, j) - yi[j] * (a[j] - oldAj) * kk(j, j);

                   if (0 < a[i] && a[i] < C)
                       b = b1;
                   else if (0 < a[j] && a[j] < C)
                       b = b2;
                   else
                       b = (b1 + b2) / 2;

                   num_changed_alphas = num_changed_alphas + 1;
double min = a[0] * a[0] * yi[0] * yi[0] * k11 + a[0] * a[1] * yi[0] * yi[1] * k12 + a[0] * a[2] * yi[0] * yi[2] * k13
+ a[1] * a[0] * yi[1] * yi[0] * k21 + a[1] * a[1] * yi[1] * yi[1] * k22 + a[1] * a[2] * yi[1] * yi[2] * k23
+ a[2] * a[0] * yi[2] * yi[0] * k31 + a[2] * a[1] * yi[2] * yi[1] * k32 + a[2] * a[2] * yi[2] * yi[2] * k33
- a[0] - a[1] - a[2];
               }
           }
           if (num_changed_alphas == 0) {
               passes++;
           } else
               passes = 0;
       }

   //   return new SVMModel(a, y, b);//返回a【】是支持向量，b是截距，w（i）=sum（1，m）a【i】y（i）*x（i），w(i)*x+b=0是分割超平面