svm.cpp

纠错输出编码SVM算法

			y[i] = +1;
		else
			y[i] = -1;

	double sum_pos = nu*l/2;
	double sum_neg = nu*l/2;

	for(i=0;i<l;i++)
		if(y[i] == +1)
		{
			alpha[i] = min(1.0,sum_pos);
			sum_pos -= alpha[i];
		}
		else
		{
			alpha[i] = min(1.0,sum_neg);
			sum_neg -= alpha[i];
		}

	double *zeros = new double[l];

	for(i=0;i<l;i++)
		zeros[i] = 0;

	Solver_NU s;
	s.Solve(l, SVC_Q(*prob,*param,y), zeros, y,
		alpha, 1.0, 1.0, param->eps, si,  param->shrinking);
	double r = si->r;

	info("C = %f\n",1/r);

	for(i=0;i<l;i++)
		alpha[i] *= y[i]/r;

	si->rho /= r;
	si->obj /= (r*r);
	si->upper_bound_p = 1/r;
	si->upper_bound_n = 1/r;

	delete[] y;
	delete[] zeros;
}

static void solve_one_class(
	const svm_problem *prob, const svm_parameter *param,
	double *alpha, Solver::SolutionInfo* si)
{
	int l = prob->l;
	double *zeros = new double[l];
	schar *ones = new schar[l];
	int i;

	int n = (int)(param->nu*prob->l);	// # of alpha's at upper bound

	for(i=0;i<n;i++)
		alpha[i] = 1;
	alpha[n] = param->nu * prob->l - n;
	for(i=n+1;i<l;i++)
		alpha[i] = 0;

	for(i=0;i<l;i++)
	{
		zeros[i] = 0;
		ones[i] = 1;
	}

	Solver s;
	s.Solve(l, ONE_CLASS_Q(*prob,*param), zeros, ones,
		alpha, 1.0, 1.0, param->eps, si, param->shrinking);

	delete[] zeros;
	delete[] ones;
}

static void solve_epsilon_svr(
	const svm_problem *prob, const svm_parameter *param,
	double *alpha, Solver::SolutionInfo* si)
{
	int l = prob->l;
	double *alpha2 = new double[2*l];
	double *linear_term = new double[2*l];
	schar *y = new schar[2*l];
	int i;

	for(i=0;i<l;i++)
	{
		alpha2[i] = 0;
		linear_term[i] = param->p - prob->y[i];
		y[i] = 1;

		alpha2[i+l] = 0;
		linear_term[i+l] = param->p + prob->y[i];
		y[i+l] = -1;
	}

	Solver s;
	s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
		alpha2, param->C, param->C, param->eps, si, param->shrinking);

	double sum_alpha = 0;
	for(i=0;i<l;i++)
	{
		alpha[i] = alpha2[i] - alpha2[i+l];
		sum_alpha += fabs(alpha[i]);
	}
	info("nu = %f\n",sum_alpha/(param->C*l));

	delete[] alpha2;
	delete[] linear_term;
	delete[] y;
}

static void solve_nu_svr(
	const svm_problem *prob, const svm_parameter *param,
	double *alpha, Solver::SolutionInfo* si)
{
	int l = prob->l;
	double C = param->C;
	double *alpha2 = new double[2*l];
	double *linear_term = new double[2*l];
	schar *y = new schar[2*l];
	int i;

	double sum = C * param->nu * l / 2;
	for(i=0;i<l;i++)
	{
		alpha2[i] = alpha2[i+l] = min(sum,C);
		sum -= alpha2[i];

		linear_term[i] = - prob->y[i];
		y[i] = 1;

		linear_term[i+l] = prob->y[i];
		y[i+l] = -1;
	}

	Solver_NU s;
	s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
		alpha2, C, C, param->eps, si, param->shrinking);

	info("epsilon = %f\n",-si->r);

	for(i=0;i<l;i++)
		alpha[i] = alpha2[i] - alpha2[i+l];

	delete[] alpha2;
	delete[] linear_term;
	delete[] y;
}

//
// decision_function
//
struct decision_function
{
	double *alpha;
	double rho;	
	int nSV;
};

