svm_hideo.h

支持向量机程序(2)Text Windows Svm

							
							
							return((int)result);
}


int solve_dual(
			   /* Solves the dual using the method of Hildreth and D'Espo. */
			   /* Can only handle problems with zero or exactly one */
			   /* equality constraints. */
			   
			   long   n,            /* number of variables */
			   long   m,            /* number of linear equality constraints */
			   double precision,    /* solve at least to this dual precision */
			   double epsilon_crit, /* stop, if KT-Conditions approx fulfilled */
			   long   maxiter,      /* stop after that many iterations */
			   double *g,
			   double *g0,          /* linear part of objective */
			   double *ce,double *ce0,     /* linear equality constraints */
			   double *low,double *up,     /* box constraints */
			   double *primal,      /* variables (with initial values) */
			   double *d,double *d0,double *ig,double *dual,double *dual_old,double *temp,       /* buffer  */
			   long goal)
{
	long i,j,k,iter;
	double sum,w,maxviol,viol,temp1,temp2,isnantest;
	double model_b,dist;
	long retrain,maxfaktor,primal_optimal=0,at_bound,scalemaxiter;
	double epsilon_a=1E-15,epsilon_hideo;
	double eq; 
	
	if((m<0) || (m>1)) 
		printe("SOLVE DUAL: inappropriate number of eq-constrains!");
	
	
	for(i=0;i<2*(n+m);i++) {
		dual[i]=0;
		dual_old[i]=0;
	}
	for(i=0;i<n;i++) {   
		for(j=0;j<n;j++) {   /* dual hessian for box constraints */
			d[i*2*(n+m)+j]=ig[i*n+j];
			d[(i+n)*2*(n+m)+j]=-ig[i*n+j];
			d[i*2*(n+m)+j+n]=-ig[i*n+j];
			d[(i+n)*2*(n+m)+j+n]=ig[i*n+j];
		}
		if(m>0) {
			sum=0;              /* dual hessian for eq constraints */
			for(j=0;j<n;j++) {
				sum+=(ce[j]*ig[i*n+j]);
			}
			d[i*2*(n+m)+2*n]=sum;
			d[i*2*(n+m)+2*n+1]=-sum;
			d[(n+i)*2*(n+m)+2*n]=-sum;
			d[(n+i)*2*(n+m)+2*n+1]=sum;
			d[(n+n)*2*(n+m)+i]=sum;
			d[(n+n+1)*2*(n+m)+i]=-sum;
			d[(n+n)*2*(n+m)+(n+i)]=-sum;
			d[(n+n+1)*2*(n+m)+(n+i)]=sum;
			
			sum=0;
			for(j=0;j<n;j++) {
				for(k=0;k<n;k++) {
					sum+=(ce[k]*ce[j]*ig[j*n+k]);
				}
			}
			d[(n+n)*2*(n+m)+2*n]=sum;
			d[(n+n)*2*(n+m)+2*n+1]=-sum;
			d[(n+n+1)*2*(n+m)+2*n]=-sum;
			d[(n+n+1)*2*(n+m)+2*n+1]=sum;
		} 
	}
	
	for(i=0;i<n;i++) {   /* dual linear component for the box constraints */
		w=0;
		for(j=0;j<n;j++) {
			w+=(ig[i*n+j]*g0[j]); 
		}
		d0[i]=up[i]+w;
		d0[i+n]=-low[i]-w;
	}
	
	if(m>0) {  
		sum=0;             /* dual linear component for eq constraints */
		for(j=0;j<n;j++) {
			for(k=0;k<n;k++) {
				sum+=(ce[k]*ig[k*n+j]*g0[j]); 
			}
		}
		d0[2*n]=ce0[0]+sum;
		d0[2*n+1]=-ce0[0]-sum;
	}
	
	maxviol=999999;
	iter=0;
	retrain=1;
	maxfaktor=1;
	scalemaxiter=maxiter/5;
	while((retrain) && (maxviol > 0) && (iter < (scalemaxiter*maxfaktor))) {
		iter++;
		
		while((maxviol > precision) && (iter < (scalemaxiter*maxfaktor))) {
			iter++;
			maxviol=0;
			for(i=0;i<2*(n+m);i++) {
				sum=d0[i];
				for(j=0;j<2*(n+m);j++) {
					sum+=d[i*2*(n+m)+j]*dual_old[j];
				}
				sum-=d[i*2*(n+m)+i]*dual_old[i];
				dual[i]=-sum/d[i*2*(n+m)+i];
				if(dual[i]<0) dual[i]=0;
				
				viol=fabs(dual[i]-dual_old[i]);
				if(viol>maxviol) 
					maxviol=viol;
				dual_old[i]=dual[i];
			}
			/*
			sprintf(temstr,"%d) maxviol=%20f precision=%f\n",iter,maxviol,precision); 
			*/
		}
		
		if(m>0) {
			for(i=0;i<n;i++) {
				temp[i]=dual[i]-dual[i+n]+ce[i]*(dual[n+n]-dual[n+n+1])+g0[i];
			}
		} 
		else {
			for(i=0;i<n;i++) {
				temp[i]=dual[i]-dual[i+n]+g0[i];
			}
		}
		for(i=0;i<n;i++) {
			primal[i]=0;             /* calc value of primal variables */
			for(j=0;j<n;j++) {
				primal[i]+=ig[i*n+j]*temp[j];
			}
			primal[i]*=-1.0;
			if(primal[i]<=(low[i])) {  /* clip conservatively */
				primal[i]=low[i];
			}
			else if(primal[i]>=(up[i])) {
				primal[i]=up[i];
			}
		}
		
		if(m>0) 
			model_b=dual[n+n+1]-dual[n+n];
		else
			model_b=0;
		
		epsilon_hideo=EPSILON_HIDEO;
		for(i=0;i<n;i++) {           /* check precision of alphas */
			isnantest+=primal[j];
			dist=-model_b*ce[i]; 
			dist+=(g0[i]+1.0);
			for(j=0;j<i;j++) {
				dist+=(primal[j]*g[j*n+i]);
			}
			for(j=i;j<n;j++) {
				dist+=(primal[j]*g[i*n+j]);
			}
			if((primal[i]<(up[i]-epsilon_hideo)) && (dist < (1.0-epsilon_crit))) {
				epsilon_hideo=(up[i]-primal[i])*2.0;
			}
			else if((primal[i]>(low[i]+epsilon_hideo)) &&(dist>(1.0+epsilon_crit))) {
				epsilon_hideo=(primal[i]-low[i])*2.0;
			}
		}
		/*sprintf(temstr,"\nEPSILON_HIDEO=%.30f\n",epsilon_hideo); */
		
		
		for(i=0;i<n;i++) {           /* clip alphas to bounds */
			if(primal[i]<=(low[i]+epsilon_hideo)) {
				primal[i]=low[i];
			}
			else if(primal[i]>=(up[i]-epsilon_hideo)) {
				primal[i]=up[i];
			}
		}
		
		retrain=0;
		primal_optimal=1;
		at_bound=0;
		for(i=0;(i<n);i++) {  /* check primal KT-Conditions */
			dist=-model_b*ce[i]; 
			dist+=(g0[i]+1.0);
			for(j=0;j<i;j++) {
				dist+=(primal[j]*g[j*n+i]);
			}
			for(j=i;j<n;j++) {
				dist+=(primal[j]*g[i*n+j]);
			}
			if((primal[i]<(up[i]-epsilon_a)) && (dist < (1.0-epsilon_crit))) {
				retrain=1;
				primal_optimal=0;
			}
			else if((primal[i]>(low[i]+epsilon_a)) && (dist > (1.0+epsilon_crit))) {
				retrain=1;
				primal_optimal=0;
			}
			if((primal[i]<=(low[i]+epsilon_a)) || (primal[i]>=(up[i]-epsilon_a))) {
				at_bound++;
			}
			/*    sprintf(temstr,"HIDEOtemp: a[%ld]=%.30f, dist=%.6f, b=%f, at_bound=%ld\n",i,primal[i],dist,model_b,at_bound);  */
		}
		if(m>0) {
			eq=-ce0[0];               /* check precision of eq-constraint */
			for(i=0;i<n;i++) { 
				eq+=(ce[i]*primal[i]);
			}
			if((EPSILON_EQ < fabs(eq)) 
			/*
			&& !((goal==PRIMAL_OPTIMAL) 
			&& (at_bound==n)) */
			) {
				retrain=1;
				primal_optimal=0;
			}
			/* sprintf(temstr,"\n eq=%.30f ce0=%f at-bound=%ld\n",eq,ce0[0],at_bound);  */
		}
		
		if(retrain) {
			precision/=10;
			if(((goal == PRIMAL_OPTIMAL) && (maxfaktor < 50000))
				|| (maxfaktor < 5)) {
				maxfaktor++;
			}
		}
  }
  
  if(!primal_optimal) {
	  for(i=0;i<n;i++) {
		  primal[i]=0;             /* calc value of primal variables */
		  for(j=0;j<n;j++) {
			  primal[i]+=ig[i*n+j]*temp[j];
		  }
		  primal[i]*=-1.0;
		  if(primal[i]<=(low[i]+epsilon_a)) {  /* clip conservatively */
			  primal[i]=low[i];
		  }
		  else if(primal[i]>=(up[i]-epsilon_a)) {
			  primal[i]=up[i];
		  }
	  }
  }
  
  isnantest=0;
  for(i=0;i<n;i++) {           /* check for isnan */
	  isnantest+=primal[i];
  }
  
  if(m>0) {
	  temp1=dual[n+n+1];   /* copy the dual variables for the eq */
	  temp2=dual[n+n];     /* constraints to a handier location */
	  for(i=n+n+1;i>=2;i--) {
		  dual[i]=dual[i-2];
	  }
	  dual[0]=temp2;
	  dual[1]=temp1;
	  isnantest+=temp1+temp2;
  }
  
  if(isnan(isnantest)) {
	  return((int)NAN_SOLUTION);
  }
  else if(primal_optimal) {
	  return((int)PRIMAL_OPTIMAL);
  }
  else if(maxviol == 0.0) {
	  return((int)DUAL_OPTIMAL);
  }
  else {
	  return((int)MAXITER_EXCEEDED);
  }
}


