svm.cpp

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

//#include "psmo_solver.h"
//#include "loqo_solver.h"
#include "util.h"

#include <stdarg.h>
#include <time.h>
#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include <mpi.h>
#include "svm.h"

#define INF HUGE_VAL
#define TAU 1e-12
#define Malloc(type,n) (type *)malloc((n)*sizeof(type))

inline double powi(double base, int times)
{
  double tmp = base, ret = 1.0;
  for(int t=times; t>0; t/=2)
    {
      if(t%2==1)
	ret *= tmp;
      tmp = tmp * tmp;
    }
  return ret;
}
inline double clocks2sec(clock_t t)
{
  return (double) t / CLOCKS_PER_SEC;
}
#if 1
void info(char *fmt,...)
{
  va_list ap;
  va_start(ap,fmt);
  vprintf(fmt,ap);
  va_end(ap);
}
void info_flush()
{
  fflush(stdout);
}
#else
void info(char *fmt,...) {}
void info_flush() {}
#endif

//
// Kernel Cache
//
// l is the number of total data items
// size is the cache size limit in bytes
//
class Cache
{
public:
  Cache(int l,int size);
  ~Cache();

  // request data [0,len)
  // return some position p where [p,len) need to be filled
  // (p >= len if nothing needs to be filled)
  int get_data(const int index, Qfloat **data, int len);
  void swap_index(int i, int j);	// future_option
  // check if row is in cache
  bool is_cached(const int index) const;
private:
  int l;
  int size;
  struct head_t
  {
    head_t *prev, *next;	// a cicular list
    Qfloat *data;
    int len;		// data[0,len) is cached in this entry
  };

  head_t *head;
  head_t lru_head;
  void lru_delete(head_t *h);
  void lru_insert(head_t *h);
};

Cache::Cache(int l_,int size_):l(l_),size(size_)
{
  head = (head_t *)calloc(l,sizeof(head_t));	// initialized to 0
  size /= sizeof(Qfloat);
  size -= l * sizeof(head_t) / sizeof(Qfloat);
  size = max(size, 2*l);	// cache must be large enough for two columns
  lru_head.next = lru_head.prev = &lru_head;
}

Cache::~Cache()
{
  for(head_t *h = lru_head.next; h != &lru_head; h=h->next)
    free(h->data);
  free(head);
}

void Cache::lru_delete(head_t *h)
{
  // delete from current location
  h->prev->next = h->next;
  h->next->prev = h->prev;
}

void Cache::lru_insert(head_t *h)
{
  // insert to last position
  h->next = &lru_head;
  h->prev = lru_head.prev;
  h->prev->next = h;
  h->next->prev = h;
}

bool Cache::is_cached(const int index) const
{
  head_t *h = &head[index];
  if(h->len > 0)
    return true;
  return false;
}

int Cache::get_data(const int index, Qfloat **data, int len)
{
  head_t *h = &head[index];
  if(h->len) lru_delete(h);
  int more = len - h->len;

  if(more > 0)
    {
      // free old space
      while(size < more)
	{
	  head_t *old = lru_head.next;
	  lru_delete(old);
	  free(old->data);
	  size += old->len;
	  old->data = 0;
	  old->len = 0;
	}

      // allocate new space
      h->data = (Qfloat *)realloc(h->data,sizeof(Qfloat)*len);
      size -= more;
      swap(h->len,len);
    }

  lru_insert(h);
  *data = h->data;
  return len;
}

void Cache::swap_index(int i, int j)
{
  if(i==j) return;

  if(head[i].len) lru_delete(&head[i]);
  if(head[j].len) lru_delete(&head[j]);
  swap(head[i].data,head[j].data);
  swap(head[i].len,head[j].len);
  if(head[i].len) lru_insert(&head[i]);
  if(head[j].len) lru_insert(&head[j]);

  if(i>j) swap(i,j);
  for(head_t *h = lru_head.next; h!=&lru_head; h=h->next)
    {
      if(h->len > i)
	{
	  if(h->len > j)
	    swap(h->data[i],h->data[j]);
	  else
	    {
	      // give up
	      lru_delete(h);
	      free(h->data);
	      size += h->len;
	      h->data = 0;
	      h->len = 0;
	    }
	}
    }
}