decision_function svm_train_one(
	const svm_problem *prob, const svm_parameter *param,
	double Cp, double Cn)
{
	double *alpha = Malloc(double,prob->l);
	Solver::SolutionInfo si;
	switch(param->svm_type)
	{
		case C_SVC:
			solve_c_svc(prob,param,alpha,&si,Cp,Cn);
			break;
		case NU_SVC:
			solve_nu_svc(prob,param,alpha,&si);
			break;
		case ONE_CLASS:
			solve_one_class(prob,param,alpha,&si);
			break;
		case EPSILON_SVR:
			solve_epsilon_svr(prob,param,alpha,&si);
			break;
		case NU_SVR:
			solve_nu_svr(prob,param,alpha,&si);
			break;
	}

	info("obj = %f, rho = %f\n",si.obj,si.rho);

	// output SVs

	int nSV = 0;
	int nBSV = 0;
	for(int i=0;i<prob->l;i++)
	{
		if(fabs(alpha[i]) > 0)
		{
			++nSV;
			if(prob->y[i] > 0)
			{
				if(fabs(alpha[i]) >= si.upper_bound_p)
					++nBSV;
			}
			else
			{
				if(fabs(alpha[i]) >= si.upper_bound_n)
					++nBSV;
			}
		}
	}

	info("nSV = %d, nBSV = %d\n",nSV,nBSV);

	decision_function f;
	f.alpha = alpha;
	f.rho = si.rho;
	f.nSV = nSV;
	return f;
}

//
// svm_model
//
struct svm_model
{
	svm_parameter param;	// parameter
	int nr_class;		// number of classes, = 2 in regression/one class svm
	int nr_binary;		// number of binary SVMs
	int **I;                // I matrix for different multiclass types
	int l;			// total #SV
	svm_node **SV;		// SVs (SV[l])
	double **sv_coef;	// coefficients for SVs in decision functions
	int **sv_ind;	        // index of SVs
	double *rho;		// constants in decision functions (rho[n*(n-1)/2])
	double *probA;          // pariwise probability information
	double *probB;

	
	// for classification only

	int *label;		// label of each class (label[n])
	int *nSV;		// number of SVs for each class (nSV[n])
				// nSV[0] + nSV[1] + ... + nSV[n-1] = l
	int *nSV_binary;		// number of SVs for each binary SVM
	// XXX
	int free_sv;		// 1 if svm_model is created by svm_load_model
				// 0 if svm_model is created by svm_train
};

// Platt's binary SVM Probablistic Output: an improvement from Lin et al.
void sigmoid_train(
	int l, const double *dec_values, const double *labels, 
	double& A, double& B)
{
	double prior1=0, prior0 = 0;
	int i;

	for (i=0;i<l;i++)
		if (labels[i] > 0) prior1+=1;
		else prior0+=1;
	
	int max_iter=100; 	// Maximal number of iterations
	double min_step=1e-10;	// Minimal step taken in line search
	double sigma=1e-3;	// For numerically strict PD of Hessian
	double eps=1e-5;
	double hiTarget=(prior1+1.0)/(prior1+2.0);
	double loTarget=1/(prior0+2.0);
	double *t=Malloc(double,l);
	double fApB,p,q,h11,h22,h21,g1,g2,det,dA,dB,gd,stepsize;
	double newA,newB,newf,d1,d2;
	int iter; 
	
	// Initial Point and Initial Fun Value
	A=0.0; B=log((prior0+1.0)/(prior1+1.0));
	double fval = 0.0;

	for (i=0;i<l;i++)
	{
		if (labels[i]>0) t[i]=hiTarget;
		else t[i]=loTarget;
		fApB = dec_values[i]*A+B;
		if (fApB>=0)
			fval += t[i]*fApB + log(1+exp(-fApB));
		else
			fval += (t[i] - 1)*fApB +log(1+exp(fApB));
	}
	for (iter=0;iter<max_iter;iter++)
	{
		// Update Gradient and Hessian (use H' = H + sigma I)
		h11=sigma; // numerically ensures strict PD
		h22=sigma;
		h21=0.0;g1=0.0;g2=0.0;
		for (i=0;i<l;i++)
		{
			fApB = dec_values[i]*A+B;
			if (fApB >= 0)
			{
				p=exp(-fApB)/(1.0+exp(-fApB));
				q=1.0/(1.0+exp(-fApB));
			}
			else
			{
				p=1.0/(1.0+exp(fApB));
				q=exp(fApB)/(1.0+exp(fApB));
			}
			d2=p*q;
			h11+=dec_values[i]*dec_values[i]*d2;
			h22+=d2;
			h21+=dec_values[i]*d2;
			d1=t[i]-p;
			g1+=dec_values[i]*d1;
			g2+=d1;
		}

		// Stopping Criteria
		if (fabs(g1)<eps && fabs(g2)<eps)
			break;

		// Finding Newton direction: -inv(H') * g
		det=h11*h22-h21*h21;
		dA=-(h22*g1 - h21 * g2) / det;
		dB=-(-h21*g1+ h11 * g2) / det;
		gd=g1*dA+g2*dB;


		stepsize = 1; 		// Line Search
		while (stepsize >= min_step)
		{
			newA = A + stepsize * dA;
			newB = B + stepsize * dB;

			// New function value
			newf = 0.0;
			for (i=0;i<l;i++)
			{
				fApB = dec_values[i]*newA+newB;
				if (fApB >= 0)
					newf += t[i]*fApB + log(1+exp(-fApB));
				else
					newf += (t[i] - 1)*fApB +log(1+exp(fApB));
			}
			// Check sufficient decrease
			if (newf<fval+0.0001*stepsize*gd)
			{
				A=newA;B=newB;fval=newf;
				break;
			}
			else
				stepsize = stepsize / 2.0;
		}

		if (stepsize < min_step)
		{
			info("Line search fails in two-class probability estimates\n");
			break;
		}
	}

	if (iter>=max_iter)
		info("Reaching maximal iterations in two-class probability estimates\n");
	free(t);
}

