loqo_solver.cpp

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

  printf("c=");
  for(int i=0; i<n; ++i)
    printf(" %g", c[i]);
  printf("\n");
  printf("low=");
  for(int i=0; i<n; ++i)
    printf(" %g", low[i]);
  printf("\n");
  printf("up=");
  for(int i=0; i<n; ++i)
    printf(" %g", up[i]);
  printf("\n");
}

double Solver_LOQO_NU::calculate_rho()
{
  int nr_free1 = 0,nr_free2 = 0;
  double ub1 = INF, ub2 = INF;
  double lb1 = -INF, lb2 = -INF;
  double sum_free1 = 0, sum_free2 = 0;

  for(int i=0;i<active_size;i++)
    {
      if(y[i]==+1)
	{
	  if(is_lower_bound(i))
	    ub1 = min(ub1,G[i]);
	  else if(is_upper_bound(i))
	    lb1 = max(lb1,G[i]);
	  else
	    {
	      ++nr_free1;
	      sum_free1 += G[i];
	    }
	}
      else
	{
	  if(is_lower_bound(i))
	    ub2 = min(ub2,G[i]);
	  else if(is_upper_bound(i))
	    lb2 = max(lb2,G[i]);
	  else
	    {
	      ++nr_free2;
	      sum_free2 += G[i];
	    }
	}
    }
  printf("nr_free1 = %d\n", nr_free1);
  printf("sum_free1 = %g\n",sum_free1);
  printf("nr_free2 = %d\n", nr_free2);
  printf("sum_freee = %g\n",sum_free2);
  double r1,r2;
  if(nr_free1 > 0)
    r1 = sum_free1/nr_free1;
  else
    r1 = (ub1+lb1)/2;
	
  if(nr_free2 > 0)
    r2 = sum_free2/nr_free2;
  else
    r2 = (ub2+lb2)/2;
	
  si->r = (r1+r2)/2;
  printf("(r1+r2)/2 = %g\n", (r1+r2)/2);
  printf("(r1+r2)/2 = %g\n", (r1-r2)/2);
  return (r1-r2)/2;
}

void Solver_LOQO::init_working_set(int *work_set, int *not_work_set)
{
  int j;
  NEXT_RAND = 1;
  for (int i=0; i<n; ++i)
  {
    do
      {
	j = next_rand_pos() % l;
      } while (work_status[j] != WORK_N);
    work_status[j] = WORK_B;
  }
  int k=0; j=0;
  for(int i=0; i<l; ++i)
    {
      if(work_status[i] == WORK_B)
	{
	  work_set[j] = i; ++j;
	  work_count[i] = 0;
	}
      else
	{
	  not_work_set[k] = i; ++k;
	  work_count[i] = -1;
	}
    }
}

int Solver_LOQO::solve_inner()
{
  Solver_NU sl;
  schar *sy = new schar[n];
  Qfloat *Q_bb_f = new Qfloat[n*n];
  Qfloat *QD_f = new Qfloat[n];
  for(int i=0; i<n; ++i)
    {
      sy[i] = (schar) a[i];
      QD_f[i] = (Qfloat) Q_bb[i*n+i];
      for(int j=0; j<n; ++j)
	Q_bb_f[i*n+j] = (Qfloat) Q_bb[i*n+j];
    }

  sl.Solve(n, SVQ_No_Cache(Q_bb_f, QD_f, n), c, sy, work_space,
	   Cp, Cn, eps, si, /* shrinking */ 0);
  delete[] sy;
  delete[] QD_f;
  return 0;
}

