svm_loqo.c

卡内基梅隆大学MaCallum开发的文本分类系统

    /* data has already been clipped/rounded by solve_qp (things within epsilon
     * are rounded, see above) */
    for (i=0; i<n; i++) {
      /* round those alpha's whose values are close to the boundaries */
      if (qd->primal[i] <= svm_epsilon_a) {
	if (a[ws[i]] > svm_epsilon_a) {
	  (*nsv)--;
	}
	a[ws[i]] = 0.0;
      } else {
	if (a[ws[i]] <= svm_epsilon_a) {
	  (*nsv)++;
	}
	if (qd->primal[i] >= qd->ubv[i]-svm_epsilon_a) {
	  a[ws[i]] = qd->ubv[i];
	} else {
	  a[ws[i]] = qd->primal[i];
	}
      }
    }
  }
}

void recompute_gradient(double *s, bow_wv **docs, int *yvect, double *weights, 
			double *old_weights, int *ws, int wss, int total) {
  int i,j;

  fprintf(stderr,"differences:");

  for (i=0; i<total; i++) {
    double tmp = 0.0;
    for (j=0; j<total; j++) {
      tmp += weights[j]*yvect[j]*svm_kernel_cache_lookup(docs[i],docs[j]);
    }
    if (s[i] - tmp > svm_epsilon_a) { 
      fprintf(stderr, "%d diff = %f", i, tmp-s[i]);
    }
    s[i] = tmp;
  }
  fprintf(stderr,"\n");
}

void update_gradient(double *s, bow_wv **docs, int *yvect, double *weights, 
		     double *old_weights, int *ws, int wss, int total) {
  int i,j,k;
  double *wdy;
  int    *wds;  /* those wdy's that are non-zero */

  wdy = (double *) alloca(sizeof(double)*wss);
  wds = (int *) alloca(sizeof(int)*wss);
  
  /* store all of the results early on, so that a potential
   * enormous s can be cycled thru in a cache friendly manner */
  for (k=i=0; i<wss; i++) {
    j = ws[i];
    if (weights[j] != old_weights[j]) {
      wdy[k] = (weights[j] - old_weights[j]) * yvect[j];
      wds[k] = j;
      k++;
    }
  }

  for (i=0; i<total; i++) {
    for (j=0; j<k; j++) {
      s[i] += wdy[j]*svm_kernel_cache(docs[i],docs[wds[j]]);
    }
  }

  kcache_age();
}

double calculate_b(double *s, int *yvect, double *a, float *cvect, int ndocs) {
  int i,j;
  double b, maxgrad, mingrad;

  mingrad = MAXDOUBLE;
  maxgrad = -1*MAXDOUBLE;

  b = 0;
  for (j=i=0; i<ndocs; i++) {
    if (a[i] > svm_epsilon_a) {
      if (a[i] < cvect[i]-svm_epsilon_a) {
	b += s[i] - yvect[i];
	j++;
      } else if (!j) {
	if ((yvect[i] == 1) && (maxgrad<s[i])) {
	  maxgrad = s[i];
	} else if ((yvect[i] == -1) && (mingrad>s[i])) {
	  mingrad = s[i];
	}
      }
    }
  }

  if (j) {
    return (b/j);
  } else {
    assert(maxgrad != MAXDOUBLE);
    return ((maxgrad+mingrad)/2);
  }
}

int check_optimality(double *s, double *a, int *y, float *cvect, double b, int n) {
  double dist, adist, max_dist;

  int i;

  max_dist = 0;

  /* sanity check 
  dist = 0.0;
  for (i=0; i<n; i++) {
    dist += y[i]*a[i];
  }
  if ((dist > svm_epsilon_crit) || (dist < -1*svm_epsilon_crit)) {
    printf("\ndist == %f\n",dist);
    abort();
    }*/

  for(i=0; i<n; i++) {
    dist = (s[i]-b)*y[i];   /* distance from hyperplane*/
    adist = fabs(dist-1.0); /* how far is it from where it should be */

    if(adist > max_dist) {
      if((a[i] < cvect[i]-svm_epsilon_a) && (dist < 1)) {
	//printf("max_dist=%f, (%f-%f)*%d\n", adist, s[i], b, y[i]);
	max_dist = adist;
      }
      if((a[i]>svm_epsilon_a) && (dist > 1)) {
	//printf("max_dist=%f, (%f-%f)*%d\n", adist, s[i], b, y[i]);
	max_dist = adist;
      }
    }
  }

  if (max_dist > svm_epsilon_crit) {  /* termination criterion */
    return (0);
  } else {
    return (1);
  }
}