double sigmoid_predict(double decision_value, double A, double B)
{
	double fApB = decision_value*A+B;
	if (fApB >= 0)
		return exp(-fApB)/(1.0+exp(-fApB));
	else
		return 1.0/(1+exp(fApB)) ;
}

// Method 2 from the multiclass_prob paper by Wu, Lin, and Weng
void multiclass_probability(int k, double **r, double *p)
{
	int t;
	int iter = 0, max_iter=100;
	double **Q=Malloc(double *,k);
	double *Qp=Malloc(double,k);
	double pQp, eps=0.001;
	
	for (t=0;t<k;t++)
	{
		p[t]=1.0/k;  // Valid if k = 1
		Q[t]=Malloc(double,k);
		Q[t][t]=0;
		for (int j=0;j<t;j++)
		{
			Q[t][t]+=r[j][t]*r[j][t];
			Q[t][j]=Q[j][t];
		}
		for (int j=t+1;j<k;j++)
		{
			Q[t][t]+=r[j][t]*r[j][t];
			Q[t][j]=-r[j][t]*r[t][j];
		}
	}
	for (iter=0;iter<max_iter;iter++)
	{
		// stopping condition, recalculate QP,pQP for numerical accuracy
		pQp=0;
		for (t=0;t<k;t++)
		{
			Qp[t]=0;
			for (int j=0;j<k;j++)
				Qp[t]+=Q[t][j]*p[j];
			pQp+=p[t]*Qp[t];
		}
		double max_error=0;
		for (t=0;t<k;t++)
		{
			double error=fabs(Qp[t]-pQp);
			if (error>max_error)
				max_error=error;
		}
		if (max_error<eps) break;
		
		for (t=0;t<k;t++)
		{
			double diff=(-Qp[t]+pQp)/Q[t][t];
			p[t]+=diff;
			pQp=(pQp+diff*(diff*Q[t][t]+2*Qp[t]))/(1+diff)/(1+diff);
			for (int j=0;j<k;j++)
			{
				Qp[j]=(Qp[j]+diff*Q[t][j])/(1+diff);
				p[j]/=(1+diff);
			}
		}
	}
	if (iter>=max_iter)
		info("Exceeds max_iter in multiclass_prob\n");
	for(t=0;t<k;t++) free(Q[t]);
	free(Q);
	free(Qp);
}

void q(double *qp, double *qn, double *b, int **I, int m, int k){
	for(int t=0; t<m; t++)
		qp[t] = qn[t] = 0.0;
		
	for(int t=0; t<m; t++)
		for(int j=0; j<k; j++)
			if(I[t][j]>0) qp[t] += exp(b[j]);
			else if(I[t][j]<0) qn[t] += exp(b[j]);
}


double Obj(double *rp, double *qp, double *qn, double *b, double mu, int m, int k){
	int i; 
	double obj = 0.0, sumexpb = 0;
	
	for(i=0; i<m; i++)
		obj -= (rp[i]*log(qp[i]) + (1-rp[i])*log(qn[i]) - log(qp[i]+qn[i]));
	
	for(i=0; i<k; i++)
		sumexpb += exp(b[i]);

	for(i=0; i<k; i++)
		obj -= mu*(b[i]-log(sumexpb));
	
	return obj;
}

