svm.cpp

支持向量分类算法在linux操作系统下的是实现

//    Solver s;
//    s.Solve(l, SVC_Q(*prob,*param,y), minus_ones, y,
// 	   alpha, Cp, Cn, param->eps, si, param->shrinking);

  double sum_alpha=0;
  for(i=0;i<l;i++)
    sum_alpha += alpha[i];

  if (Cp==Cn)
    info("nu = %f\n", sum_alpha/(Cp*prob->l));

  for(i=0;i<l;i++)
    alpha[i] *= y[i];

  delete[] minus_ones;
  delete[] y;
}

static void solve_nu_svc(const svm_problem *prob, const svm_parameter *param,
			 double *alpha, Solver::SolutionInfo* si)
{
  int i;
  int l = prob->l;
  double nu = param->nu;

  schar *y = new schar[l];

  for(i=0;i<l;i++)
    if(prob->y[i]>0)
      {
	y[i] = +1;
      }
    else
      {
	y[i] = -1;
      }

  double sum_pos = nu*l/2;
  double sum_neg = nu*l/2;
  
  for(i=0;i<l;i++)
     if(y[i] == +1)
       {
 	 alpha[i] = min(1.0,sum_pos);
	 sum_pos -= alpha[i];
       }
     else
       {
 	 alpha[i] = min(1.0,sum_neg);
	 sum_neg -= alpha[i];
       }
  double *zeros = new double[l];
  
  for(i=0;i<l;i++)
    {
      zeros[i] = 0;
    }

   int size;
   MPI_Comm_size(MPI_COMM_WORLD, &size);

#ifdef SOLVER_PSMO
   Solver_Parallel_SMO_NU sps(param->o, param->q, MPI_COMM_WORLD);
   sps.Solve(l, SVC_Q(*prob,*param,y), zeros, y,
	     alpha, 1.0, 1.0, param->eps, si,  param->shrinking);
#endif
#ifdef SOLVER_LOQO
   Solver_Parallel_LOQO_NU spl(param->o, param->q, 2, MPI_COMM_WORLD, 
 			      size, 1, param->o/size, nu);
    spl.Solve(l, SVC_Q(*prob,*param,y), zeros, y,
	      alpha, 1.0, 1.0, param->eps, si,  param->shrinking);
#endif

    // Serial solvers:
    //   Solver_LOQO_NU sl(param->o, param->q, 2, param->nu);
    //   sl.Solve(l, SVC_Q(*prob,*param,y), zeros, y,
    //  	     alpha, 1.0, 1.0, param->eps, si,  param->shrinking);
    //   Solver_NU s;
    //   s.Solve(l, SVC_Q(*prob,*param,y), zeros, y,
    // 	  alpha, 1.0, 1.0, param->eps, si,  param->shrinking);
    //
    //   Solver_GPM_NU sgpm(param->o, param->q);
    //   sgpm.Solve(l, SVC_Q_LOQO(*prob,*param,y), zeros, y,
    // 	     alpha, 1.0, 1.0, param->eps, si,  param->shrinking);

  double r = si->r;
  printf("si->r %g\n",si->r);
  info("C = %f\n",1/r);

//   printf("alphan ");
  for(i=0;i<l;i++)
    {
      alpha[i] *= y[i]/r;
//       printf(" %g", alpha[i]);
    }
//   printf("\n");


  si->rho /= r;
  si->obj /= (r*r);
  si->upper_bound_p = 1/r;
  si->upper_bound_n = 1/r;

  printf("si->rho = %g", si->rho);
  printf("si->r = %g", si->r);

  delete[] y;
  delete[] zeros;
}

static void solve_one_class(const svm_problem *prob, const svm_parameter *param,
			    double *alpha, Solver::SolutionInfo* si)
{
  int l = prob->l;
  double *zeros = new double[l];
  schar *ones = new schar[l];
  int i;

  int n = (int)(param->nu*prob->l);	// # of alpha's at upper bound

  for(i=0;i<n;i++)
    alpha[i] = 1;
  if(n<prob->l)
    alpha[n] = param->nu * prob->l - n;
  for(i=n+1;i<l;i++)
    alpha[i] = 0;

  for(i=0;i<l;i++)
    {
      zeros[i] = 0;
      ones[i] = 1;
    }

   int size;
   MPI_Comm_size(MPI_COMM_WORLD, &size);

   // Note: The following has not been tested but
   // should work.
#ifdef SOLVER_PSMO
  Solver_Parallel_SMO sps(param->o, param->q, MPI_COMM_WORLD);
  sps.Solve(l, ONE_CLASS_Q(*prob,*param), zeros, ones,
	     alpha, 1.0, 1.0, param->eps, si,  param->shrinking);
#endif
#ifdef SOLVER_LOQO
   Solver_Parallel_LOQO spl(param->o, param->q, 2, MPI_COMM_WORLD, 
 			   size, 1, param->o/size);
   spl.Solve(l, ONE_CLASS_Q(*prob,*param), zeros, ones,
	     alpha, 1.0, 1.0, param->eps, si,  param->shrinking);
#endif

   // Serial solver:
   //   Solver s;
   //   s.Solve(l, ONE_CLASS_Q(*prob,*param), zeros, ones,
   //  	  alpha, 1.0, 1.0, param->eps, si, param->shrinking);

  delete[] zeros;
  delete[] ones;
}