// int Solver_LOQO::solve_inner()
// {
//   int result = 0;
//   double sigdig;
//   double epsilon_loqo = 1e-10;
//   double margin;
//   int iteration;
//   for(margin=init_margin, iteration=init_iter; 
//       margin <= 0.9999999 && result != OPTIMAL_SOLUTION;)
//     {
//       sigdig = -log10(opt_precision);
//       // run pr_loqo
//       result = pr_loqo(n, m, c, Q_bb, a, d, low, up, work_space,
// 		       &work_space[3*n], dist, 3, sigdig, iteration, margin, 
// 		       up[0]/4, 0);
//       // Note: No need to check for choldc problem, as we employ
//       // Manteuffel shifting.
//       if(result == -1)
// 	{
// 	  info("NOTICE: Restarting PR_LOQO with more conservative");
// 	  info("parameters.\n"); info_flush();
// 	  if(init_margin < 0.80)
// 	    init_margin = (4.0*margin+1.0)/5.0;
// 	  margin = (margin+1.0)/2.0;
// 	  opt_precision *= 10.0;
// 	  info("NOTICE: Reducing precision of PR_LOQO.\n");
// 	}
//       else if(result != OPTIMAL_SOLUTION)
// 	{
// 	  // increase number of iterations
// 	  iteration += 2000;
// 	  init_iter += 10;
// 	  // reduce precision
// 	  opt_precision *= 10.0;
// 	  info("NOTICE: Reducing precision of PR_LOQO!\n"); info_flush();
// 	}
//     }
//   // Check precision of alphas
//   for(int i=0; i<n; ++i)
//     {
//       if(work_space[i] < up[i]-epsilon_loqo && dist[i] < -eps)
// 	epsilon_loqo = 2*(up[i]-work_space[i]);
//       else if(work_space[i] > low[i]+epsilon_loqo && dist[i] > eps)
// 	epsilon_loqo = 2*work_space[i];
//     }
//   info("Using epsilon_loqo= %g\n", epsilon_loqo); info_flush();
//   // Clip alphas to bounds
//   for(int i=0; i<n; ++i)
//     {
//       if(fabs(work_space[i]) < epsilon_loqo)
// 	{
// 	  work_space[i] = 0;
// 	}
//       if(fabs(work_space[i]-get_C(i)) < epsilon_loqo)
// 	{
// 	  work_space[i] = get_C(i);
// 	}
//     }
//   // Compute obj after optimization
//   double obj_after = 0.0;
//   for(int i=0; i<n; ++i)
//     {
//       obj_after += work_space[i]*c[i];
//       obj_after += 0.5*work_space[i]*work_space[i]*Q_bb[i*n+i];
//       for(int j=0; j<i; ++j)
// 	{
// 	  obj_after += work_space[j]*work_space[i]*Q_bb[j*n+i];
// 	}
//     }
//   printf("obj_after/before = %g / %g\n", obj_after, obj_before); fflush(stdout);
//   printf("delta a/b = %g\n", obj_after-obj_before); fflush(stdout);
//   // Check for progress
//   if(obj_after >= obj_before)
//     {
//       // Increase precision
//       opt_precision /= 100.0;
//       ++precision_violations;
//       info("NOTICE: Increasing Precision of PR_LOQO.\n"); info_flush();
//     }
//   if(precision_violations > 500)
//     {
//       // Relax stopping criterion
//       eps *= 10.0;
//       precision_violations = 0;
//       info("WARNING: Relaxing epsilon on KKT-Conditions.\n"); info_flush();
//     }
//   return result;
// }