int build_svm_guts(bow_wv **docs, int *yvect, double *weights, double *b, 
		   double **W, int ndocs, double *s, float *cvect, int *nsv) {
  double       tb;
  int          cwss;        /* current working set size */
  int          n2inc_prec;  /* # of iterations before we try to increase 
			     * the prec. of the solver  */
  double       original_eps_crit; /* global epsilon_crit gets altered, this 
				   * is to set it back */
  double      *original_weights;  /* address of the vector passed in */
  double      *old_weights; /* lagrange multipliers */
  int         *old_ws; /* just for debugging... */
  struct svm_qp qdata;       
  int          qp_cnt;
  struct di   *scratch;     /* scratch area for 2*bsize doubles */
  int         *ws;          /* bsize of these - the current working set */

#ifdef GCSJPRC
  int old_digits=-1;
#endif

  int i,j;


  //recompute_gradient(s, docs, yvect, weights, old_weights, ws, cwss, ndocs);
  
  npr_loqo_failures=0;

  original_eps_crit = svm_epsilon_crit;

  scratch = (struct di *) alloca(sizeof(struct di)*ndocs);
  old_weights = (double *) alloca(sizeof(double)*ndocs);
  ws = (int *) alloca(sizeof(int)*svm_bsize);
  old_ws = (int *) alloca(sizeof(int)*svm_bsize);

  qdata.init_a = (double *) alloca(sizeof(double)*svm_bsize);
  qdata.ce = (double *) alloca(sizeof(double)*svm_bsize);
  qdata.ce0 = (double *) alloca(sizeof(double)); /* only 1 constant in 1 constraint */
  qdata.g = (double *) alloca(sizeof(double)*svm_bsize*svm_bsize); /* hessian */
  qdata.g0 = (double *) alloca(sizeof(double)*svm_bsize);      /* qbn */

  qdata.primal = (double *) alloca(sizeof(double)*svm_bsize*3);
  qdata.dual = (double *) alloca(sizeof(double)*(svm_bsize*2+1));
  qdata.ubv = (double *) alloca(sizeof(double)*svm_bsize/* should be m */);
  qdata.lbv = (double *) alloca(sizeof(double)*svm_bsize);
  
  /* initialize lbv to non-restricting values */
  /* also hit the bottom triangle of the hessian */
  for (i=0; i<svm_bsize; i++) {
    for (j=i;j<svm_bsize;j++) {
      qdata.g[i*svm_bsize+j] = 0.0;
    }
    qdata.lbv[i] = 0.0;
  }

  /* this is what svmlight does, i'm not sure what the bound is used for */
  qdata.bound = svm_C/4.0;
  qdata.digits = INIT_SIGDIGIT;
  qdata.margin = 0.15;
  qdata.init_iter = 500;
  
  for (i=0; i<ndocs; i++) {
    old_weights[i] = weights[i];
  }

  if (svm_weight_style == WEIGHTS_PER_MODEL) {
    kcache_init(ndocs);
  }

  n2inc_prec = LOOSE2LIVE;
  original_weights = NULL;
  qp_cnt = 0;
  cwss = 0;

  while (1) {
    /* the optimality check is first so that when active learning is happening,
     * it becomes a lot quicker - since a update_gradient may not need to be
     * called for a good number of iterations. */
    /* update b */

    tb = calculate_b(s, yvect, weights, cvect, ndocs);

    /* check optimality */
    if (check_optimality(s, weights, yvect, cvect, tb, ndocs)) {
      break;
    }

    qp_cnt++;
    if (svm_verbosity > 1) {
      fprintf(stderr,"%dth iteration of solve_qp\n", qp_cnt);
    } else {
      if (!(qp_cnt % 200)) {
	fprintf(stderr,"\r\t\t\t\t\t\t%dth iteration", qp_cnt);
	fflush(stdout);
      }
    }

    /* put a working set together */
    for (i=0; i<cwss; i++) {
      old_ws[i] = ws[i];
    }
    cwss = get_ws(ws, yvect, weights, s, cvect, ndocs, cwss, svm_bsize, scratch);    

    for (i=j=0; i<cwss; i++) {
      if (old_weights[ws[i]] == weights[ws[i]] && ws[i] == old_ws[i]) {
	j++;
      }
      old_weights[ws[i]] = weights[ws[i]];
      qdata.ubv[i] = cvect[ws[i]];
    }

    /* this detects infinite loops - which shouldn't happen - but... */
#if 0
    if (j == cwss && qdata.digits == old_digits) {
      fprintf(stderr, "Uh-oh - old weights identical to new weights");
      system("echo \"rainbow did a boo-boo - stopping!\" | /usr/sbin/sendmail [e-mail]");
      svm_verbosity = 4;
      fflush(stderr);
      kill(getpid(),SIGSTOP);
    }
#endif

#ifdef GCSJPRC
    old_digits = qdata.digits;
#endif

    /* using the working set, solve the subproblem */
    setup_solve_sub_qp(ws, yvect, weights, docs, &qdata, cwss, nsv);

    /* update s(t) */
    update_gradient(s, docs, yvect, weights, old_weights, ws, cwss, ndocs);

    if (qdata.digits < INIT_SIGDIGIT) {
      if (n2inc_prec) {
	n2inc_prec --;
      } else {
	n2inc_prec = LOOSE2LIVE;
	qdata.digits = INIT_SIGDIGIT;
	/* fprintf(stderr, "LOOSE2LIVE reached... Increasing precision\n"); */
      }
    } else {
      n2inc_prec = LOOSE2LIVE;
    }
  }

  /* make a hyperplane if we can, since they're so fast :) */
  if (svm_kernel_type == 0) {
    int num_words = bow_num_words();
    int i,j,k;

    *W = (double *) malloc(sizeof(double)*num_words);
    for (i=j=0; i<num_words; i++) {
      (*W)[i] = 0.0;
    }

    for (i=j=0; j<*nsv; i++) {
      if (weights[i] != 0.0) {
	for (k=0; k<docs[i]->num_entries; k++) {
	  (*W)[docs[i]->entry[k].wi] += weights[i]*yvect[i]*docs[i]->entry[k].weight;
	}
	j++;
      }
    }
  }

  if (svm_weight_style == WEIGHTS_PER_MODEL) {
    kcache_clear();
  }

  svm_epsilon_crit = original_eps_crit;

  *b = tb;

  return qp_cnt;
}

/* if this gets called, the weight must have been bound */
inline double svm_loqo_tval_to_err(double si, double b, int y) {
  double ytest = si-b;
  if (y == 1) {
    if (ytest < y) {
      return(y - ytest);
    }
  } else {
    if (ytest > y) {
      return(ytest - y);
    }
  }
  return (0.0);
}