📄 smo_cpp.htm
字号:
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<!-- saved from url=(0050)http://bbs.matwav.com/user/download/205363/smo.cpp -->
<HTML><HEAD>
<META http-equiv=Content-Type content="text/html; charset=gb2312">
<META content="MSHTML 6.00.2800.1276" name=GENERATOR></HEAD>
<BODY><PRE>#include "iostream.h"
#include "fstream.h"
#include "ctype.h"
#include "math.h"
#include "stdlib.h"
#include "time.h"
#define N 14104//样本数总数
#define end_support_i 12090
#define first_test_i 12090
#define d 340//维数
//////////global variables
//int N= //样本数
//int d= //维数
float C=100;//惩罚参数
float tolerance=0.001;//ui与边界0,C之间的公差范围
float eps=1e-3;//一个近似0的小数
float two_sigma_squared=10;//RBF核函数中的参数
float alph[end_support_i];//Lagrange multipiers
float b;//threshold
float error_cache[end_support_i];//存放non-bound样本误差
int target[N];//训练与测试样本的目标值
float precomputed_self_dot_product[N];//预存dot_product_func(i,i)的值，以减少计算量
int dense_points[N][d];//存放训练与测试样本，0-end_support_i-1训练;first_test_i-N测试
//函数的申明
int takeStep(int,int);
int examineNonBound(int);
int examineBound(int);
int examineFirstChoice(int,float);
int examineExample(int);
float kernel_func(int,int);
float learned_func(int);
float dot_product_func(int,int);
float error_rate();
void setX();
void setT();
void initialize();
///////////////////////////////////////////////////////
void main()
{
ofstream os("d:\\smo1.txt");//存放训练的后的参数
int numChanged=0,examineAll=1;
//srand((unsigned int)time(NULL));
initialize();
//以下两层循环，开始时检查所有样本，选择不符合KKT条件的两个乘子进行优化，选择成功，返回1，否则，返回0
//所以成功了，numChanged必然>0，从第二遍循环时，就不从整个样本中去寻找不符合KKT条件的两个乘子进行优化，
//而是从边界的乘子中去寻找，因为非边界样本需要调整的可能性更大，边界样本往往不被调整而始终停留在边界上。
//如果，没找到，再从整个样本中去找，直到整个样本中再也找不到需要改变的乘子为止，此时，算法结束。
while(numChanged>0||examineAll)
{
numChanged=0;
if(examineAll)
{
for(int k=0;k<end_support_i;k++)
numChanged+=examineExample(k);//examin all exmples
}
else{
for(int k=0;k<end_support_i;k++)
if(alph[k]!=0&&alph[k]!=C)
numChanged+=examineExample(k);//loop k over all non-bound lagrange multipliers
}
if(examineAll==1)
examineAll=0;
else if(numChanged==0)
examineAll=1;
}
//存放训练后的参数
{
os<<"d="<<d<<endl;
os<<"b="<<b<<endl;
os<<"two_sigma_squared="<<two_sigma_squared<<endl;
int n_support_vectors=0;
for(int i=0;i<end_support_i;i++)
if(alph[i]>0)//此地方是否可以改为alph[i]>tolerance?????????????????
n_support_vectors++;
os<<"n_support_vectors="<<n_support_vectors<<"rate="<<(float)n_support_vectors/first_test_i<<endl;
for(i=0;i<end_support_i;i++)
if(alph[i]>0)
os<<"alph["<<i<<"]="<<alph[i]<<" ";
os<<endl;
}
//输出，以及测试
cout<<error_rate()<<endl;
}
/////////examineExample程序
//假定第一个乘子ai(位置为i1)，examineExample(i1)首先检查，如果它超出tolerance而违背KKT条件，那么它就成为第一个乘子；
//然后，寻找第二个乘子(位置为i2),通过调用takeStep(i1,i2)来优化这两个乘子。
int examineExample(int i1)
{
float y1,alph1,E1,r1;
y1=target[i1];
alph1=alph[i1];
if(alph1>0&&alph1<C)
E1=error_cache[i1];
else
E1=learned_func(i1)-y1;//learned_func为计算输出函数
r1=y1*E1;
//////违反KKT条件的判断
if((r1>tolerance&&alph1>0)||(r1<-tolerance&&alph1<C))
{
/////////////使用三种方法选择第二个乘子
//1：在non-bound乘子中寻找maximum fabs(E1-E2)的样本
//2：如果上面没取得进展,那么从随机位置查找non-boundary 样本
//3：如果上面也失败，则从随机位置查找整个样本,改为bound样本
if (examineFirstChoice(i1,E1))//1
{
return 1;
}
if (examineNonBound(i1))//2
{
return 1;
}
if (examineBound(i1))//3
{
return 1;
}
}
///没有进展
return 0;
}
////////////1：在non-bound乘子中寻找maximum fabs(E1-E2)的样本
int examineFirstChoice(int i1,float E1)
{
int k,i2;
float tmax;
for(i2=-1,tmax=0.0,k=0;k<end_support_i;k++)//*******************************end_support_i
if(alph[k]>0&&alph[k]<C)
{
float E2,temp;
E2=error_cache[k];
temp=fabs(E1-E2);
if(temp>tmax)
{
tmax=temp;
i2=k;
}
}
if(i2>=0)
{
if(takeStep(i1,i2))
return 1;
}
return 0;
}
/////////////// 2：如果上面没取得进展,那么从随机位置查找non-boundary 样本
int examineNonBound(int i1)
{
int k,k0 = rand()%end_support_i;
int i2;
for (k = 0; k < end_support_i; k++)
{
i2 = (k + k0) % end_support_i;//从随机位开始
if (alph[i2] > 0.0 && alph[i2] < C)
{
if (takeStep(i1, i2))
{
return 1;
}
}
}
return 0;
}
////////////3：如果上面也失败，则从随机位置查找整个样本,(改为bound样本)
int examineBound(int i1)
{
int k,k0 = rand()%end_support_i;
int i2;
for (k = 0; k < end_support_i; k++)
{
i2 = (k + k0) % end_support_i;//从随机位开始
//if (alph[i2]= 0.0 || alph[i2]=C)//修改******************
{
if (takeStep(i1, i2))
{
return 1;
}
}
}
return 0;
}
///////////takeStep()
//用于优化两个乘子，成功，返回1，否则，返回0
int takeStep(int i1,int i2)
{
int y1,y2,s;
float alph1,alph2;//两个乘子的旧值
float a1,a2;//两个乘子的新值
float E1,E2,L,H,k11,k22,k12,eta,Lobj,Hobj,delta_b;
if(i1==i2) return 0;//不会优化两个同一样本
//给变量赋值
alph1=alph[i1];
alph2=alph[i2];
y1=target[i1];
y2=target[i2];
if(alph1>0&&alph1<C)
E1=error_cache[i1];
else
E1=learned_func(i1)-y1;//learned_func(int)为非线性的评价函数，即输出函数
if(alph2>0&&alph2<C)
E2=error_cache[i2];
else
E2=learned_func(i2)-y2;
s=y1*y2;
//计算乘子的上下限
if(y1==y2)
{
float gamma=alph1+alph2;
if(gamma>C)
{
L=gamma-C;
H=C;
}
else
{
L=0;
H=gamma;
}
}
else
{
float gamma=alph1-alph2;
if(gamma>0)
{
L=0;
H=C-gamma;
}
else
{
L=-gamma;
H=C;
}
}//计算乘子的上下限
if(L==H) return 0;
//计算eta
k11=kernel_func(i1,i1);//kernel_func(int,int)为核函数
k22=kernel_func(i2,i2);
k12=kernel_func(i1,i2);
eta=2*k12-k11-k22;
if(eta<0)
{
a2=alph2+y2*(E2-E1)/eta;//计算新的alph2
//调整a2，使其处于可行域
if(a2<L) a2=L;
if(a2>H) a2=H;
}
else//此时，得分别从端点H,L求目标函数值Lobj,Hobj，然后设a2为求得最大目标函数值的端点值
{
float c1=eta/2;
float c2=y2*(E2-E1)-eta*alph2;
Lobj=c1*L*L+c2*L;
Hobj=c1*H*H+c2*H;
if(Lobj>Hobj+eps)//eps****************
a2=L;
else if(Lobj<Hobj-eps)
a2=H;
else
a2=alph2;//加eps的目的在于，使得Lobj与Hobj尽量分开，如果，很接近，就认为没有改进(make progress)
}
if(fabs(a2-alph2)<eps*(a2+alph2+eps))
return 0;
a1=alph1-s*(a2-alph2);//计算新的a1
if(a1<0)//调整a1,使其符合条件*****??????????????????????????????????????????
{
a2+=s*a1;
a1=0;
}
else if(a1>C)
{
float t=a1-C;
a2+=s*t;
a1=C;
}
//更新阀值b
{
float b1,b2,bnew;
if(a1>0&&a1<C)
bnew=b+E1+y1*(a1-alph1)*k11+y2*(a2-alph2)*k12;
else{
if(a2>0&&a2<C)
bnew=b+E2+y1*(a1-alph1)*k12+y2*(a2-alph2)*k22;
else{
b1=b+E1+y1*(a1-alph1)*k11+y2*(a2-alph2)*k12;
b2=b+E2+y1*(a1-alph1)*k12+y2*(a2-alph2)*k22;
bnew=(b1+b2)/2;
}
}
delta_b=bnew-b;
b=bnew;
}
//对于线性情况，要更新权向量，这里不用了
//更新error_cache，对取得进展的a1,a2,所对应的i1,i2的error_cache[i1]=error_cache[i2]=0
{
float t1=y1*(a1-alph1);
float t2=y2*(a2-alph2);
for(int i=0;i<end_support_i;i++)
if(0<alph[i]&&alph[i]<C)
error_cache[i]+=t1*kernel_func(i1,i)+t2*(kernel_func(i2,i))-delta_b;
error_cache[i1]=0.0;
error_cache[i2]=0.0;
}
alph[i1]=a1;//store a1,a2 in the alpha array
alph[i2]=a2;
return 1;//说明已经取得进展
}
//////////learned_func(int)
//评价分类学习函数
float learned_func(int k)
{
float s=0.0;
for(int i=0;i<end_support_i;i++)
if(alph[i]>0)
s+=alph[i]*target[i]*kernel_func(i,k);
s-=b;
return s;
}
//计算点积函数dot_product_func(int,int)
float dot_product_func(int i1,int i2)
{
float dot=0;
for(int i=0;i<d;i++)
dot+=dense_points[i1][i]*dense_points[i2][i];
return dot;
}
//径向机核函数RBF:kernel_func(int,int)
float kernel_func(int i1,int i2)
{
float s=dot_product_func(i1,i2);
s*=-2;
s+=precomputed_self_dot_product[i1]+precomputed_self_dot_product[i2];
return exp(-s/two_sigma_squared);
}
//初始化initialize()
void initialize()
{
//初始化阀值b为0
b=0.0;
//初始化alph[]为0
for(int i=0;i<end_support_i;i++)
alph[i]=0.0;
//设置样本值矩阵
setX();
//设置目标值向量
setT();
//设置预计算点积
{for(int i=0;i<N;i++)
precomputed_self_dot_product[i]=dot_product_func(i,i);//注意，这是对训练样本的设置，对于测试样本也应考虑?????????????????
}
}
//计算误差率error_rate()
float error_rate()
{ ofstream to("d:\\smo_test.txt");
int tp=0,tn=0,fp=0,fn=0;
float ming=0,te=0,total_q=0,temp=0;
for(int i=first_test_i;i<N;i++)
{ temp=learned_func(i);
if(temp>0&&target[i]>0)
tp++;
else if(temp>0&&target[i]<0)
fp++;
else if(temp<0&&target[i]>0)
fn++;
else if(temp<0&&target[i]<0)
tn++;
to<<i<<" 实际输出"<<temp<<endl;
}
total_q=(float)(tp+tn)/(float)(tp+tn+fp+fn);//总精度
ming=(float)tp/(float)(tp+fn);
te=(float)tp/(float)(tp+fp);
to<<"---------------测试结果-----------------"<<endl;
to<<"tp="<<tp<<" tn="<<tn<<" fp="<<fp<<" fn="<<fn<<endl;
to<<"ming="<<ming<<" te="<<te<<" total_q="<<total_q<<endl;
return (1-total_q);
}
//计算样本X[]
void setX(){
ifstream ff("D:\\17_smo.txt");////////////////////////
if(!ff)
{
cerr<<"error!!"<<endl;
exit(1);
}
int i=0,j=0;
char ch;
while(ff.get(ch))
{
if(isspace(ch)) continue;
if(isdigit(ch))
{
dense_points[i][j]=(int)ch-48;
j++;
if(j==d)
{
j=0;
i++;
}
}
}
}
//set targetT[]
void setT()
{
for(int i=0;i<6045;i++)
target[i]=1;
for(i=6045;i<12090;i++)
target[i]=-1;
for(i=12090;i<13097;i++)
target[i]=1;
for(i=13097;i<14104;i++)
target[i]=-1;
}</PRE></BODY></HTML>
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -