svmregression.java

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

/*
 *  YALE - Yet Another Learning Environment
 *  Copyright (C) 2001-2004
 *      Simon Fischer, Ralf Klinkenberg, Ingo Mierswa, 
 *          Katharina Morik, Oliver Ritthoff
 *      Artificial Intelligence Unit
 *      Computer Science Department
 *      University of Dortmund
 *      44221 Dortmund,  Germany
 *  email: [e-mail]
 *  web:   http://yale.cs.uni-dortmund.de/
 *
 *  This program is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License as 
 *  published by the Free Software Foundation; either version 2 of the
 *  License, or (at your option) any later version. 
 *
 *  This program is distributed in the hope that it will be useful, but
 *  WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 *  General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 *  USA.
 */
package edu.udo.cs.mySVM.SVM;

import edu.udo.cs.mySVM.Kernel.Kernel;
import edu.udo.cs.mySVM.Examples.ExampleSet;

import edu.udo.cs.yale.operator.Operator;

/**
 * Class for regression SVM
 * @author Stefan R?ping, Ingo Mierswa (only Yale additions)
 * @version $Id: SVMregression.java,v 1.5 2004/08/27 11:57:30 ingomierswa Exp $
 */
public class SVMregression extends SVM {  

    public SVMregression() {}

    public SVMregression(Operator paramOperator, Kernel kernel, ExampleSet exampleSet) {
	super(paramOperator, kernel, exampleSet);
    };
  