// Gradient decent
void grad(int multiclass_type, int m, int k, double *rp, double *p, int **I){
	int iter=0,max_iter=10000;
	double mu=0.0, eps=0.001, sigma = 0.01, delta = 1,
	       *qp=Malloc(double, m),
	       *qn=Malloc(double, m);
	double *b = Malloc(double, k);
	double *bt = Malloc(double, k);
	double *grad = Malloc(double, k);
	double *d = Malloc(double, k);
	double sumexpb = 0.0, gradnorm, obj, max_error, newobj;
	
	if(multiclass_type == DENSE || multiclass_type == SPARSE)
		mu = 0.001;

	// initialize beta
	for(int i=0; i<k; i++)
		b[i] = log(1.0/k);
	
	q(qp,qn,b,I,m,k);
        obj = Obj(rp, qp, qn, b, mu, m, k);	
	for(iter=0; iter<max_iter; iter++)
	{
		// Calculate gradient
		sumexpb = 0;
		for(int p=0; p<k; p++)
			sumexpb += exp(b[p]);
		
		for(int i=0; i<k; i++)
		{
			double delta_d = 0.0, delta_n = 0.0; 	
			for(int j=0; j<m; j++)
			{
				if(I[j][i]>0)
				{
					delta_n += rp[j]/qp[j];
					delta_d += 1.0/(qp[j]+qn[j]);
				}
				else if(I[j][i]<0)
				{
					delta_n += (1.0-rp[j])/qn[j];
					delta_d += 1.0/(qp[j]+qn[j]);
				}
			}
			
			grad[i] = -(delta_n*exp(b[i])+mu) + (delta_d+mu*k/sumexpb)*exp(b[i]); 
			d[i] = (delta_n+mu/exp(b[i]))/(delta_d+mu*k/sumexpb);
			
		}
	
 		max_error = 0;
		for(int i=0; i<k; i++)
			if(fabs(d[i]-1) > max_error)
				max_error = fabs(d[i]-1);

		if(max_error < eps)
			break;
		
		gradnorm = 0;
		for(int i=0; i<k; i++)
			gradnorm += grad[i]*grad[i]; 	
		
		// Search for step size
		delta = 1;
		while(1)
		{
			// Calculate bt
			//printf("%f\n",delta);
			for(int i=0; i<k; i++)
				bt[i] = b[i] - delta * grad[i];
			q(qp,qn,bt,I,m,k);
			newobj = Obj(rp,qp,qn,bt,mu,m,k);
			if(newobj-obj+sigma*delta*gradnorm <= 0)
				break;
			delta = delta*0.1;
			
		}
	
		memmove(b,bt,k*sizeof(double));	
		obj = newobj;
	
	}
		
	sumexpb = 0;
	for(int i=0; i<k; i++)
		sumexpb += exp(b[i]);
	
	for(int i=0; i<k; i++)
		p[i] = exp(b[i])/sumexpb;
	
	if (iter>=max_iter)
		info("Exceeds max_iter in GBT_multiclass_prob: %lf\navg iter: %d\n", max_error,iter*k);
	else
		info("avg iter: %d\n", iter*k);
	
	free(qp);
	free(qn);
	free(b);
	free(bt);
	free(grad);
	free(d);
}


// General Bradley-Terry Model, update all k variables in one iteration

void gbtall(int multiclass_type, int m, int k, double *rp, double *p, int **I)
{
	int iter=0,max_iter=10000;
	double mu=0.0, eps=0.001,
	       *qp=Malloc(double, m),
	       *qn=Malloc(double, m);
	
	if(multiclass_type == DENSE || multiclass_type == SPARSE)
		mu = 0.001;

	// initialize p
	for(int i=0; i<k; i++)
		p[i] = 1.0/k;

	// Algorithm 2
	double *delta=Malloc(double, k);
	for(iter=0; iter<max_iter; iter++)
	{
		for(int t=0; t<m; t++)
			qp[t] = qn[t] = 0.0;
		
		for(int t=0; t<m; t++)
			for(int j=0; j<k; j++)
				if(I[t][j]>0) qp[t] += p[j];
				else if(I[t][j]<0) qn[t] += p[j];
		
		for(int i=0; i<k; i++)
		{
			double delta_d = 0.0, delta_n = 0.0; 	
			for(int j=0; j<m; j++)
			{
				if(I[j][i]>0)
				{
					delta_n += rp[j]/qp[j];
					delta_d += 1.0/(qp[j]+qn[j]);
				}
				else if(I[j][i]<0)
				{
					delta_n += (1.0-rp[j])/qn[j];
					delta_d += 1.0/(qp[j]+qn[j]);
				}
			}
			delta[i] = (delta_n+mu/p[i])/(delta_d+mu*k);
			p[i] = delta[i]*p[i];
		}
		
		
		double sump = 0.0;
		for(int t=0; t<k; t++)
			sump += p[t];
		for(int t=0; t<k; t++)
			p[t] = p[t]/sump;
		
		double max_error = 0.0;
		for(int i=0; i<k; i++)
			if(fabs(delta[i]-1.0)>max_error)
				max_error = fabs(delta[i]-1.0);
		
		if(max_error < eps) break;
	}
	
	if (iter>=max_iter)
		info("Exceeds max_iter in GBT_multiclass_prob\navg iter: %d\n",iter*k);
	else
		info("avg iter: %d\n", iter*k);
	
	free(qp);
	free(qn);
	free(delta);
}

// Generalized Bradley-Terry Model

void gbtone(int multiclass_type, int m, int k, double *rp, double *p, int **I)
{
	int iter=0,max_iter=10000;
	double mu=0.0, eps=0.001,
	       *qp=Malloc(double, m),
	       *qn=Malloc(double, m);
	
	if(multiclass_type == DENSE || multiclass_type == SPARSE)
		mu = 0.001;

	// initialize p