void linvert_matrix(
					double *matrix,
					long depth,
					double *inverse,double lindep_sensitivity,
					long *lin_dependent)  /* indicates the active parts of matrix on 
					input and output*/
{
	long i,j,k;
	double factor;
	
	for(i=0;i<depth;i++) {
		/*    lin_dependent[i]=0; */
		for(j=0;j<depth;j++) {
			inverse[i*depth+j]=0.0;
		}
		inverse[i*depth+i]=1.0;
	}
	for(i=0;i<depth;i++) {
		if(lin_dependent[i] || (fabs(matrix[i*depth+i])<lindep_sensitivity)) {
			lin_dependent[i]=1;
		}
		else {
			for(j=i+1;j<depth;j++) {
				factor=matrix[j*depth+i]/matrix[i*depth+i];
				for(k=i;k<depth;k++) {
					matrix[j*depth+k]-=(factor*matrix[i*depth+k]);
				}
				for(k=0;k<depth;k++) {
					inverse[j*depth+k]-=(factor*inverse[i*depth+k]);
				}
			}
		}
	}
	for(i=depth-1;i>=0;i--) {
		if(!lin_dependent[i]) {
			factor=1/matrix[i*depth+i];
			for(k=0;k<depth;k++) {
				inverse[i*depth+k]*=factor;
			}
			matrix[i*depth+i]=1;
			for(j=i-1;j>=0;j--) {
				factor=matrix[j*depth+i];
				matrix[j*depth+i]=0;
				for(k=0;k<depth;k++) {
					inverse[j*depth+k]-=(factor*inverse[i*depth+k]);
				}
			}
		}
	}
}

