svm.cpp

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

	  // reset active set size and check
	  active_size = l;
	  // info("*"); info_flush();
	  if(select_working_set(i,j)!=0)
	    break;
	  else
	    counter = 1;	// do shrinking next iteration
	}
		
      ++iter;

      // update alpha[i] and alpha[j], handle bounds carefully
		
      const Qfloat *Q_i = Q.get_Q(i,active_size);
      const Qfloat *Q_j = Q.get_Q(j,active_size);

      double C_i = get_C(i);
      double C_j = get_C(j);

      double old_alpha_i = alpha[i];
      double old_alpha_j = alpha[j];

      if(y[i]!=y[j])
	{
	  double quad_coef = Q_i[i]+Q_j[j]+2*Q_i[j];
	  if (quad_coef <= 0)
	    quad_coef = TAU;
	  double delta = (-G[i]-G[j])/quad_coef;
	  double diff = alpha[i] - alpha[j];
	  alpha[i] += delta;
	  alpha[j] += delta;
			
	  if(diff > 0)
	    {
	      if(alpha[j] < 0)
		{
		  alpha[j] = 0;
		  alpha[i] = diff;
		}
	    }
	  else
	    {
	      if(alpha[i] < 0)
		{
		  alpha[i] = 0;
		  alpha[j] = -diff;
		}
	    }
	  if(diff > C_i - C_j)
	    {
	      if(alpha[i] > C_i)
		{
		  alpha[i] = C_i;
		  alpha[j] = C_i - diff;
		}
	    }
	  else
	    {
	      if(alpha[j] > C_j)
		{
		  alpha[j] = C_j;
		  alpha[i] = C_j + diff;
		}
	    }
	}
      else
	{
	  double quad_coef = Q_i[i]+Q_j[j]-2*Q_i[j];
	  if (quad_coef <= 0)
	    quad_coef = TAU;
	  double delta = (G[i]-G[j])/quad_coef;
	  double sum = alpha[i] + alpha[j];
	  alpha[i] -= delta;
	  alpha[j] += delta;

	  if(sum > C_i)
	    {
	      if(alpha[i] > C_i)
		{
		  alpha[i] = C_i;
		  alpha[j] = sum - C_i;
		}
	    }
	  else
	    {
	      if(alpha[j] < 0)
		{
		  alpha[j] = 0;
		  alpha[i] = sum;
		}
	    }
	  if(sum > C_j)
	    {
	      if(alpha[j] > C_j)
		{
		  alpha[j] = C_j;
		  alpha[i] = sum - C_j;
		}
	    }
	  else
	    {
	      if(alpha[i] < 0)
		{
		  alpha[i] = 0;
		  alpha[j] = sum;
		}
	    }
	}
//       printf("alpha=");
//       for(int k=0; k<active_size; ++k)
// 	printf("%g ", alpha[k]);
//       printf("\n");

      // update G

      double delta_alpha_i = alpha[i] - old_alpha_i;
      double delta_alpha_j = alpha[j] - old_alpha_j;
		
      for(int k=0;k<active_size;k++)
	{
	  G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;
	}

      // update alpha_status and G_bar

      {
	bool ui = is_upper_bound(i);
	bool uj = is_upper_bound(j);
	update_alpha_status(i);
	update_alpha_status(j);
	int k;
	if(ui != is_upper_bound(i))
	  {
	    Q_i = Q.get_Q(i,l);
	    if(ui)
	      for(k=0;k<l;k++)
		G_bar[k] -= C_i * Q_i[k];
	    else
	      for(k=0;k<l;k++)
		G_bar[k] += C_i * Q_i[k];
	  }

	if(uj != is_upper_bound(j))
	  {
	    Q_j = Q.get_Q(j,l);
	    if(uj)
	      for(k=0;k<l;k++)
		G_bar[k] -= C_j * Q_j[k];
	    else
	      for(k=0;k<l;k++)
		G_bar[k] += C_j * Q_j[k];
	  }
      }
    }

  // 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
  {
//     info("alpha = ");
    for(int i=0;i<l;i++)
      {
	alpha_[active_set[i]] = alpha[i];
// 	info(" %g",alpha[i]);
      }
//     info("\n"); info_flush();
  }

  // juggle everything back
  /*{
    for(int i=0;i<l;i++)
    while(active_set[i] != i)
    swap_index(i,active_set[i]);
    // or Q.swap_index(i,active_set[i]);
    }*/

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

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

  delete[] b;
  delete[] y;
  delete[] alpha;
  delete[] alpha_status;
  delete[] active_set;
  delete[] G;
  delete[] G_bar;
}

