svm.java

一个很好的LIBSVM的JAVA源码。对于要研究和改进SVM算法的学者。可以参考。来自数据挖掘工具YALE工具包。

    /**
     * init the optimizer
     */
    protected void init_optimizer()
    {
	which_alpha = new boolean[working_set_size];
	primal = new double[working_set_size];
	sum = new double[examples_total];
	at_bound = new int[examples_total];
	
	// init variables
	if(working_set_size>examples_total) working_set_size = examples_total;
	
	qp = new quadraticProblemSMO(is_zero/100,convergence_epsilon/100,
				     working_set_size*working_set_size);
	qp.set_n(working_set_size);
	
	// reserve workspace for calculate_working_set
	working_set = new int[working_set_size];
	heap_max = new MaxHeap(0);    
	heap_min = new MinHeap(0);    
	
//  	if(the_examples.Dev != null){
//  	    double dev = the_examples.Dev[the_examples.get_dim()];
//  	    if(dev != 0){
//  		epsilon_pos /= dev;
//  		epsilon_neg /= dev;
//  	    };
//  	};
	
	int i;
	for(i=0;i<working_set_size;i++){
	    qp.l[i] = 0; //-is_zero; 
	};
	
	if(quadraticLossPos){
	    Cpos = Double.MAX_VALUE;
	};
	if(quadraticLossNeg){ 
	    Cneg = Double.MAX_VALUE;
	};

	for(i=0;i<examples_total;i++){
	    alphas[i] = 0.0;
	};

	lambda_WS = 0;
	to_shrink=0;
	
	qp.set_n(working_set_size);
    };


    /**
     * exit the optimizer
     */
    protected void exit_optimizer()
    {
	qp = null;
    };


    /**
     * shrink the variables
     */
    protected void shrink()
    {
	// move shrinked examples to back
	if(to_shrink > examples_total/10){
	    int i;
	    int last_pos=examples_total;
	    if(last_pos > working_set_size){
		for(i=0;i<last_pos;i++){
		    if(at_bound[i] >= shrink_const){
			// shrink i
			sum_alpha += alphas[i];
			last_pos--;
			the_examples.swap(i,last_pos);
			the_kernel.swap(i,last_pos);
			sum[i] = sum[last_pos];
			at_bound[i] = at_bound[last_pos];
			if(last_pos <= working_set_size){
			    break;
			};
		    };
		};
		to_shrink=0;
		shrinked = true;
		if(last_pos < examples_total){
		    examples_total = last_pos;
		    the_kernel.set_examples_size(examples_total);
		};
	    };
	    logln(4,"shrinked to "+examples_total+" variables");
	};
    };

    
    /**
     * reset the shrinked variables
     */
    protected void reset_shrinked()
    {
	int old_ex_tot=examples_total;
	target_count=0;
	examples_total = the_examples.count_examples();
	the_kernel.set_examples_size(examples_total);
	// unshrink, recalculate sum for all variables 
	int i,j;
	// reset all sums
	double Kij;
	for(i=old_ex_tot;i<examples_total;i++){
	    sum[i] = 0;
	    at_bound[i] = 0;
	};
	double alpha;
	double[] kernel_row;
	for(i=0;i<examples_total;i++){
	    alpha = alphas[i];
	    if(alpha != 0){
		kernel_row = the_kernel.get_row(i);
		for(j=old_ex_tot;j<examples_total;j++){
		    sum[j]+=alpha*kernel_row[j];
		};
	    };
	};
	sum_alpha=0;
	shrinked = false;
	logln(5,"Resetting shrinked from "+old_ex_tot+" to "+examples_total);
    };
    
    
    /**
     * Project variables to constraints
     */
    protected void project_to_constraint()
    {
      // project alphas to match the constraint
      double alpha_sum = sum_alpha;
      int SVcount=0;
      double alpha;
      int i;
      for(i=0;i<examples_total;i++){
	alpha = alphas[i];
	alpha_sum += alpha;
	if(((alpha>is_zero) && (alpha-Cneg < -is_zero)) ||
	   ((alpha<-is_zero) && (alpha+Cpos > is_zero))){
	  SVcount++;
	};
      };
      if(SVcount > 0){
	// project
	alpha_sum /= (double)SVcount;
	for(i=0;i<examples_total;i++){
	  alpha = alphas[i];
	  if(((alpha>is_zero) && (alpha-Cneg < -is_zero)) ||
	     ((alpha<-is_zero) && (alpha+Cpos > is_zero))){
	    alphas[i] -= alpha_sum;
	  };
	};
      };
    }; 
    
    
    /**
     * Calculates the working set
     * @exception Exception on any error
     */
    protected void calculate_working_set()
    {
      // reset WSS
      if(working_set_size < parameters_working_set_size){
	working_set_size = parameters_working_set_size;
	if(working_set_size>examples_total) working_set_size = examples_total;
      };

      heap_min.init(working_set_size/2);
      heap_max.init(working_set_size/2+working_set_size%2);

      int i=0;
      double sort_value;

      double the_nabla;
      boolean is_feasible;
      int atbound;
      int j;

      while(i<examples_total) {
	is_feasible = feasible(i);
	if(is_feasible){
	  the_nabla = nabla(i);
	  if(is_alpha_neg(i)){
	    sort_value = -the_nabla;  // - : maximum inconsistency approach
	  }
	  else{
	    sort_value = the_nabla;
	  };
	  // add to heaps
	  heap_min.add(sort_value,i);
	  heap_max.add(sort_value,i);
	};
	i++;
      };

      int[] new_ws = heap_min.get_values();
      working_set_size = 0;
      int pos;
      j = heap_min.size();
      for(i=0;i<j;i++){
	working_set[working_set_size] = new_ws[i];
	working_set_size++;
      };
      pos = working_set_size;
      new_ws = heap_max.get_values();
      j = heap_max.size();
      for(i=0;i<j;i++){
	working_set[working_set_size] = new_ws[i];
	working_set_size++;
      };
      if((! heap_min.empty()) &&
	 (! heap_max.empty())){
	if(heap_min.top_value() >= heap_max.top_value()){
	  // there could be the same values in the min- and maxheap,
	  // sort them out (this is very unlikely)
	  j=0;
	  i=0;
	  while(i<pos){
	    // working_set[i] also in max-heap?
	    j=pos;
	    while((j<working_set_size) && 
		  (working_set[j] != working_set[i])){
	      j++;
	    };
	    if(j<working_set_size){
	      // working_set[i] equals working_set[j]
	      // remove j from WS
	      working_set[j] = working_set[working_set_size-1];
	      working_set_size--;
	    }
	    else{
	      i++;
	    };
	  };
	};
      };

      if(target_count > 0){
	// convergence error on last iteration?
	// some more tests on WS
	// unlikely to happen, so speed isn't so important
	// are all variables at the bound?
	int pos_abs;
	boolean bounded_pos=true;
	boolean bounded_neg=true;
	pos=0;
	double alpha;
	while((pos<working_set_size) && (bounded_pos || bounded_neg)){
	  pos_abs = working_set[pos];
	  alpha = alphas[pos_abs];
	  if(is_alpha_neg(pos_abs)){
	    if(alpha-Cneg < -is_zero){
	      bounded_pos = false;
	    };
	    if(alpha > is_zero){
	      bounded_neg = false;
	    };
	  }
	  else{
	    if(alpha+Cneg > is_zero){
	      bounded_neg = false;
	    };
	    if(alpha < -is_zero){
	      bounded_pos = false;
	    };
	  };
	  pos++;
	};
	if(bounded_pos){
	  // all alphas are at upper bound
	  // need alpha that can be moved upward
	  // use alpha with smallest lambda
	  double max_lambda = Double.MAX_VALUE;
	  int max_pos=examples_total;
	  for(pos_abs=0;pos_abs<examples_total;pos_abs++){
	    alpha = alphas[pos_abs];
	    if(is_alpha_neg(pos_abs)){
	      if(alpha-Cneg < -is_zero){
		if(lambda(pos_abs) < max_lambda){
		  max_lambda = lambda(pos_abs);
		  max_pos = pos_abs;
		};
	      };
	    }
	    else{
	      if(alpha < -is_zero){
		if(lambda(pos_abs) < max_lambda){
		  max_lambda = lambda(pos_abs);
		  max_pos = pos_abs;
		};
	      };
	    };
	  };
	  if(max_pos<examples_total){
	    if(working_set_size<parameters_working_set_size){
	      working_set_size++;
	    };
	    working_set[working_set_size-1] = max_pos;
	  };
	}
	else if(bounded_neg){
	  // all alphas are at lower bound
	  // need alpha that can be moved downward
	  // use alpha with smallest lambda
	  double max_lambda = Double.MAX_VALUE;
	  int max_pos=examples_total;
	  for(pos_abs=0;pos_abs<examples_total;pos_abs++){
	    alpha = alphas[pos_abs];
	    if(is_alpha_neg(pos_abs)){
	      if(alpha > is_zero){
		if(lambda(pos_abs) < max_lambda){
		  max_lambda = lambda(pos_abs);
		  max_pos = pos_abs;
		};
	      };
	    }
	    else{
	      if(alpha+Cneg > is_zero){
		if(lambda(pos_abs) < max_lambda){
		  max_lambda = lambda(pos_abs);
		  max_pos = pos_abs;
		};
	      };
	    };
	  };
	  if(max_pos<examples_total){
	    if(working_set_size<parameters_working_set_size){
	      working_set_size++;
	    };
	    working_set[working_set_size-1] = max_pos;
	  };
	};
      };

      if((working_set_size<parameters_working_set_size) &&
	 (working_set_size<examples_total)){
	// use full working set
	pos = (int)(Math.random()*(double)examples_total);
	int ok;
	while((working_set_size<parameters_working_set_size) &&
	      (working_set_size<examples_total)){
	  // add pos into WS if it isn't already
	  ok = 1;
	  for(i=0;i<working_set_size;i++){
	    if(working_set[i] == pos){
	      ok=0;
	      i = working_set_size;
	    };
	  };
	  if(1 == ok){
	    working_set[working_set_size] = pos;
	    working_set_size++;
	  };
	  pos = (pos+1)%examples_total;
	};
      };
      
      int ipos;
      for(ipos=0;ipos<working_set_size;ipos++){
	which_alpha[ipos] = is_alpha_neg(working_set[ipos]);
      };
    };
    

    /**
     * Updates the working set
     */
    protected void update_working_set()
    {
      // setup subproblem
      int i,j;
      int pos_i, pos_j;
      double[] kernel_row;
      double sum_WS;
      boolean[] my_which_alpha = which_alpha;
      for(pos_i=0;pos_i<working_set_size;pos_i++){
	i = working_set[pos_i];
	
	// put row sort_i in hessian 
	kernel_row = the_kernel.get_row(i);

	sum_WS=0;
	//    for(pos_j=0;pos_j<working_set_size;pos_j++){
	for(pos_j=0;pos_j<pos_i;pos_j++){
	  j = working_set[pos_j];
	  // put all elements K(i,j) in hessian, where j in WS
	  if(((! my_which_alpha[pos_j]) && (! my_which_alpha[pos_i])) ||
	     ((my_which_alpha[pos_j]) && (my_which_alpha[pos_i]))){
	    // both i and j positive or negative
	    (qp.H)[pos_i*working_set_size+pos_j] = kernel_row[j];
	    (qp.H)[pos_j*working_set_size+pos_i] = kernel_row[j];
	  }
	  else{
	    // one of i and j positive, one negative
	    (qp.H)[pos_i*working_set_size+pos_j] = -kernel_row[j];
	    (qp.H)[pos_j*working_set_size+pos_i] = -kernel_row[j];
	  };
	};
	for(pos_j=0;pos_j<working_set_size;pos_j++){
	  j = working_set[pos_j];
	  sum_WS+=alphas[j]*kernel_row[j];
	};
	// set main diagonal 
	(qp.H)[pos_i*working_set_size+pos_i] = kernel_row[i];
	
	// linear and box constraints
	if(! my_which_alpha[pos_i]){
	  // alpha
	  (qp.A)[pos_i] = -1;
	  // lin(alpha) = y_i+eps-sum_{i not in WS} alpha_i K_{ij}
	  //            = y_i+eps-sum_i+sum_{i in WS}
	  (qp.c)[pos_i] = ys[i]+epsilon_pos-sum[i]+sum_WS;
	  primal[pos_i] = -alphas[i];
	  (qp.u)[pos_i] = Cpos;
	}
	else{
	  // alpha^*
	  (qp.A)[pos_i] = 1;
	  (qp.c)[pos_i] = -ys[i]+epsilon_neg+sum[i]-sum_WS;
	  primal[pos_i] = alphas[i];
	  (qp.u)[pos_i] = Cneg;