SVM支持向量机分类算法C++实现 零调库

代码中的alpha和有些教材中的lamda是同意义的值，在此记录说明

bits/stdc++.h是标准C和标准C++库的汇总库，不属于标准C/C++，但是省事。很多编译器是支持的，VS不支持

样本等均为随机生成

可以通过更改K函数修改核函数，此处为线性核

#include <bits/stdc++.h>
#define vd vector<double>
#define vvd vector<vd>
#define ull long long unsigned int
#define MY_RAND_MAX ((RAND_MAX<<15)|RAND_MAX)
using namespace std;

//#define DEBUG
//向量运算
template <class T>
T operator* (const vector<T>&AT,const vector<T>&B){
    T sum=0;
    for(ull i=0;i<AT.size()&&i<B.size();i++){
        sum+=AT[i]*B[i];
    }
    return sum;
}
template <class T>
vector<T> operator*(const vector<T>&A,T k){
    vector<T> kA(A);
    for(auto&x:kA)
        x*=k;
    return kA;
}

template <class T>
vector<T> operator*(T k,const vector<T>&A){
    vector<T> kA(A);
    for(auto&x:kA)
        x*=k;
    return kA;
}

//生成随机数
inline double my_rand(){
    return (rand()<<15)|rand();
}

inline double my_rand_double(){
    return my_rand() / double(MY_RAND_MAX);
}

inline double my_rand_double(double st,double ed){
    return my_rand_double()*(ed-st)+st;
}

//esp
const double esp=1e-7;
const double inf=DBL_MAX;
//α表 E表 K表
vd alpha,E;
vvd K;
//训练集
vvd x;
vd y;
//测试集
vvd x_test;
vd y_test;
//生成训练集样本的分类标准
double b;
vd omega;

//训练集样本容量与维度，测试集容量
int n,m;
int n_test;

//训练结果
double b_star=0;

//返回由alpha计算得到的omega_star
void omega_star(vd&res){
    res.clear();
    res.resize(m,0);
    for(int i=0;i<n;i++)if(alpha[i]>esp){
        for(int j=0;j<m;j++) {
            res[j]+=alpha[i]*y[i]*x[i][j];
        }
    }
}


//初始化K表
void init_K(){
    K.resize(n,vd(n));
    for(int i=0;i<n;i++)for(int j=0;j<n;j++){
        K[i][j]=x[i]*x[j];
    }
}

//g函数
double g(int i){
//E[i]=g(i)-y[i];
    return E[i]+y[i];
}

//更新E表
void set_E(){
    for(int i=0;i<n;i++){
        double res=b_star-y[i];
        for(int j=0;j<n;j++)if(alpha[j]>esp){
            if(abs(y[j]-1.0)<esp)
                res+=alpha[j]*K[i][j];
            else
                res-=alpha[j]*K[i][j];
        }
        E[i]=res;
    }
}

//一步更新
void step(int pos1,int pos2){
    double alpha_p1_new;
    double alpha_p2_new;
    double alpha_p1_old=alpha[pos1];
    double alpha_p2_old=alpha[pos2];

    double L,H;

    double K11=K[pos1][pos1];
    double K12=K[pos1][pos2];
    double K22=K[pos2][pos2];

    double y1=y[pos1];
    double y2=y[pos2];
    double E1=E[pos1];
    double E2=E[pos2];

    //计算L和H
    if(y1==y2){
        L=0;
        H=alpha_p1_old+alpha_p2_old;
    }
    else {
        L=max(0.0,alpha_p2_old-alpha_p1_old);
        H=inf;
    }

    //更新两个alpha
    alpha_p2_new=alpha_p2_old+y2*(E1-E2)/(K11+K22-K12-K12);

    if(alpha_p2_new>H)
        alpha_p2_new=H;
    else if(alpha_p2_new<L)
        alpha_p2_new=L;
    alpha_p1_new=alpha_p1_old+y1*y2*(alpha_p2_old-alpha_p2_new);

    alpha[pos1]=alpha_p1_new;
    alpha[pos2]=alpha_p2_new;

    //更新b_star
    double b_p1_new=-E1-y1*K11*(alpha_p1_new-alpha_p1_old)-y2*K12*(alpha_p2_new-alpha_p2_old)+b_star;
    double b_p2_new=-E2-y1*K12*(alpha_p1_new-alpha_p1_old)-y2*K22*(alpha_p2_new-alpha_p2_old)+b_star;
    b_star=(b_p1_new+b_p2_new)/2.0;

    //更新E
    set_E();

}

bool KKT(int i){
    if(abs(alpha[i])<=esp){
        if(y[i]==1)
            return g(i)>=1-esp;
        else
            return g(i)<=-1+esp;
    }
    else {
        return abs(y[i]-g(i))<=esp;
    }
}