//
// Kernel evaluation
//
// the static method k_function is for doing single kernel evaluation
// the constructor of Kernel prepares to calculate the l*l kernel matrix
// the member function get_Q is for getting one column from the Q Matrix
//
class QMatrix {
public:
  virtual Qfloat *get_Q(int column, int len) const = 0;
  virtual Qfloat *get_QD() const = 0;
  virtual Qfloat *get_Q_subset(int i, int *idxs, int n) const = 0;
  virtual Qfloat get_non_cached(int i, int j) const = 0;
  virtual bool is_cached(int i) const = 0;
  virtual void swap_index(int i, int j) const = 0;
};  

class Kernel: public QMatrix {
public:
  Kernel(int l, Xfloat **x, int **nz_idx, 
	 const int *x_len, const int max_idx, 
	 const svm_parameter& param);
  virtual ~Kernel();

  static double k_function(const Xfloat *x, const int *nz_x, const int lx,
			   Xfloat *y, int *nz_y, int ly, 
			   const svm_parameter& param);
  virtual Qfloat *get_Q(int column, int len) const = 0;
  virtual Qfloat *get_QD() const = 0;
  virtual Qfloat *get_Q_subset(int i, int *idxs, int n) const = 0;
  virtual Qfloat get_non_cached(int i, int j) const = 0;
  virtual bool is_cached(int i) const = 0;
  virtual void swap_index(int i, int j) const	// no so const...
  {
    swap(x[i],x[j]);
    swap(nz_idx[i], nz_idx[j]);
    swap(x_len[i], x_len[j]);
    if(unrolled == i)
      unrolled = j;
    else if(unrolled == j)
      unrolled = i;
    if(x_square) swap(x_square[i],x_square[j]);
  }
protected:

  double (Kernel::*kernel_function)(int i, int j) const;

private:
  Xfloat **x;
  int **nz_idx;
  int *x_len;
  double *x_square;
  // dense unrolled sparse vector
  mutable Xfloat *v; 
  // index of currently unrolled vector
  mutable int unrolled;
  int max_idx;

  // svm_parameter
  const int kernel_type;
  const int degree;
  const double gamma;
  const double coef0;

  static double dot(const Xfloat *x, const int *nz_x, const int lx, 
		    const Xfloat *y, const int *nz_y, const int ly);
  double dot(const int i, const int j) const
  {
    register int k;
    register double sum;
    if(i != unrolled)
      {
	for(k=0; k<x_len[unrolled]; ++k)
	  v[nz_idx[unrolled][k]] = 0;
	unrolled = i;
	for(k=0; k<x_len[i]; ++k)
	  v[nz_idx[i][k]] = x[i][k];
      }
    sum = 0;
    for(k=0; k<x_len[j]; ++k)
      sum += v[nz_idx[j][k]] * x[j][k];
    return sum;
  }
  double kernel_linear(int i, int j) const
  {
    return dot(i,j);
  }
  double kernel_poly(int i, int j) const
  {
    return powi(gamma*dot(i,j)+coef0,degree);
  }
  double kernel_rbf(int i, int j) const
  {
    return exp(-gamma*(x_square[i]+x_square[j]-2*dot(i,j)));
  }
  double kernel_sigmoid(int i, int j) const
  {
    return tanh(gamma*dot(i,j)+coef0);
  }
};