void Solver_LOQO::Solve(int l, const QMatrix& Q, const double *b_, 
			const schar *y_, double *alpha_, double Cp, 
			double Cn, double eps, SolutionInfo* si, 
			int shrinking)
{
  this->l = l;
  this->si = si;
  this->Q = &Q;
  QD = Q.get_QD(); // we need diagonal for working_set_selection
  clone(b,b_,l);
  clone(y,y_,l);
  clone(alpha,alpha_,l);
  this->Cp = Cp;
  this->Cn = Cn;
  this->eps = eps;
  this->active_size = l;
  this->lmn = l - n;
  int *work_set = new int[n];
  int *not_work_set = new int[l];
  double *delta_alpha = new double[n];
  // init alpha status and work status
  {
    alpha_status = new char[l];
    work_status = new char[l];
    work_count = new int[l];
    for(int i=0; i<l; ++i)
      {
	update_alpha_status(i);
	work_status[i] = WORK_N;
	work_count[i] = -1;
      }
  }
  // init gradient
  info("initializing gradient..."); info_flush();
  {
    G = new double[l];
    for(int i=0; i<l; ++i)
      {
	G[i] = b[i];
      }
    for(int i=0; i<l; ++i)
      {
 	if(!is_lower_bound(i))
 	  {
	    const Qfloat *Q_i = Q.get_Q(i,l); 
	    double alpha_i = alpha[i];
	    for(int j=0; j<l; ++j)
	      {
		G[j] += alpha_i * Q_i[j];
	      }
 	  }
      }
  }
  info("done.\n"); info_flush();

  // Allocate space for pr_loqo
  Q_bb = new LOQOfloat[n*n];
  c = new LOQOfloat[n];
  up = new LOQOfloat[n];
  low = new LOQOfloat[n];
  dist = new LOQOfloat[n];
  allocate_a();
  allocate_d();
  allocate_work_space();
  iter=0;
  // optimization loop
  while(1)
    {
      // select working set and check for optimality
      if(iter > 0)
	{
	  if(select_working_set(work_set, not_work_set) != 0)
	    break;
	}
      else
	{
	  info("initializing working set..."); info_flush();
   	  init_working_set(work_set, not_work_set);
	  info("done.\n"); info_flush();
 	}

      lmn = l - n;
//       printf("working_set=");
//       for(int i=0; i<n; ++i)
// 	{
// 	  printf("%d ",work_set[i]);
// 	}
//       printf("\n");
//       printf("not_working_set=");
//       for(int i=0; i<lmn; ++i)
// 	{
// 	  printf("%d ",not_work_set[i]);
// 	}
//       printf("\n");

      ++iter;
      // setup problem for pr_loqo
      info("setting up problem for pr_loqo..."); info_flush();
      setup_problem(work_set, not_work_set);
      setup_up(work_set);
      setup_low();
//       print_problem();
      info("done.\n");
      // run inner solver
      for(int i=0; i<n; ++i)
	work_space[i]=alpha[work_set[i]];

      // Compute obj before optimization
      obj_before = 0.0;
      for(int i=0; i<n; ++i)
	{
	  obj_before += alpha[work_set[i]]*c[i];
	  obj_before += 0.5*alpha[work_set[i]]*alpha[work_set[i]]*Q_bb[i*n+i];
	  for(int j=0; j<i; ++j)
	    {
	      obj_before += alpha[work_set[j]]*alpha[work_set[i]]*Q_bb[j*n+i];
	    }
	}
      printf("obj_before = %g\n", obj_before); fflush(stdout);

      int status = solve_inner();
      printf("pr_loqo status = %d\n",status);

      // Restore Q_bb, lower triangle, overwritten by pr_loqo
      for(int i=0; i<n; ++i)
	{
	  for(int j=i+1; j<n; ++j)
	    {
	      Q_bb[n*j+i] = Q_bb[n*i+j];
	    }
	}
      // update gradient, compute G_b += Q_bb*delta_alpha and 
      // G_n += Q_nb*delta_alpha.
      int *nz = new int[n];
      double sum_delta_alpha_pos=0;
      double sum_delta_alpha_neg=0;
      for(int i=0; i<n; ++i)
	{
	  delta_alpha[i] = work_space[i] - alpha[work_set[i]];
	  if(y[work_set[i]]==+1)
	    sum_delta_alpha_pos += delta_alpha[i];
	  else
	    sum_delta_alpha_neg += delta_alpha[i];
   	  if(fabs(delta_alpha[i]) > TOL_ZERO)
	    nz[i] = 1;
  	  else
  	    nz[i] = 0;
	}
      printf("sum_delta_alpha_pos = %g\n", sum_delta_alpha_pos);
      printf("sum_delta_alpha_neg = %g\n", sum_delta_alpha_neg);
      for(int i=0; i<n; ++i)
	for(int j=0; j<n; ++j) // G_b
	  {
	    if(nz[j])
	      G[work_set[i]] += Q_bb[n*i+j]*delta_alpha[j];
	  }

      for(int i=0; i<n; ++i) // G_n
	{
	  if(nz[i])
	    {
 	      const Qfloat *Q_i = 
 		Q.get_Q_subset(work_set[i],not_work_set,lmn);
//  	      const Qfloat *Q_i = Q.get_Q(work_set[i],l);
	      for(int j=0; j<lmn; ++j)
		G[not_work_set[j]] += Q_i[not_work_set[j]] * delta_alpha[i];
	    }
	}
      delete[] nz;
      // update alpha
      for(int i=0; i<n; ++i)
	{
 	  alpha[work_set[i]] = work_space[i];
	  update_alpha_status(work_set[i]);
	}
//       double sum_alpha=0;
//       printf("alpha = ");
//       double sum_alpha_pos=0;
//       double sum_alpha_neg=0;
//       for(int i=0; i<l; ++i)
// 	{
// 	  sum_alpha += alpha[i];
// 	  if(y[i]==+1)
// 	    sum_alpha_pos += alpha[i];
// 	  else
// 	    sum_alpha_neg += alpha[i];
// 	  printf(" %g", alpha[i]);
// 	}
//       printf("\n");
//       printf("sum_alpha = %g\n",sum_alpha);
//       printf("sum_alpha_pos = %g\n",sum_alpha_pos);
//       printf("sum_alpha_neg = %g\n",sum_alpha_neg);
//       double delta_alpha_pos = (sum_alpha/2)-sum_alpha_pos;
//       double delta_alpha_neg = (sum_alpha/2)-sum_alpha_neg;
      // Restore consistency
//       if(fabs(delta_alpha_pos) > 0)
// 	{
// 	  for(int i=0; i<l; ++i)
// 	    {
// 	      if(y[i] == +1)
// 		{
// 		  if(alpha[i]+delta_alpha_pos >= low[i] &&
// 		     alpha[i]+delta_alpha_pos <= up[i])
// 		    {
// 		      alpha[i] += delta_alpha_pos;
// 		      break;
// 		    }
// 		}
// 	    }
// 	}
//       if(fabs(delta_alpha_neg) > 0)
// 	{
// 	  for(int i=0; i<l; ++i)
// 	    {
// 	      if(y[i] == -1)
// 		{
// 		  if(alpha[i]+delta_alpha_neg >= low[i] &&
// 		     alpha[i]+delta_alpha_neg <= up[i])
// 		    {
// 		      alpha[i] += delta_alpha_neg;
// 		      break;
// 		    }
// 		}
// 	    }
// 	}
      // Recheck
//       sum_alpha_pos = 0;
//       sum_alpha_neg = 0;
//       sum_alpha = 0;
//       for(int i=0; i<l; ++i)
// 	{
// 	  sum_alpha += alpha[i];
// 	  if(y[i] == +1)
// 	    sum_alpha_pos += alpha[i];
// 	  else
// 	    sum_alpha_neg += alpha[i];
// 	}
//       printf("sum_alpha = %g\n",sum_alpha);
//       printf("sum_alpha_pos = %g\n",sum_alpha_pos);
//       printf("sum_alpha_neg = %g\n",sum_alpha_neg);
    } // while(1)

  // calculate rho
  si->rho=calculate_rho();

  // calculate objective value
  {
    double v = 0;
    int i;
    for(i=0;i<l;i++)
      v += alpha[i] * (G[i] + b[i]);

    si->obj = v/2;
  }

  // put back the solution
  {
    for(int i=0;i<l;i++)
      alpha_[i] = alpha[i];
  }

  si->upper_bound_p = Cp;
  si->upper_bound_n = Cn;

  info("\noptimization finished, #iter = %d\n",iter);

  // clean up
  delete[] b;
  delete[] y;
  delete[] alpha;
  delete[] alpha_status;
  delete[] delta_alpha;
  delete[] work_status;
  delete[] work_count;
  delete[] work_set;
  delete[] not_work_set;
  delete[] G;
  delete[] Q_bb;
  delete[] c;
  delete[] up;
  delete[] low;
  delete[] dist;
  delete[] a;
  delete[] d;
  delete[] work_space;
}