void lprint_matrix(
				   double *matrix,
				   long depth)
{
	long i,j;
	for(i=0;i<depth;i++) {
		for(j=0;j<depth;j++) {
			sprintf(temstr,"%5.2f ",(double)(matrix[i*depth+j]));
		}
		sprintf(temstr,"\n");
	}
	sprintf(temstr,"\n");
}

void ladd_matrix(
				 double *matrix,
				 long depth,
				 double scalar)
{
	long i,j;
	for(i=0;i<depth;i++) {
		for(j=0;j<depth;j++) {
			matrix[i*depth+j]+=scalar;
		}
	}
}

void lcopy_matrix(
				  double *matrix,
				  long depth,
				  double *matrix2)
{
	long i;
	
	for(i=0;i<(depth)*(depth);i++) {
		matrix2[i]=matrix[i];
	}
}

void lswitch_rows_matrix(
						 double *matrix,
						 long depth,long r1,long r2)
{
	long i;
	double temp;
	
	for(i=0;i<depth;i++) {
		temp=matrix[r1*depth+i];
		matrix[r1*depth+i]=matrix[r2*depth+i];
		matrix[r2*depth+i]=temp;
	}
}

void lswitchrk_matrix(
					  double *matrix,
					  long depth,long rk1,long rk2)
{
	long i;
	double temp;
	
	for(i=0;i<depth;i++) {
		temp=matrix[rk1*depth+i];
		matrix[rk1*depth+i]=matrix[rk2*depth+i];
		matrix[rk2*depth+i]=temp;
	}
	for(i=0;i<depth;i++) {
		temp=matrix[i*depth+rk1];
		matrix[i*depth+rk1]=matrix[i*depth+rk2];
		matrix[i*depth+rk2]=temp;
	}
}

double calculate_qp_objective(
							  long opt_n,
							  double *opt_g,double *opt_g0,double *alpha)
{
	double obj;
	long i,j;
	obj=0;  /* calculate objective  */
	for(i=0;i<opt_n;i++) {
		obj+=(opt_g0[i]*alpha[i]);
		obj+=(0.5*alpha[i]*alpha[i]*opt_g[i*opt_n+i]);
		for(j=0;j<i;j++) {
			obj+=(alpha[j]*alpha[i]*opt_g[j*opt_n+i]);
		}
	}
	return(obj);
}