// return 1 if already optimal, return 0 otherwise
int Solver::select_working_set(int &out_i, int &out_j)
{
  // return i,j such that
  // i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)
  // j: minimizes the decrease of obj value
  //    (if quadratic coefficeint <= 0, replace it with tau)
  //    -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)
	
  double Gmax = -INF;
  int Gmax_idx = -1;
  int Gmin_idx = -1;
  double obj_diff_min = INF;

//   printf("G = ");
  for(int t=0;t<active_size;t++)
    {
      if(y[t]==+1)	
	{
	  if(!is_upper_bound(t))
	    if(-G[t] > Gmax)
	      {
		Gmax = -G[t];
		Gmax_idx = t;
	      }
	}
      else
	{
	  if(!is_lower_bound(t))
	    if(G[t] > Gmax)
	      {
		Gmax = G[t];
		Gmax_idx = t;
	      }
	}
//       printf("%g ",G[t]);
    }
//   printf("\nGmax = %g\n", Gmax);
  int i = Gmax_idx;
  const Qfloat *Q_i = NULL;
  if(i != -1) // NULL Q_i not accessed: Gmax=-INF if i=-1
    Q_i = Q->get_Q(i,active_size);
//   printf("quad_coef=");
  for(int j=0;j<active_size;j++)
    {
      if(y[j]==+1)
	{
	  if (!is_lower_bound(j))
	    {
	      double grad_diff=Gmax+G[j];
	      if (grad_diff >= eps)
		{
		  double obj_diff; 
		  double quad_coef=Q_i[i]+QD[j]-2*y[i]*Q_i[j];
		  //	  printf("%g ", quad_coef);
		  if (quad_coef > 0)
		    obj_diff = -(grad_diff*grad_diff)/quad_coef;
		  else
		    obj_diff = -(grad_diff*grad_diff)/TAU;

		  if (obj_diff <= obj_diff_min)
		    {
		      Gmin_idx=j;
		      obj_diff_min = obj_diff;
		    }
		}
	    }
	}
      else
	{
	  if (!is_upper_bound(j))
	    {
	      double grad_diff= Gmax-G[j];
	      if (grad_diff >= eps)
		{
		  double obj_diff; 
		  double quad_coef=Q_i[i]+QD[j]+2*y[i]*Q_i[j];
// 		  printf("%g ", quad_coef);
		  if (quad_coef > 0)
		    obj_diff = -(grad_diff*grad_diff)/quad_coef;
		  else
		    obj_diff = -(grad_diff*grad_diff)/TAU;

		  if (obj_diff <= obj_diff_min)
		    {
		      Gmin_idx=j;
		      obj_diff_min = obj_diff;
		    }
		}
	    }
	}
    }
//   printf("\n");
  if(Gmin_idx == -1)
    return 1;
//   printf("i=%d j=%d\n",Gmax_idx, Gmin_idx);

  out_i = Gmax_idx;
  out_j = Gmin_idx;
  return 0;
}

// return 1 if already optimal, return 0 otherwise
int Solver::max_violating_pair(int &out_i, int &out_j)
{
  // return i,j: maximal violating pair

  double Gmax1 = -INF;		// max { -y_i * grad(f)_i | i in I_up(\alpha) }
  int Gmax1_idx = -1;

  double Gmax2 = -INF;		// max { y_i * grad(f)_i | i in I_low(\alpha) }
  int Gmax2_idx = -1;

  for(int i=0;i<active_size;i++)
    {
      if(y[i]==+1)	// y = +1
	{
	  if(!is_upper_bound(i))	// d = +1
	    {
	      if(-G[i] >= Gmax1)
		{
		  Gmax1 = -G[i];
		  Gmax1_idx = i;
		}
	    }
	  if(!is_lower_bound(i))	// d = -1
	    {
	      if(G[i] >= Gmax2)
		{
		  Gmax2 = G[i];
		  Gmax2_idx = i;
		}
	    }
	}
      else		// y = -1
	{
	  if(!is_upper_bound(i))	// d = +1
	    {
	      if(-G[i] >= Gmax2)
		{
		  Gmax2 = -G[i];
		  Gmax2_idx = i;
		}
	    }
	  if(!is_lower_bound(i))	// d = -1
	    {
	      if(G[i] >= Gmax1)
		{
		  Gmax1 = G[i];
		  Gmax1_idx = i;
		}
	    }
	}
    }
  //printf("Gmax1+Gmax2 = %g, i=%d, j=%d\n", Gmax1+Gmax2, Gmax1_idx, Gmax2_idx);
  if(Gmax1+Gmax2 < eps)
    return 1;

  out_i = Gmax1_idx;
  out_j = Gmax2_idx;
  return 0;
}

void Solver::do_shrinking()
{
  int i,j,k;
  if(max_violating_pair(i,j)!=0) return;
  double Gm1 = -y[j]*G[j];
  double Gm2 = y[i]*G[i];

  // shrink
	
  for(k=0;k<active_size;k++)
    {
      if(is_lower_bound(k))
	{
	  if(y[k]==+1)
	    {
	      if(-G[k] >= Gm1) continue;
	    }
	  else	if(-G[k] >= Gm2) continue;
	}
      else if(is_upper_bound(k))
	{
	  if(y[k]==+1)
	    {
	      if(G[k] >= Gm2) continue;
	    }
	  else	if(G[k] >= Gm1) continue;
	}
      else continue;

      --active_size;
      swap_index(k,active_size);
      --k;	// look at the newcomer
    }

  // unshrink, check all variables again before final iterations

  if(unshrinked || -(Gm1 + Gm2) > eps*10) return;
	
  unshrinked = true;
  reconstruct_gradient();

  for(k=l-1;k>=active_size;k--)
    {
      if(is_lower_bound(k))
	{
	  if(y[k]==+1)
	    {
	      if(-G[k] < Gm1) continue;
	    }
	  else	if(-G[k] < Gm2) continue;
	}
      else if(is_upper_bound(k))
	{
	  if(y[k]==+1)
	    {
	      if(G[k] < Gm2) continue;
	    }
	  else	if(G[k] < Gm1) continue;
	}
      else continue;

      swap_index(k,active_size);
      active_size++;
      ++k;	// look at the newcomer
    }
}