Kernel::Kernel(int l, Xfloat **x_, int **nz_idx_, 
	       const int *x_len_, const int max_idx_,
	       const svm_parameter& param)
  :kernel_type(param.kernel_type), degree(param.degree),
   gamma(param.gamma), coef0(param.coef0)
{
  switch(kernel_type)
    {
    case LINEAR:
      kernel_function = &Kernel::kernel_linear;
      break;
    case POLY:
      kernel_function = &Kernel::kernel_poly;
      break;
    case RBF:
      kernel_function = &Kernel::kernel_rbf;
      break;
    case SIGMOID:
      kernel_function = &Kernel::kernel_sigmoid;
      break;
    }

  clone(x,x_,l);
  clone(nz_idx,nz_idx_,l);
  clone(x_len,x_len_,l);
  max_idx = max_idx_;
  v = new Xfloat[max_idx];
  unrolled = 0;
  for(int k=0; k<x_len[unrolled]; ++k)
    v[nz_idx[unrolled][k]] = x[unrolled][k];

  if(kernel_type == RBF)
    {
      x_square = new double[l];
      for(int i=0;i<l;i++)
	x_square[i] = dot(i,i);
    }
  else
    x_square = 0;
}

Kernel::~Kernel()
{
  delete[] x;
  delete[] nz_idx;
  delete[] x_len;
  delete[] v;
  delete[] x_square;
}

double Kernel::dot(const Xfloat *x, const int *nz_x, const int lx, 
		   const Xfloat *y, const int *nz_y, const int ly)
{
  register double sum = 0;
  register int i = 0; 
  register int j = 0;
  while(i < lx && j < ly)
    {
      if(nz_x[i] == nz_y[j])
	{
	  sum += x[i] * y[j];
	  ++i; ++j;
	}
      else if(nz_x[i] > nz_y[j])
	++j;
      else if(nz_x[i] < nz_y[j])
	++i;
    }
  return sum;
}

double Kernel::k_function(const Xfloat *x, const int *nz_x, const int lx,
			  Xfloat *y, int *nz_y, int ly, 
			  const svm_parameter& param)
{
  switch(param.kernel_type)
    {
    case LINEAR:
      return dot(x, nz_x, lx, y, nz_y, ly);
    case POLY:
      return powi(param.gamma*dot(x, nz_x, lx, y, nz_y, ly)
		  +param.coef0,param.degree);
    case RBF:
      {
	return exp(-param.gamma*(
				 dot(x, nz_x, lx, x, nz_x, lx)+
				 dot(y, nz_y, ly, y, nz_y, ly)-
				 2*dot(x, nz_x, lx, y, nz_y, ly)));
      }
    case SIGMOID:
      return tanh(param.gamma*dot(x, nz_x, lx, y, nz_y, ly)
		  +param.coef0);
    default:
      return 0;	/* Unreachable */
    }
}

// Generalized SMO+SVMlight algorithm
// Solves:
//
//	min 0.5(\alpha^T Q \alpha) + b^T \alpha
//
//		y^T \alpha = \delta
//		y_i = +1 or -1
//		0 <= alpha_i <= Cp for y_i = 1
//		0 <= alpha_i <= Cn for y_i = -1
//
// Given:
//
//	Q, b, y, Cp, Cn, and an initial feasible point \alpha
//	l is the size of vectors and matrices
//	eps is the stopping criterion
//
// solution will be put in \alpha, objective value will be put in obj
//
class Solver {
public:
  Solver() {};
  virtual ~Solver() {};

  struct SolutionInfo {
    double obj;
    double rho;
    double upper_bound_p;
    double upper_bound_n;
    double r;	// for Solver_NU
  };

  void Solve(int l, const QMatrix& Q, const double *b_, const schar *y_,
	     double *alpha_, double Cp, double Cn, double eps,
	     SolutionInfo* si, int shrinking);
protected:
  int active_size;
  schar *y;
  double *G;		// gradient of objective function
  enum { LOWER_BOUND, UPPER_BOUND, FREE };
  char *alpha_status;	// LOWER_BOUND, UPPER_BOUND, FREE
  double *alpha;
  const QMatrix *Q;
  const Qfloat *QD;
  double eps;
  double Cp,Cn;
  double *b;
  int *active_set;
  double *G_bar;		// gradient, if we treat free variables as 0
  int l;
  bool unshrinked;	// XXX

