svm_smo.c

机器学习作者tom mitchell的书上代码

/* Copyright (C) 1999 Greg Schohn - [e-mail] */

/* ********************* svm_smo.c **********************
 * John Platt's Sequential Minimal Optimization algorithm 
 * (http://www.research.microsoft.com/~jplatt/smo-book.pdf).
 * This version is based on the modifications proposed by 
 * Keerthi, Shevade, Bhattacharyya & Murthy 
 * (http://guppy.mpe.nus.edu.sg/~mpessk/smo_mod.shtml) 
 */
 
#include <bow/svm.h>
/* the debugging information is out of date since it uses svm_C as its bound
 * instead of values in cvect */
//#define DEBUG

inline int make_set(int bmap_size, int max_list_size, struct set *s) {
  int i;
  s->items = (int *) malloc(sizeof(int)*max_list_size);
  s->loc_map = (int *) malloc(sizeof(int)*bmap_size);

  if (!(s->items && s->loc_map)) {
    return 0;
  }

  for (i=0; i<bmap_size; i++) {
    s->loc_map[i] = -1;
  }

  s->ilength = 0;
  return 1;
}

void free_set(struct set *s) {
  free(s->items);
  free(s->loc_map);
}

inline int set_insert(int a, struct set *s) {
  if (s->loc_map[a] != -1) {
    return 0;
  }
  s->loc_map[a] = s->ilength;
  s->items[s->ilength] = a;
  s->ilength ++;
  return 1;
}

inline int set_delete(int a, struct set *s) {
  if (s->loc_map[a] == -1) {
    return 0;
  }
  if (!(s->ilength == s->loc_map[a] + 1)) {
    /* we need to swap */
    s->items[s->loc_map[a]] = s->items[s->ilength-1];
    s->loc_map[s->items[s->ilength-1]] = s->loc_map[a];
  }
  s->loc_map[a] = -1;
  s->ilength --;
  return 1;
}

inline int set_lookup(int a, struct set *s) {
  return (s->loc_map[a] != -1);
}

#ifdef DEBUG
int c0(double t, int y) { if (t < svm_epsilon_a || t > svm_C - svm_epsilon_a) return 1; else return 0; }
int c1(double t, int y) { if (t < svm_epsilon_a && y == 1) return 0; else return 1; }
int c2(double t, int y) { if (t > svm_C - svm_epsilon_a && y == -1) return 0; else return 1; }
int c3(double t, int y) { if (t > svm_C - svm_epsilon_a && y == 1) return 0; else return 1; }
int c4(double t, int y) { if (t < svm_epsilon_a  && y == -1) return 0; else return 1; }

void check_s(struct set *s, int(*check_weight)(double,int), struct svm_smo_model *m) {
  int i;
  for (i=0; i<s->ilength; i++) {
    if (s->loc_map[s->items[i]] != i) {
      printf("loc_map inconsistent with item list!\n");
      abort();
    }
    if (check_weight(m->weights[s->items[i]],m->yvect[s->items[i]])) {
      printf("wrong set!\n");
      abort();
    }
  }
  return;
}

void check_inv(struct svm_smo_model *m,int ndocs) {
  check_s(&(m->I0), c0, m);
  check_s(&(m->I1), c1, m);
  check_s(&(m->I2), c2, m);
  check_s(&(m->I3), c3, m);
  check_s(&(m->I4), c4, m);

  if (m->I0.ilength + m->I1.ilength + m->I2.ilength + 
      m->I3.ilength + m->I4.ilength != ndocs) {
    abort();
  }
  return;
}
#endif

static int m1=0,m2=0,m3=0,m4=0;

#define PRINT_SMO_PROGRESS(f,ms) (fprintf((f),                          \
       "\r\t\t\t\t\t\tmajor: %d   opt_single: %d/%d   opt_pair: %d/%d     ",\
       (ms)->n_outer, (ms)->n_single_suc, (ms)->n_single_tot,           \
       (ms)->n_pair_suc, (ms)->n_pair_tot))


/* does NOT call kcache_age */
double smo_evaluate_error(struct svm_smo_model *model, int ex) {
  /* do the hyperplane calculation... */
  if (svm_kernel_type == 0) {
    return (evaluate_model_hyperplane(model->W, 0.0, model->docs[ex]) - model->yvect[ex]);
  } else {
    bow_wv **docs;
    int      ndocs;
    double   sum;
    double  *weights;
    int     *yvect;
    int i;

    docs = model->docs;
    ndocs = model->ndocs;
    weights = model->weights;
    yvect = model->yvect;
 
    for (i=0, sum=0.0; i<ndocs; i++) {
      if (weights[i] != 0.0) {
	sum += yvect[i]*weights[i]*svm_kernel_cache(docs[i], docs[ex]);
      }
    }
    return (sum - yvect[ex]);
  }
}