    /**
     * Calls the optimizer
     */
    protected void optimize()
    {
	// optimizer-specific call

	//	qp.n = working_set_size;
	
	int i;
	int j;

	double[] my_primal = primal;	
	// equality constraint
	qp.b[0]=0;
	for(i=0;i<working_set_size;i++){
	    qp.b[0] += alphas[working_set[i]];
	};
	
	// set initial optimization parameters
	double new_target=0;
	double old_target=0;
	double target_tmp;
	for(i=0;i<working_set_size;i++){
	    target_tmp = my_primal[i]*qp.H[i*working_set_size+i]/2;
	    for(j=0;j<i;j++){
		target_tmp+=my_primal[j]*qp.H[j*working_set_size+i];
	    };
	    target_tmp+=qp.c[i];
	    old_target+=target_tmp*my_primal[i];
	};
	
	double new_constraint_sum=0;
	double my_is_zero = is_zero;
	int sv_count=working_set_size;
	
	// optimize
	boolean KKTerror=true; // KKT not yet satisfied
	boolean convError=false; // KKT can still be satisfied
	
	qp.max_allowed_error = convergence_epsilon;
	
        qp.x = primal;
        qp.solve();
        primal = qp.x;
        lambda_WS = qp.lambda_eq;
	my_primal = primal;

	// loop while some KKT condition is not valid (alpha=0)


	int it=3;

	double lambda_lo;
	while(KKTerror && (it>0)){
	    // iterate optimization 3 times with changed sign on variables, if KKT conditions are not satisfied
	    KKTerror = false;
	    it--;
	    for(i=0;i<working_set_size;i++){
		if(my_primal[i]<is_zero){
		    lambda_lo =  epsilon_neg + epsilon_pos - qp.c[i];
		    for(j=0;j<working_set_size;j++){
			lambda_lo -= my_primal[j]*qp.H[i*working_set_size+j];
		    };
		    if(qp.A[i] > 0){
			lambda_lo -= lambda_WS;
		    }
		    else{
			lambda_lo += lambda_WS;
		    };

		    if(lambda_lo<-convergence_epsilon){
			// change sign of i
			KKTerror=true;
			qp.A[i] = -qp.A[i];
			which_alpha[i] = ! which_alpha[i];
			my_primal[i] = -my_primal[i];
			qp.c[i] = epsilon_neg + epsilon_pos - qp.c[i];
			if(qp.A[i]>0){
			    qp.u[i] = Cneg;
			}
			else{
			    qp.u[i] = Cpos;
			};
			for(j=0;j<working_set_size;j++){
			    qp.H[i*working_set_size+j] = -qp.H[i*working_set_size+j];
			    qp.H[j*working_set_size+i] = -qp.H[j*working_set_size+i];
			};
			if(quadraticLossNeg){
			    if(which_alpha[i]){
				(qp.H)[i*(working_set_size+1)] += 1/Cneg;
				(qp.u)[i] = Double.MAX_VALUE;
			    }
			    else{
				// previous was neg
				(qp.H)[i*(working_set_size+1)] -= 1/Cneg;
			    };
			};
			if(quadraticLossPos){
			    if(! which_alpha[i]){
				(qp.H)[i*(working_set_size+1)] += 1/Cpos;
				(qp.u)[i] = Double.MAX_VALUE;
			    }
			    else{
				//previous was pos
				(qp.H)[i*(working_set_size+1)] -= 1/Cpos;
			    };
			};
		    };
		};
	    };
	    qp.x = my_primal;
	    qp.solve();
	    my_primal = qp.x;
	    lambda_WS = qp.lambda_eq;
	};
	
	KKTerror = true;
	
	while(KKTerror){
	    // clip
	    sv_count=working_set_size;
	    new_constraint_sum=qp.b[0];
	    for(i=0;i<working_set_size;i++){
		// check if at bound
		if(my_primal[i] <= my_is_zero){
		    // at lower bound
		    my_primal[i] = qp.l[i];
		    sv_count--;
		}
		else if(qp.u[i]-my_primal[i] <= my_is_zero){
		    // at upper bound
		    my_primal[i] = qp.u[i];
		    sv_count--;
		};
		new_constraint_sum -= qp.A[i]*my_primal[i];
	    };
	    
	    // enforce equality constraint
	    if(sv_count>0){
		new_constraint_sum /= (double)sv_count;
		logln(5,"adjusting "+sv_count+" alphas by "+new_constraint_sum);
		for(i=0;i<working_set_size;i++){
		    if((my_primal[i] > qp.l[i]) && 
		       (my_primal[i] < qp.u[i])){
			// real sv
			my_primal[i] += qp.A[i]*new_constraint_sum;
		    };
		};
	    }
	    else if(Math.abs(new_constraint_sum)>(double)working_set_size*is_zero){
		// error, can't get feasible point
		logln(5,"WARNING: No SVs, constraint_sum = "+new_constraint_sum);
		old_target = -Double.MIN_VALUE; 
		convError=true;
	    };
	    // test descend
	    new_target=0;
	    for(i=0;i<working_set_size;i++){
		// attention: optimizer changes one triangle of H!
		target_tmp = my_primal[i]*qp.H[i*working_set_size+i]/2.0;
		for(j=0;j<i;j++){
		    target_tmp+=my_primal[j]*qp.H[j*working_set_size+i];
		};
		target_tmp+=qp.c[i];
		new_target+=target_tmp*my_primal[i];
	    };
	    
	    if(new_target < old_target){
		KKTerror = false;
		if(descend < old_target - new_target){
		    target_count=0;
		}
		else{
		    convError=true;
		};
		logln(5,"descend = "+(old_target-new_target));
	    }
	    else if(sv_count > 0){
		// less SVs
		// set my_is_zero to min_i(primal[i]-qp.l[i], qp.u[i]-primal[i])
		my_is_zero = Double.MAX_VALUE;
		for(i=0;i<working_set_size;i++){
		    if((my_primal[i] > qp.l[i]) && (my_primal[i] < qp.u[i])){
			if(my_primal[i] - qp.l[i] < my_is_zero){
			    my_is_zero = my_primal[i]-qp.l[i];
			};
			if(qp.u[i]  - my_primal[i]  < my_is_zero){
			    my_is_zero = qp.u[i] - my_primal[i];
			};
		    };
		};
		if(target_count == 0){
		    my_is_zero *= 2;
		};
		logln(5,"WARNING: no descend ("+(old_target-new_target)
		      +" <= "+descend
		      +"), adjusting is_zero to "+my_is_zero);
		logln(5,"new_target = "+new_target);
	    }
	    else{
		// nothing we can do
		logln(5,"WARNING: no descend ("+(old_target-new_target)
		      +" <= "+descend+"), stopping.");
		KKTerror=false;
		convError=true;
	    };
	};
	
	if(convError){
	    target_count++;
	    if(old_target < new_target){
		for(i=0;i<working_set_size;i++){
		    my_primal[i] = qp.A[i]*alphas[working_set[i]];
		};                              
		logln(5,"WARNING: Convergence error, restoring old primals");
	    };                                          
	};
	
	if(target_count>50){
	    // non-recoverable numerical error
	    convergence_epsilon*=2;
	    feasible_epsilon = convergence_epsilon;
	    logln(1,"WARNING: reducing KKT precision to "+convergence_epsilon);
	    target_count=0;
	};
    };
    


    protected final boolean is_alpha_neg(int i)
    {
	boolean result;
	double alpha = alphas[i];
	if(alpha > 0){
	    result = true;
	}
	else if(alpha == 0){
	    if(sum[i] - ys[i] + lambda_eq>0){
		result = false;
	    }
	    else{
		result = true;
	    };
	}
	else{
	    result = false;
	};
	return result;
    };


    protected final double nabla(int i)
    {
	double alpha =  alphas[i];
	double y = ys[i];
	double result;
	if(alpha > 0){
	    result = ( sum[i] - y + epsilon_neg);
	}
	else if(alpha == 0){
	    if(is_alpha_neg(i)){
		result = ( sum[i] - y + epsilon_neg);
	    }
	    else{
		result = (-sum[i] + y + epsilon_pos);
	    };
	}
	else{
	    result = (-sum[i] + y + epsilon_pos);
	};
	return result;
    };
  
};