double Solver::calculate_rho()
{
  double r;
  int nr_free = 0;
  double ub = INF, lb = -INF, sum_free = 0;
  for(int i=0;i<active_size;i++)
    {
      double yG = y[i]*G[i];

      if(is_lower_bound(i))
	{
	  if(y[i] > 0)
	    ub = min(ub,yG);
	  else
	    lb = max(lb,yG);
	}
      else if(is_upper_bound(i))
	{
	  if(y[i] < 0)
	    ub = min(ub,yG);
	  else
	    lb = max(lb,yG);
	}
      else
	{
	  ++nr_free;
	  sum_free += yG;
	}
    }

  if(nr_free>0)
    r = sum_free/nr_free;
  else
    r = (ub+lb)/2;
  return r;
}

//
// Solver using LOQO, added by D. Brugger 15.08.06
// $Id: svm.cpp 1 2006-12-16 14:03:21Z beeblbrox $

class SVQ_No_Cache : public QMatrix
{
public:
  SVQ_No_Cache(Qfloat *Q_bb_, Qfloat *QD_, const int n_);
  Qfloat *get_Q(int i, int len) const;
  Qfloat *get_Q_subset(int i, int *idxs, int n) const;
  Qfloat *get_QD() const;
  Qfloat get_non_cached(int i, int j) const;
  inline bool is_cached(const int i) const;
  void swap_index(int i, int j) const;
  virtual ~SVQ_No_Cache();
private:
  Qfloat *Q_bb;
  Qfloat *QD;
  int *idx;
  int n;
};

SVQ_No_Cache::SVQ_No_Cache(Qfloat *Q_bb_, Qfloat *QD_, const int n_)
{
  this->Q_bb = Q_bb_;
  this->QD = QD_;
  this->n = n_;
  this->idx = new int[n];
  for(int i=0; i<n; ++i)
    idx[i] = i;
}

Qfloat *SVQ_No_Cache::get_Q(int i, int len) const
{
  return &Q_bb[n*idx[i]];
}

Qfloat *SVQ_No_Cache::get_Q_subset(int i, int *idxs, int n) const
{
  printf("Error: Not implemented yet!\n"); exit(1);
  return NULL;
}

Qfloat *SVQ_No_Cache::get_QD() const
{
  return QD;
}

Qfloat SVQ_No_Cache::get_non_cached(int i, int j) const
{
  return Q_bb[i*n+j];
}
inline bool SVQ_No_Cache::is_cached(const int i) const
{
  return false;
}

void SVQ_No_Cache::swap_index(int i, int j) const
{
  swap(idx[i], idx[j]);
  swap(QD[i], QD[j]);
}

SVQ_No_Cache::~SVQ_No_Cache()
{
  delete[] idx;
}

#ifdef SOLVER_PSMO
#include "psmo/psmo_solver.cpp"
#endif
//#include "gpm/gpm_solver.cpp"
#ifdef SOLVER_LOQO
#include "loqo/loqo_solver.cpp"
#endif

//
// Q matrices for various formulations
//
class SVC_Q: public Kernel
{ 
public:
  SVC_Q(const svm_problem& prob, const svm_parameter& param, const schar *y_)
    :Kernel(prob.l, prob.x, prob.nz_idx, prob.x_len, prob.max_idx, param)
  {
    clone(y,y_,prob.l);
    this->l = prob.l;
    cache = new Cache(prob.l,(int)(param.cache_size*(1<<20)));
    QD = new Qfloat[prob.l];
    for(int i=0;i<prob.l;i++)
      QD[i]= (Qfloat)(this->*kernel_function)(i,i);
  }
	
  Qfloat *get_Q(int i, int len) const
  {
    Qfloat *data;
    int start;
    if((start = cache->get_data(i,&data,len)) < len)
      {
	for(int j=start;j<len;j++)
	  data[j] = (Qfloat)(y[i]*y[j]*(this->*kernel_function)(i,j));
      }
    return data;
  }

  Qfloat *get_QD() const
  {
    return QD;
  }

  Qfloat *get_Q_subset(int i, int *idxs, int n) const
  {
    Qfloat *data;
    int start = cache->get_data(i,&data,l);
    if(start == 0) // Initialize cache row
      {
	for(int j=0; j<l; ++j)
	  data[j] = NAN;
      }
    for(int j=0; j<n; ++j)
      {
	if(isnan(data[idxs[j]]))
	  data[idxs[j]] = (Qfloat)(y[i]*y[idxs[j]]*
				   (this->*kernel_function)(i,idxs[j]));
      }