/* in this case we need to compute the obj. function when a2 is at the endpoints */
double calc_eta_hi(int ex1, int ex2, double L, double H, double k11, double k12, 
		   double k22, struct svm_smo_model *ms) {
  bow_wv **docs;
  int ndocs;
  double *weights;
  int *yvect;

  double Lf, Hf,s, gamma;
  double tmp;
  int i;

  docs = ms->docs;
  ndocs = ms->ndocs;
  weights = ms->weights;
  yvect = ms->yvect;

  /* this need not be horribly efficient, since it is only called,
   * for every time eta is non-negative (which is alledgedly rare) */
  {
    double v1,v2;
    double a1, a2;
    int y1, y2;

    for(i=0,v1=v2=0.0; i<ndocs; i++) {
      if (i==ex1 || i==ex2) {
	continue;
      }
      if (weights[i] == 0.0) {
	continue;
      }
      tmp = yvect[i]*weights[i];
      v1 += tmp*svm_kernel_cache(docs[ex1],docs[i]);
      v2 += tmp*svm_kernel_cache(docs[ex2],docs[i]);
    }

    a1 = weights[ex1];
    a2 = weights[ex2];

    y1 = yvect[ex1];
    y2 = yvect[ex2];

#define CALC_W(gamma,s,a2) ((tmp=gamma-s*a2) + a2 - .5*(k11*(tmp*tmp) - k22*a2*a2) \
			    - s*k12*tmp*a2 - y1*tmp*v1 - y2*a2*v2)

    s = y1*y2;
    gamma = a1-s*a2;
    Lf = CALC_W(gamma,s,L);
    Hf = CALC_W(gamma,s,H);
  }

  if (Lf > Hf + svm_epsilon_crit) {
    return L;
  } else if (Lf < Hf - svm_epsilon_crit) {
    return H;
  } else {
    return MAXDOUBLE;
  }
}

/* "tries" to jointly optimize a pair of lagrange weights ...
 * can't always succeed - in those cases, 0 is returned, 1 on success */

/* INVARIANTS: both error[ex1] & error[ex2] must be valid, though 
 * it doesn't matter what set they belong to 

 * all of the weights are feasible & obey the lin. equality
 * constraint when they come in & only this fn plays with the weights */