static void solve_epsilon_svr(const svm_problem *prob, const svm_parameter *param,
			      double *alpha, Solver::SolutionInfo* si)
{
  int l = prob->l;
  double *alpha2 = new double[2*l];
  double *linear_term = new double[2*l];
  schar *y = new schar[2*l];
  int i;

  for(i=0;i<l;i++)
    {
      alpha2[i] = 0;
      linear_term[i] = param->p - prob->y[i];
      y[i] = 1;

      alpha2[i+l] = 0;
      linear_term[i+l] = param->p + prob->y[i];
      y[i+l] = -1;
    }

  int size;
  MPI_Comm_size(MPI_COMM_WORLD, &size);

#ifdef SOLVER_PSMO
  Solver_Parallel_SMO sps(param->o, param->q, MPI_COMM_WORLD);
  sps.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
	    alpha2, param->C, param->C, param->eps, si, param->shrinking);
#endif
#ifdef SOLVER_LOQO
   Solver_Parallel_LOQO spl(param->o, param->q, 1, MPI_COMM_WORLD, 
 			   size, 1, param->o/size);
   spl.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
	     alpha2, param->C, param->C, param->eps, si, param->shrinking);
#endif

   // Serial solver
   //   Solver s;
   //   s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
   // 	  alpha2, param->C, param->C, param->eps, si, param->shrinking);
  

  double sum_alpha = 0;
  for(i=0;i<l;i++)
    {
      alpha[i] = alpha2[i] - alpha2[i+l];
      sum_alpha += fabs(alpha[i]);
    }
  info("nu = %f\n",sum_alpha/(param->C*l));

  delete[] alpha2;
  delete[] linear_term;
  delete[] y;
}

static void solve_nu_svr(const svm_problem *prob, const svm_parameter *param,
			 double *alpha, Solver::SolutionInfo* si)
{
  int l = prob->l;
  double C = param->C;
  double *alpha2 = new double[2*l];
  double *linear_term = new double[2*l];
  schar *y = new schar[2*l];
  int i;

  double sum = C * param->nu * l / 2;
  for(i=0;i<l;i++)
    {
      alpha2[i] = alpha2[i+l] = min(sum,C);
      sum -= alpha2[i];

      linear_term[i] = - prob->y[i];
      y[i] = 1;

      linear_term[i+l] = prob->y[i];
      y[i+l] = -1;
    }
//   printf("alpha = ");
//   for(i = 0; i<2*l; ++i)
//     printf(" %g", alpha2[i]);
//   printf("\n");

  int size;
  MPI_Comm_size(MPI_COMM_WORLD, &size);

#ifdef SOLVER_PSMO
  Solver_Parallel_SMO_NU sps(param->o, param->q, MPI_COMM_WORLD);
  sps.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
	    alpha2, C, C, param->eps, si, param->shrinking);
#endif
#ifdef SOLVER_LOQO
  Solver_Parallel_LOQO_NU spl(param->o, param->q, 2, MPI_COMM_WORLD, 
 			      size, 1, param->o/size, C*param->nu/2);
  spl.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
	    alpha2, C, C, param->eps, si,  param->shrinking);
#endif

  // Serial solvers:
  //   Solver_LOQO_NU sl(param->o, param->q, 2, C*param->nu/2);
  //   sl.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
  // 	   alpha2, C, C, param->eps, si,  param->shrinking);
  //   Solver_NU s;
  //   s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
  // 	  alpha2, C, C, param->eps, si, param->shrinking);

  info("epsilon = %f\n",-si->r);

  for(i=0;i<l;i++)
    alpha[i] = alpha2[i] - alpha2[i+l];

  delete[] alpha2;
  delete[] linear_term;
  delete[] y;
}

//
// decision_function
//
struct decision_function
{
  double *alpha;
  double rho;	
};