  inline double get_C(int i)
  {
    return (y[i] > 0)? Cp : Cn;
  }
  inline void update_alpha_status(int i)
  {
    if(alpha[i] >= get_C(i))
      alpha_status[i] = UPPER_BOUND;
    else if(alpha[i] <= 0)
      alpha_status[i] = LOWER_BOUND;
    else alpha_status[i] = FREE;
  }
  bool is_upper_bound(int i) { return alpha_status[i] == UPPER_BOUND; }
  bool is_lower_bound(int i) { return alpha_status[i] == LOWER_BOUND; }
  bool is_free(int i) { return alpha_status[i] == FREE; }
  void swap_index(int i, int j);
  void reconstruct_gradient();
  virtual int select_working_set(int &i, int &j);
  virtual int max_violating_pair(int &i, int &j);
  virtual double calculate_rho();
  virtual void do_shrinking();
};

//
// Solver for nu-svm classification and regression
//
// additional constraint: e^T \alpha = constant
//
class Solver_NU : public Solver
{
public:
  Solver_NU() {}
  void Solve(int l, const QMatrix& Q, const double *b, const schar *y,
	     double *alpha, double Cp, double Cn, double eps,
	     SolutionInfo* si, int shrinking);
private:
  SolutionInfo *si;
  int select_working_set(int &i, int &j);
  double calculate_rho();
  void do_shrinking();
};

void Solver_NU::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->si = si;
  Solver::Solve(l,Q,b,y,alpha,Cp,Cn,eps,si,shrinking);
}

void Solver::swap_index(int i, int j)
{
  Q->swap_index(i,j);
  swap(y[i],y[j]);
  swap(G[i],G[j]);
  swap(alpha_status[i],alpha_status[j]);
  swap(alpha[i],alpha[j]);
  swap(b[i],b[j]);
  swap(active_set[i],active_set[j]);
  swap(G_bar[i],G_bar[j]);
}

void Solver::reconstruct_gradient()
{
  // reconstruct inactive elements of G from G_bar and free variables

  if(active_size == l) return;

  int i;
  for(i=active_size;i<l;i++)
    G[i] = G_bar[i] + b[i];
	
  for(i=0;i<active_size;i++)
    if(is_free(i))
      {
	const Qfloat *Q_i = Q->get_Q(i,l);
	double alpha_i = alpha[i];
	for(int j=active_size;j<l;j++)
	  G[j] += alpha_i * Q_i[j];
      }
}

void Solver::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->Q = &Q;
  QD=Q.get_QD();
  clone(b, b_,l);
  clone(y, y_,l);
  clone(alpha,alpha_,l);
  this->Cp = Cp;
  this->Cn = Cn;
  this->eps = eps;
  unshrinked = false;

  // initialize alpha_status
  {
    alpha_status = new char[l];
    for(int i=0;i<l;i++)
      update_alpha_status(i);
  }

  // initialize active set (for shrinking)
  {
    active_set = new int[l];
    for(int i=0;i<l;i++)
      active_set[i] = i;
    active_size = l;
  }

  //  info("initializing gradient..."); info_flush();
  // initialize gradient
  {
    G = new double[l];
    G_bar = new double[l];
    int i;
    for(i=0;i<l;i++)
      {
	G[i] = b[i];
	G_bar[i] = 0;
      }
    for(i=0;i<l;i++)
      if(!is_lower_bound(i))
	{
	  const Qfloat *Q_i = Q.get_Q(i,l);
	  double alpha_i = alpha[i];
	  int j;
	  for(j=0;j<l;j++)
	    G[j] += alpha_i*Q_i[j];
	  if(is_upper_bound(i))
	    for(j=0;j<l;j++)
	      G_bar[j] += get_C(i) * Q_i[j];
	}
  }
  //  info("done.\n"); info_flush();
  // optimization step

  int iter = 0;
  int counter = min(l,1000)+1;

  while(1)
    {
      // show progress and do shrinking

      if(--counter == 0)
	{
	  counter = min(l,1000);
	  if(shrinking) do_shrinking();
	  //  info("."); info_flush();
	}

      int i,j;
      if(select_working_set(i,j)!=0)
	{
	  // reconstruct the whole gradient
	  reconstruct_gradient();