int opt_pair(int ex1, int ex2, struct svm_smo_model *ms) {
  double   a1, a2;
  double   ao1, ao2;
  float    C1, C2, C_min;
  bow_wv **docs;
  double   diff1, diff2;
  double   e1, e2;
  double   eta;    /* the value of the second deriv. of the obj */
  double   k11, k12, k22;
  int      ndocs;
  double   L, H;
  double  *weights;
  int     *yvect;
  int      y1, y2;

  int i;

  //printf("opt_pair(%d, %d)\n",ex1,ex2);
  //printV("",ms->error,ms->ndocs,"\n");

  if (ex1 == ex2) {
    m1 ++;
    return 0;
  }

  ms->n_pair_tot ++;

  weights = ms->weights;
  yvect = ms->yvect;

  C1 = ms->cvect[ex1];
  C2 = ms->cvect[ex2];
  C_min = MIN(C1, C2);

  y1 = yvect[ex1];
  y2 = yvect[ex2];
  a1 = weights[ex1];
  a2 = weights[ex2];

  if (y1 == y2) {
    H = a1 + a2;
    L = H - C1;
    L = (0 > L) ? 0 : L;
    H = (C2 < H) ? C2 : H;
  } else {
    L = a2 - a1;
    H = L + C1;
    L = (0 > L) ? 0 : L;
    H = (C2 < H) ? C2 : H;
  }

  if (L >= H) {
    m2++;
    return 0;
  }

  docs = ms->docs;
  ndocs = ms->ndocs;

  k12 = svm_kernel_cache(docs[ex1],docs[ex2]);
  k11 = svm_kernel_cache(docs[ex1],docs[ex1]);
  k22 = svm_kernel_cache(docs[ex2],docs[ex2]);

  eta = 2*k12 - k11 - k22;

  //printf("k11,k12,k22,eta:(%f,%f,%f,%f)\n",k11,k12,k22,eta);

  e1 = ms->error[ex1];
  e2 = ms->error[ex2];

  ao2 = a2;

  if (eta < 0) {
    /* a2 still holds weights[j] */
    a2 = a2 - y2*(e1-e2)/eta;
    if (a2 < L) a2 = L;
    else if (a2 > H) a2 = H;

    if (a2 < svm_epsilon_a) {
      a2 = 0;
    } else if (a2 > C2 - svm_epsilon_a) {
      a2 = C2;
    }
  } else {
    a2 = calc_eta_hi(ex1, ex2, L, H, k11, k12, k22, ms);
    if (a2 == MAXDOUBLE)
      return 0;
  }

  if (fabs(a2 - ao2) < svm_epsilon_a) { //*(a2 + ao2 + svm_epsilon_crit)) {
    m4 ++;
    return 0;
  }

  ao1 = weights[ex1];
  a1 = ao1 + y1*y2*(ao2 - a2);

  /* we know that a2 can't be out of the feasible range since we expilicitly
   * tested for this (by clipping) - however, due to prec. problems - a1
   * could be out of range - if it is, we need to make it feasible (to the
   * alpha constraints), since the number is bogus anyway & was caused by
   * precision problems - there's no reason to alter a2 */
  if (a1 < svm_epsilon_a) {
    a1 = 0.0;
  } else if (a1 > C1 - svm_epsilon_a) {
    a1 = C1;
  }

  weights[ex1] = a1;
  weights[ex2] = a2;

  diff1 = y1*(a1 - ao1);
  diff2 = y2*(a2 - ao2);

  /* update the hyperplane */
  if (svm_kernel_type == 0) {
    double *W = ms->W;

    for (i=0; i<docs[ex1]->num_entries; i++) {
      W[docs[ex1]->entry[i].wi] += diff1 * docs[ex1]->entry[i].weight;
    }
    
    for (i=0; i<docs[ex2]->num_entries; i++) {
      W[docs[ex2]->entry[i].wi] += diff2 * docs[ex2]->entry[i].weight;
    }
  }

  /* update the sets (& start to re-evaluate bup & blow) */
  { 
    int j, i, y;
    double a, aold, C, e;
    struct set *s;

    ms->bup = MAXDOUBLE;
    ms->blow = -1*MAXDOUBLE;

    for (j=0, i=ex1, a=a1, aold=ao1, y=y1, C=C1, e=e1; 
	 j<2; 
	 j++, i=ex2, a=a2, aold=ao2, y=y2, C=C2, e=e2) {
      /* the following block also sets bup & blow to preliminary values.
       * this is so that we don't need to repeat these checks when we're 
       * trying to figure out whether or not some */
      if (a < svm_epsilon_a) {
	if (y == 1)  {
	  s = &(ms->I1);
	  if (ms->bup > e) { ms->bup = e; ms->iup = i; }
	} else {
	  s = &(ms->I4);
	  if (ms->blow < e) { ms->blow = e; ms->ilow = i; }
	}
      } else if (a > C - svm_epsilon_a) {
	if (y == 1)  {
	  s = &(ms->I3);
	  if (ms->blow < e) { ms->blow = e; ms->ilow = i; }
	} else {
	  s = &(ms->I2);
	  if (ms->bup > e) { ms->bup = e; ms->iup = i; }
	}
      } else {
	s = &(ms->I0);
	if (ms->blow < e) { ms->blow = e; ms->ilow = i; }
	if (ms->bup > e) { ms->bup = e; ms->iup = i; }	
      }

      if (set_insert(i, s)) { /* if this was actually inserted, 
				 the state of sets has changed, something needs deleted */
	int deleted=0;
	if (aold < svm_epsilon_a) {
	  ms->nsv ++;
	} else if (a < svm_epsilon_a) { /* if this a changed & its zero now, it used to be an SV */
	  ms->nsv --;
	}

	/* there's 12 different possible ways for the sets to change, 
	 * I believe this to be a pretty simple & efficient way to do it... */
	if (y == 1) {
	  if (s != &(ms->I1))
	    deleted = set_delete(i,&(ms->I1));
	  if (!deleted && s != &(ms->I3))
	    deleted = set_delete(i,&(ms->I3));
	} else if (y == -1) {
	  if (s != &(ms->I2)) 
	    deleted = set_delete(i,&(ms->I2));
	  if (!deleted && s != &(ms->I4))
	    deleted = set_delete(i,&(ms->I4));
	}
	if (!deleted) {
	  set_delete(i,&(ms->I0));
	}
      }
    }
  }

  ms->n_pair_suc ++;