decision_function svm_train_one(const svm_problem *prob, 
				const svm_parameter *param,
				double Cp, double Cn)
{
  double *alpha = Malloc(double,prob->l);
  Solver::SolutionInfo si;
  switch(param->svm_type)
    {
    case C_SVC:
      solve_c_svc(prob,param,alpha,&si,Cp,Cn);
      break;
    case NU_SVC:
      solve_nu_svc(prob,param,alpha,&si);
      break;
    case ONE_CLASS:
      solve_one_class(prob,param,alpha,&si);
      break;
    case EPSILON_SVR:
      solve_epsilon_svr(prob,param,alpha,&si);
      break;
    case NU_SVR:
      solve_nu_svr(prob,param,alpha,&si);
      break;
    }

  info("obj = %f, rho = %f\n",si.obj,si.rho);

  // output SVs

  int nSV = 0;
  int nBSV = 0;
  for(int i=0;i<prob->l;i++)
    {
      if(fabs(alpha[i]) > 0)
	{
	  ++nSV;
	  if(prob->y[i] > 0)
	    {
	      if(fabs(alpha[i]) >= si.upper_bound_p)
		++nBSV;
	    }
	  else
	    {
	      if(fabs(alpha[i]) >= si.upper_bound_n)
		++nBSV;
	    }
	}
    }

  info("nSV = %d, nBSV = %d\n",nSV,nBSV);

  decision_function f;
  f.alpha = alpha;
  f.rho = si.rho;
  return f;
}

//
// svm_model
//
struct svm_model
{
  svm_parameter param;	// parameter
  int nr_class;		// number of classes, = 2 in regression/one class svm
  int l;			// total #SV
  Xfloat **SV;
  int **nz_sv;
  int *sv_len;
  int max_idx;
  double **sv_coef;	// coefficients for SVs in decision functions (sv_coef[n-1][l])
  double *rho;		// constants in decision functions (rho[n*(n-1)/2])
  double *probA;          // pariwise probability information
  double *probB;

  // for classification only

  int *label;		// label of each class (label[n])
  int *nSV;		// number of SVs for each class (nSV[n])
  // nSV[0] + nSV[1] + ... + nSV[n-1] = l
  // XXX
  int free_sv;		// 1 if svm_model is created by svm_load_model
  // 0 if svm_model is created by svm_train
};

// Platt's binary SVM Probablistic Output: an improvement from Lin et al.
void sigmoid_train(int l, const double *dec_values, const double *labels, 
		   double& A, double& B)
{
  double prior1=0, prior0 = 0;
  int i;

  for (i=0;i<l;i++)
    if (labels[i] > 0) prior1+=1;
    else prior0+=1;
	
  int max_iter=100; 	// Maximal number of iterations
  double min_step=1e-10;	// Minimal step taken in line search
  double sigma=1e-3;	// For numerically strict PD of Hessian
  double eps=1e-5;
  double hiTarget=(prior1+1.0)/(prior1+2.0);
  double loTarget=1/(prior0+2.0);
  double *t=Malloc(double,l);
  double fApB,p,q,h11,h22,h21,g1,g2,det,dA,dB,gd,stepsize;
  double newA,newB,newf,d1,d2;
  int iter; 
	
  // Initial Point and Initial Fun Value
  A=0.0; B=log((prior0+1.0)/(prior1+1.0));
  double fval = 0.0;

  for (i=0;i<l;i++)
    {
      if (labels[i]>0) t[i]=hiTarget;
      else t[i]=loTarget;
      fApB = dec_values[i]*A+B;
      if (fApB>=0)
	fval += t[i]*fApB + log(1+exp(-fApB));
      else
	fval += (t[i] - 1)*fApB +log(1+exp(fApB));
    }
  for (iter=0;iter<max_iter;iter++)
    {
      // Update Gradient and Hessian (use H' = H + sigma I)
      h11=sigma; // numerically ensures strict PD
      h22=sigma;
      h21=0.0;g1=0.0;g2=0.0;
      for (i=0;i<l;i++)
	{
	  fApB = dec_values[i]*A+B;
	  if (fApB >= 0)
	    {
	      p=exp(-fApB)/(1.0+exp(-fApB));
	      q=1.0/(1.0+exp(-fApB));
	    }
	  else
	    {
	      p=1.0/(1.0+exp(fApB));
	      q=exp(fApB)/(1.0+exp(fApB));
	    }
	  d2=p*q;
	  h11+=dec_values[i]*dec_values[i]*d2;
	  h22+=d2;
	  h21+=dec_values[i]*d2;
	  d1=t[i]-p;
	  g1+=dec_values[i]*d1;
	  g2+=d1;
	}

      // Stopping Criteria
      if (fabs(g1)<eps && fabs(g2)<eps)
	break;

      // Finding Newton direction: -inv(H') * g
      det=h11*h22-h21*h21;
      dA=-(h22*g1 - h21 * g2) / det;
      dB=-(-h21*g1+ h11 * g2) / det;
      gd=g1*dA+g2*dB;


      stepsize = 1; 		// Line Search
      while (stepsize >= min_step)
	{
	  newA = A + stepsize * dA;
	  newB = B + stepsize * dB;

	  // New function value
	  newf = 0.0;
	  for (i=0;i<l;i++)
	    {
	      fApB = dec_values[i]*newA+newB;
	      if (fApB >= 0)
		newf += t[i]*fApB + log(1+exp(-fApB));
	      else
		newf += (t[i] - 1)*fApB +log(1+exp(fApB));
	    }
	  // Check sufficient decrease
	  if (newf<fval+0.0001*stepsize*gd)
	    {
	      A=newA;B=newB;fval=newf;
	      break;
	    }
	  else
	    stepsize = stepsize / 2.0;
	}

      if (stepsize < min_step)
	{
	  info("Line search fails in two-class probability estimates\n");
	  break;
	}
    }

  if (iter>=max_iter)
    info("Reaching maximal iterations in two-class probability estimates\n");
  free(t);
}

