svm.java

一个java程序编写的svm支持向量机小程序

	  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-my_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+my_sum[i]-sum_WS;
	  primal[pos_i] = alphas[i];
	  (qp.u)[pos_i] = Cneg;
	};
      };
      if(quadraticLossNeg){
	for(pos_i=0;pos_i<working_set_size;pos_i++){
	  if(my_which_alpha[pos_i]){
	    (qp.H)[pos_i*(working_set_size+1)] += 1/Cneg;
	    (qp.u)[pos_i] = Double.MAX_VALUE;
	  };
	};
      };
      if(quadraticLossPos){
	for(pos_i=0;pos_i<working_set_size;pos_i++){
	  if(! my_which_alpha[pos_i]){
	    (qp.H)[pos_i*(working_set_size+1)] += 1/Cpos;
	    (qp.u)[pos_i] = Double.MAX_VALUE;
	  };
	};
      };
    };


    /**
     * Initialises the working set
     * @exception Exception on any error
     */
    protected void init_working_set()
      throws Exception
    {
      // calculate sum
      int i,j;
      
      project_to_constraint();
      // check bounds!
      if(the_container.initialised_alpha()){
	logln(2,"Initialising variables, this may take some time.");
	for(i=0; i<examples_total;i++){
	  sum[i] = 0;
	  at_bound[i] = 0;
	  for(j=0; j<examples_total;j++){
	    sum[i] += the_container.get_alpha(j)*the_kernel.calculate_K(i,j);
	  };
	};
      }
      else{
	// skip kernel calculation as all alphas = 0
	for(i=0; i<examples_total;i++){
	  sum[i] = 0;
	  at_bound[i] = 0;
	};    
      };
      
      if(the_container.initialised_alpha()){
	calculate_working_set();
      }
      else{
	// first working set is random
	j=0;
	i=0;

	while((i<working_set_size) && (j < examples_total)){
	  working_set[i] = j;
	  if(is_alpha_neg(j)){
	    which_alpha[i] = true;
	  }
	  else{
	    which_alpha[i] = false;
	  };
	  i++;
	  j++;
	};
      };   
      update_working_set();
    };

    
    /**
     * Calls the optimizer
     */
    protected abstract void optimize() throws Exception;
    
    
    /**
     * Stores the optimizer results
     */
    protected void put_optimizer_values()
	throws Exception
    {
      // update nabla, sum, examples.
      // sum[i] += (primal_j^*-primal_j-alpha_j^*+alpha_j)K(i,j)
      // check for |nabla| < is_zero (nabla <-> nabla*)
      //  cout<<"put_optimizer_values()"<<endl;
      int i=0; 
      int j=0;
      int pos_i;
      double the_new_alpha;
      double[] kernel_row;
      double alpha_diff;
      double my_sum[] = sum;
      double[] alphas = the_container.get_alphas();

      pos_i=working_set_size;
      while(pos_i>0){
	pos_i--;
	if(which_alpha[pos_i]){
	  the_new_alpha = primal[pos_i];
	}
	else{
	  the_new_alpha = -primal[pos_i];
	};
	// next three statements: keep this order!
	i = working_set[pos_i];
	alpha_diff = the_new_alpha-alphas[i];
	alphas[i] = the_new_alpha;
	
	if(alpha_diff != 0){
	  // update sum ( => nabla)
	  kernel_row = the_kernel.get_row(i);
	  for(j=examples_total-1;j>=0;j--){
	      my_sum[j] += alpha_diff*kernel_row[j];
	  };
	};
      };
    };
    
    
    /**
     * Checks if the optimization converged
     * @return boolean true optimzation if converged
     */
    protected boolean convergence()
	throws Exception
    {
      double the_lambda_eq = 0;
      int total = 0;
      double alpha_sum=0;
      double alpha=0;
      int i;
      boolean result = true;
      double[] alphas = the_container.get_alphas();
      // actual convergence-test
      total = 0; alpha_sum=0;
      
      for(i=0;i<examples_total;i++){
	alpha = alphas[i];
	alpha_sum += alpha;
	if((alpha>is_zero) && (alpha-Cneg < -is_zero)){
	  // alpha^* = - nabla
	  the_lambda_eq += -nabla(i); //all_ys[i]-epsilon_neg-sum[i];
	  total++;
	}
	else if((alpha<-is_zero) && (alpha+Cpos > is_zero)){
	  // alpha = nabla
	  the_lambda_eq += nabla(i); //all_ys[i]+epsilon_pos-sum[i];
	  total++;
	};
      };

      logln(4,"lambda_eq = "+(the_lambda_eq/total));
      if(total>0){
	lambda_eq = the_lambda_eq / total;
      }
      else{
	// keep WS lambda_eq
	lambda_eq = lambda_WS; //(lambda_eq+4*lambda_WS)/5;
	logln(4,"*** no SVs in convergence(), lambda_eq = "+lambda_eq+".");
      };

      if(target_count>2){
	// estimate lambda from WS
	if(target_count>20){
	  // desperate attempt to get good lambda!
	  lambda_eq = ((40-target_count)*lambda_eq + (target_count-20)*lambda_WS)/20;
	  logln(5,"Re-Re-calculated lambda from WS: "+lambda_eq);
	  if(target_count>40){
	    // really desperate, kick one example out!
	    i = working_set[target_count%working_set_size];
	    if(is_alpha_neg(i)){
	      lambda_eq = -nabla(i);
	    }
	    else{
	      lambda_eq = nabla(i);
	    };
	    logln(5,"set lambda_eq to nabla("+i+"): "+lambda_eq);
	  };
	}
	else{
	  lambda_eq = lambda_WS;
	  logln(5,"Re-calculated lambda_eq from WS: "+lambda_eq);
	};
      };

      // check linear constraint
      if(java.lang.Math.abs(alpha_sum+sum_alpha) > convergence_epsilon){
	// equality constraint violated
	logln(4,"No convergence: equality constraint violated: |"+(alpha_sum+sum_alpha)+"| >> 0");
	project_to_constraint();
	result = false;  
      };
      
      i=0;
      while((i<examples_total) && (result != false)){
	if(lambda(i)>=-convergence_epsilon){
	  i++;
	}
	else{
	    result = false;
	};
      };

      return result;
    };
    

    protected abstract double nabla(int i);


    /**
     * lagrangion multiplier of variable i
     * @param int variable index
     * @return lambda
     */
    protected double lambda(int i)
    {
      double alpha;
      double result;
      if(is_alpha_neg(i)){
	result = - java.lang.Math.abs(nabla(i)+lambda_eq);
      }
      else{
	result = - java.lang.Math.abs(nabla(i)-lambda_eq);
      };
      // default = not at bound
      
      alpha=the_container.get_alpha(i);
      
      if(alpha>is_zero){
	// alpha*
	if(alpha-Cneg >= - is_zero){
	  // upper bound active
	  result = -lambda_eq-nabla(i);
	};
      }
      else if(alpha >= -is_zero){
	// lower bound active
	if(is_alpha_neg(i)){
	  result = nabla(i) + lambda_eq;
	}
	else{
	  result = nabla(i)-lambda_eq;
	};
      }
      else if(alpha+Cpos <= is_zero){
	// upper bound active
	result = lambda_eq - nabla(i);
      };
      return result;
    };


    protected boolean feasible(int i)
	throws Exception
    {
      boolean is_feasible = true;
      double alpha = the_container.get_alpha(i);
      double the_lambda = lambda(i);
      
      if(alpha-Cneg >= - is_zero){
	// alpha* at upper bound
	if(the_lambda >= 0){
	  at_bound[i]++;
	  if(at_bound[i] == shrink_const) to_shrink++;
	}
	else{
	  at_bound[i] = 0;
	};
      }
      else if((alpha<=is_zero) && (alpha >= -is_zero)){
	// lower bound active
	if(the_lambda >= 0){
	  at_bound[i]++;
	  if(at_bound[i] == shrink_const) to_shrink++;
	}
	else{
	  at_bound[i] = 0;
	};
      }
      else if(alpha+Cpos <= is_zero){
	// alpha at upper bound
	if(the_lambda >= 0){
	  at_bound[i]++;
	  if(at_bound[i] == shrink_const) to_shrink++;
	}
	else{
	  at_bound[i] = 0;
	};
      }
      else{
	// not at bound
	at_bound[i] = 0;
      };
      if((the_lambda >= feasible_epsilon) || (at_bound[i] >= shrink_const)){
	is_feasible = false; 
      };
      return is_feasible;
    };


    protected abstract boolean is_alpha_neg(int i);


    /**
     * logs the outputs
     * @param int warning level
     * @param String Message test
     */
    protected void log(int level, String message)
    {
	if(level <= verbosity){
	    System.out.print(message);
	    System.out.flush();
	};
    };

    
    /**
     * log the output plus newline
     * @param int warning level
     * @param String Message test
     */
    protected void logln(int level, String message)
    {
	if(level <= verbosity){
	    System.out.println(message);
	};
    };

  
    /**
     * predict values on the testset with model
     */
    public void predict()
	throws Exception
    {
	int i;
	double prediction;
	double[] example;
	int size = the_container.count_test_examples();
	for(i=0;i<size;i++){
	    example = the_container.get_test_example(i);
	    prediction = predict(example);
	    the_container.set_test_y(i,prediction);
	};
	logln(4,"Prediction generated");
    };


    /**
     * predict a single example
     */
    protected double predict(double[] example)
	throws Exception
    {
	int i;
	double[] sv;
	double the_sum=the_container.get_b();
	double alpha;
	for(i=0;i<examples_total;i++){
	    alpha = the_container.get_alpha(i);
	    if(alpha != 0){
		sv = the_container.get_example(i);
		the_sum += alpha*the_kernel.calculate_K(sv,example);
	    };
	};
	if(the_container.Exp != null){
	    // unscale prediction
	    int dim = the_container.get_dim();
	    double exp = the_container.Exp[dim];
	    double dev = the_container.Dev[dim];
	    if(dev == 0){
		dev = 1;
	    };
	    the_sum = the_sum*dev+exp;
	};
	return the_sum;
    };


    /**
     * check internal variables
     */
    protected void check()
	throws Exception
    {
	double tsum;
	int i,j;
	double s=0;
	for(i=0; i<examples_total;i++){
	  s += the_container.get_alpha(i);
	  tsum = 0;
	  for(j=0; j<examples_total;j++){
	    tsum += the_container.get_alpha(j)*the_kernel.calculate_K(i,j);
	  };
	  if(Math.abs(tsum-sum[i]) > is_zero){
	      logln(1,"ERROR: sum["+i+"] off by "+(tsum-sum[i]));
	      throw(new Exception("ERROR: sum["+i+"] off by "+(tsum-sum[i])));
	  };
	};
	if(Math.abs(s+sum_alpha) > is_zero){
	    logln(1,"ERROR: sum_alpha is off by "+(s+sum_alpha));
	    throw(new Exception("ERROR: sum_alpha is off by "+(s+sum_alpha)));
	};
    };

};