bool test(int i,int j){
    return E[i]*E[j]<0;
    /*return abs(E[i]-E[j])>0.1;*/
}

void run(){
    bool flag=true;
    int ct=0;
    #ifdef DEBUG
    int loop = 0;
    while(flag){
        cout<<loop++<<endl;
    #else
    while(flag){
    #endif // DEBUG
        flag=false;
        for(int i=0;i<n;i++)if(!KKT(i)){
            for(int j=0;j<n;j++)if((!KKT(i)&&i!=j)&&((ct!=2&&test(i,j))||(alpha[j]>esp&&alpha[i]>esp))){
                flag=true;
                step(i,j);
#ifdef DEBUG
                cout<<"(i,j) = "<<'('<<i<<','<<j<<")\n";
                for(ull i=0;i<alpha.size();i++){
                    cout<<alpha[i]<<' ';
                }
                cout<<'\n';
                vd o_star;
                omega_star(o_star);
                cout<<o_star[0]<<' '<<o_star[1]<<' '<<b_star<<"\n";
                cout<<"E:\n";
                for(ull i=0;i<E.size();i++)
                    cout<<E[i]<<' ';
                cout<<"\n\n";
#endif // DEBUG
                int cntf0=0;
                for(int k=0; k<n; k++) if(alpha[k]>esp) cntf0++;
                if(cntf0>=n*0.2&&ct==0) ct=1;
                if(ct==1&&cntf0<=n*0.1) ct=2;

//                if(KKT(i))
//                    break;
            }
        }
    }
    #ifdef DEBUG
    cout<<"loop="<<loop<<endl;
    #endif
}

int main(){
#ifdef DEBUG
    //freopen("log.txt","w",stdout);
#endif
    srand(time(NULL));
    n=70;
    m=3;
    omega.resize(m);
    omega[0]=2;
    omega[1]=-0.5;
    omega[2]=-1.5;
    b=1;
    x.resize(n,vd(m));
    for(int i=0;i<n;i++){
        for(int j=0;j<m;j++)
            x[i][j]=my_rand_double(-10,10);
    }

    y.resize(n);
    for(int i=0;i<n;i++){
        double sum=0;
        for(int j=0;j<m;j++)
            sum+=x[i][j]*omega[j];
        if(sum+b>0)
            y[i]=1;
        else
            y[i]=-1;
    }

    init_K();
    alpha.resize(n);
    E.resize(n);

    set_E();

#ifdef DEBUG
    cout<<"train_set:"<<endl;
    for(int i=0;i<n;i++){
        for(int j=0;j<m;j++)
            cout<<x[i][j]<<' ';
        cout<<':';
        cout<<y[i]<<endl;
    }
    cout<<endl;
    cout<<"K:\n"<<endl;
    for(ull i=0;i<K.size();i++){
        for(ull j=0;j<K[i].size();j++)
            cout<<K[i][j]<<' ';
        cout<<endl;
    }

#endif //DEBUG





    run();

    cout<<"alpha:\n";
    for(ull i=0;i<alpha.size();i++){
        cout<<alpha[i]<<' ';
    }
    cout<<"\n\n";

    cout<<"ans:\n";
    vd o_star;
    omega_star(o_star);
    cout<<"ans_omega:\n";
    for(ull i=0;i<o_star.size();i++)
        cout<<o_star[i]<<' ';
    cout<<"\nb_star: "<<b_star<<"\n\n";



#ifdef DEBUG
    cout<<"E:\n";
    for(ull i=0;i<E.size();i++){
        cout<<E[i]<<' ';
    }
    cout<<"\n\n";

    cout<<"G:\n";
    for(int i=0;i<n;i++){
        cout<<g(i)<<' ';
    }
    cout<<"\n\n";
#endif //DEBUG

//测试集
    n_test=1000;
    x_test.resize(n_test,vd(m));
    for(int i=0;i<n_test;i++){
        for(int j=0;j<m;j++)
            x_test[i][j]=my_rand_double(-10,10);
    }
    y_test.resize(n_test);
    for(int i=0;i<n_test;i++){
        double sum=0;
        for(int j=0;j<m;j++)
            sum+=x_test[i][j]*omega[j];
        if(sum+b>0)
            y_test[i]=1;
        else
            y_test[i]=-1;
    }
    int t=0;
    int f=0;
    for(int i=0;i<n_test;i++){
        double sum=0;
        for(int j=0;j<m;j++)
            sum+=x_test[i][j]*o_star[j];
        double ans;
        if(sum+b_star>0)
            ans=1;
        else
            ans=-1;
        if(ans==y_test[i])
            t++;
        else {
            f++;
            cout<<"false: ";
            for(int j=0;j<m;j++)
                cout<<x_test[i][j]<<' ';
            cout<<endl;
        }
    }
    cout<<"number of true/false:"<<t<<"/"<<f<<endl;
}