double sigmoid_predict(double decision_value, double A, double B)
{
  double fApB = decision_value*A+B;
  if (fApB >= 0)
    return exp(-fApB)/(1.0+exp(-fApB));
  else
    return 1.0/(1+exp(fApB)) ;
}

// Method 2 from the multiclass_prob paper by Wu, Lin, and Weng
void multiclass_probability(int k, double **r, double *p)
{
  int t,j;
  int iter = 0, max_iter=100;
  double **Q=Malloc(double *,k);
  double *Qp=Malloc(double,k);
  double pQp, eps=0.005/k;
	
  for (t=0;t<k;t++)
    {
      p[t]=1.0/k;  // Valid if k = 1
      Q[t]=Malloc(double,k);
      Q[t][t]=0;
      for (j=0;j<t;j++)
	{
	  Q[t][t]+=r[j][t]*r[j][t];
	  Q[t][j]=Q[j][t];
	}
      for (j=t+1;j<k;j++)
	{
	  Q[t][t]+=r[j][t]*r[j][t];
	  Q[t][j]=-r[j][t]*r[t][j];
	}
    }
  for (iter=0;iter<max_iter;iter++)
    {
      // stopping condition, recalculate QP,pQP for numerical accuracy
      pQp=0;
      for (t=0;t<k;t++)
	{
	  Qp[t]=0;
	  for (j=0;j<k;j++)
	    Qp[t]+=Q[t][j]*p[j];
	  pQp+=p[t]*Qp[t];
	}
      double max_error=0;
      for (t=0;t<k;t++)
	{
	  double error=fabs(Qp[t]-pQp);
	  if (error>max_error)
	    max_error=error;
	}
      if (max_error<eps) break;
		
      for (t=0;t<k;t++)
	{
	  double diff=(-Qp[t]+pQp)/Q[t][t];
	  p[t]+=diff;
	  pQp=(pQp+diff*(diff*Q[t][t]+2*Qp[t]))/(1+diff)/(1+diff);
	  for (j=0;j<k;j++)
	    {
	      Qp[j]=(Qp[j]+diff*Q[t][j])/(1+diff);
	      p[j]/=(1+diff);
	    }
	}
    }
  if (iter>=max_iter)
    info("Exceeds max_iter in multiclass_prob\n");
  for(t=0;t<k;t++) free(Q[t]);
  free(Q);
  free(Qp);
}

// Cross-validation decision values for probability estimates
void svm_binary_svc_probability(const svm_problem *prob, 
				const svm_parameter *param,
				double Cp, double Cn, 
				double& probA, double& probB)
{
  int i;
  int nr_fold = 5;
  int *perm = Malloc(int,prob->l);
  double *dec_values = Malloc(double,prob->l);

  // random shuffle
  for(i=0;i<prob->l;i++) perm[i]=i;
  for(i=0;i<prob->l;i++)
    {
      int j = i+rand()%(prob->l-i);
      swap(perm[i],perm[j]);
    }
  for(i=0;i<nr_fold;i++)
    {
      int begin = i*prob->l/nr_fold;
      int end = (i+1)*prob->l/nr_fold;
      int j,k;
      struct svm_problem subprob;

      subprob.l = prob->l-(end-begin);
      subprob.max_idx = prob->max_idx;
      subprob.x = Malloc(Xfloat *,subprob.l);
      subprob.nz_idx = Malloc(int *,subprob.l);
      subprob.x_len = Malloc(int,subprob.l);
      subprob.y = Malloc(double,subprob.l);
			
      k=0;
      for(j=0;j<begin;j++)
	{
	  subprob.x[k] = prob->x[perm[j]];
	  subprob.nz_idx[k] = prob->nz_idx[perm[j]];
	  subprob.x_len[k] = prob->x_len[perm[j]];
	  subprob.y[k] = prob->y[perm[j]];
	  ++k;
	}
      for(j=end;j<prob->l;j++)
	{