wolffd@0: /***********************************************************************/
wolffd@0: /*                                                                     */
wolffd@0: /*   svm_loqo.c                                                        */
wolffd@0: /*                                                                     */
wolffd@0: /*   Interface to the PR_LOQO optimization package for SVM.            */
wolffd@0: /*                                                                     */
wolffd@0: /*   Author: Thorsten Joachims                                         */
wolffd@0: /*   Date: 19.07.99                                                    */
wolffd@0: /*                                                                     */
wolffd@0: /*   Copyright (c) 1999  Universitaet Dortmund - All rights reserved   */
wolffd@0: /*                                                                     */
wolffd@0: /*   This software is available for non-commercial use only. It must   */
wolffd@0: /*   not be modified and distributed without prior permission of the   */
wolffd@0: /*   author. The author is not responsible for implications from the   */
wolffd@0: /*   use of this software.                                             */
wolffd@0: /*                                                                     */
wolffd@0: /***********************************************************************/
wolffd@0: 
wolffd@0: # include <math.h>
wolffd@0: # include "pr_loqo/pr_loqo.h"
wolffd@0: # include "svm_common.h"
wolffd@0: 
wolffd@0: /* Common Block Declarations */
wolffd@0: 
wolffd@0: long verbosity;
wolffd@0: 
wolffd@0: /* /////////////////////////////////////////////////////////////// */
wolffd@0: 
wolffd@0: # define DEF_PRECISION_LINEAR    1E-8
wolffd@0: # define DEF_PRECISION_NONLINEAR 1E-14
wolffd@0: 
wolffd@0: double *optimize_qp();
wolffd@0: double *primal=0,*dual=0;
wolffd@0: double init_margin=0.15;
wolffd@0: long   init_iter=500,precision_violations=0;
wolffd@0: double model_b;
wolffd@0: double opt_precision=DEF_PRECISION_LINEAR;
wolffd@0: 
wolffd@0: /* /////////////////////////////////////////////////////////////// */
wolffd@0: 
wolffd@0: void *my_malloc();
wolffd@0: 
wolffd@0: double *optimize_qp(qp,epsilon_crit,nx,threshold,learn_parm)
wolffd@0: QP *qp;
wolffd@0: double *epsilon_crit;
wolffd@0: long nx; /* Maximum number of variables in QP */
wolffd@0: double *threshold;
wolffd@0: LEARN_PARM *learn_parm;
wolffd@0: /* start the optimizer and return the optimal values */
wolffd@0: {
wolffd@0:   register long i,j,result;
wolffd@0:   double margin,obj_before,obj_after;
wolffd@0:   double sigdig,dist,epsilon_loqo;
wolffd@0:   int iter;
wolffd@0:  
wolffd@0:   if(!primal) { /* allocate memory at first call */
wolffd@0:     primal=(double *)my_malloc(sizeof(double)*nx*3);
wolffd@0:     dual=(double *)my_malloc(sizeof(double)*(nx*2+1));
wolffd@0:   }
wolffd@0:   
wolffd@0:   if(verbosity>=4) { /* really verbose */
wolffd@0:     printf("\n\n");
wolffd@0:     for(i=0;i<qp->opt_n;i++) {
wolffd@0:       printf("%f: ",qp->opt_g0[i]);
wolffd@0:       for(j=0;j<qp->opt_n;j++) {
wolffd@0: 	printf("%f ",qp->opt_g[i*qp->opt_n+j]);
wolffd@0:       }
wolffd@0:       printf(": a%ld=%.10f < %f",i,qp->opt_xinit[i],qp->opt_up[i]);
wolffd@0:       printf(": y=%f\n",qp->opt_ce[i]);
wolffd@0:     }
wolffd@0:     for(j=0;j<qp->opt_m;j++) {
wolffd@0:       printf("EQ-%ld: %f*a0",j,qp->opt_ce[j]);
wolffd@0:       for(i=1;i<qp->opt_n;i++) {
wolffd@0: 	printf(" + %f*a%ld",qp->opt_ce[i],i);
wolffd@0:       }
wolffd@0:       printf(" = %f\n\n",-qp->opt_ce0[0]);
wolffd@0:     }
wolffd@0: }
wolffd@0: 
wolffd@0:   obj_before=0; /* calculate objective before optimization */
wolffd@0:   for(i=0;i<qp->opt_n;i++) {
wolffd@0:     obj_before+=(qp->opt_g0[i]*qp->opt_xinit[i]);
wolffd@0:     obj_before+=(0.5*qp->opt_xinit[i]*qp->opt_xinit[i]*qp->opt_g[i*qp->opt_n+i]);
wolffd@0:     for(j=0;j<i;j++) {
wolffd@0:       obj_before+=(qp->opt_xinit[j]*qp->opt_xinit[i]*qp->opt_g[j*qp->opt_n+i]);
wolffd@0:     }
wolffd@0:   }
wolffd@0: 
wolffd@0:   result=STILL_RUNNING;
wolffd@0:   qp->opt_ce0[0]*=(-1.0);
wolffd@0:   /* Run pr_loqo. If a run fails, try again with parameters which lead */
wolffd@0:   /* to a slower, but more robust setting. */
wolffd@0:   for(margin=init_margin,iter=init_iter;
wolffd@0:       (margin<=0.9999999) && (result!=OPTIMAL_SOLUTION);) {
wolffd@0:     sigdig=-log10(opt_precision);
wolffd@0: 
wolffd@0:     result=pr_loqo((int)qp->opt_n,(int)qp->opt_m,
wolffd@0: 		   (double *)qp->opt_g0,(double *)qp->opt_g,
wolffd@0: 		   (double *)qp->opt_ce,(double *)qp->opt_ce0,
wolffd@0: 		   (double *)qp->opt_low,(double *)qp->opt_up,
wolffd@0: 		   (double *)primal,(double *)dual, 
wolffd@0: 		   (int)(verbosity-2),
wolffd@0: 		   (double)sigdig,(int)iter, 
wolffd@0: 		   (double)margin,(double)(qp->opt_up[0])/4.0,(int)0);
wolffd@0: 
wolffd@0:     if(isnan(dual[0])) {     /* check for choldc problem */
wolffd@0:       if(verbosity>=2) {
wolffd@0: 	printf("NOTICE: Restarting PR_LOQO with more conservative parameters.\n");
wolffd@0:       }
wolffd@0:       if(init_margin<0.80) { /* become more conservative in general */
wolffd@0: 	init_margin=(4.0*margin+1.0)/5.0;
wolffd@0:       }
wolffd@0:       margin=(margin+1.0)/2.0;
wolffd@0:       (opt_precision)*=10.0;   /* reduce precision */
wolffd@0:       if(verbosity>=2) {
wolffd@0: 	printf("NOTICE: Reducing precision of PR_LOQO.\n");
wolffd@0:       }
wolffd@0:     }
wolffd@0:     else if(result!=OPTIMAL_SOLUTION) {
wolffd@0:       iter+=2000; 
wolffd@0:       init_iter+=10;
wolffd@0:       (opt_precision)*=10.0;   /* reduce precision */
wolffd@0:       if(verbosity>=2) {
wolffd@0: 	printf("NOTICE: Reducing precision of PR_LOQO due to (%ld).\n",result);
wolffd@0:       }      
wolffd@0:     }
wolffd@0:   }
wolffd@0: 
wolffd@0:   if(qp->opt_m)         /* Thanks to Alex Smola for this hint */
wolffd@0:     model_b=dual[0];
wolffd@0:   else
wolffd@0:     model_b=0;
wolffd@0: 
wolffd@0:   /* Check the precision of the alphas. If results of current optimization */
wolffd@0:   /* violate KT-Conditions, relax the epsilon on the bounds on alphas. */
wolffd@0:   epsilon_loqo=1E-10;
wolffd@0:   for(i=0;i<qp->opt_n;i++) {
wolffd@0:     dist=-model_b*qp->opt_ce[i]; 
wolffd@0:     dist+=(qp->opt_g0[i]+1.0);
wolffd@0:     for(j=0;j<i;j++) {
wolffd@0:       dist+=(primal[j]*qp->opt_g[j*qp->opt_n+i]);
wolffd@0:     }
wolffd@0:     for(j=i;j<qp->opt_n;j++) {
wolffd@0:       dist+=(primal[j]*qp->opt_g[i*qp->opt_n+j]);
wolffd@0:     }
wolffd@0:     /*  printf("LOQO: a[%d]=%f, dist=%f, b=%f\n",i,primal[i],dist,dual[0]); */
wolffd@0:     if((primal[i]<(qp->opt_up[i]-epsilon_loqo)) && (dist < (1.0-(*epsilon_crit)))) {
wolffd@0:       epsilon_loqo=(qp->opt_up[i]-primal[i])*2.0;
wolffd@0:     }
wolffd@0:     else if((primal[i]>(0+epsilon_loqo)) && (dist > (1.0+(*epsilon_crit)))) {
wolffd@0:       epsilon_loqo=primal[i]*2.0;
wolffd@0:     }
wolffd@0:   }
wolffd@0: 
wolffd@0:   for(i=0;i<qp->opt_n;i++) {  /* clip alphas to bounds */
wolffd@0:     if(primal[i]<=(0+epsilon_loqo)) {
wolffd@0:       primal[i]=0;
wolffd@0:     }
wolffd@0:     else if(primal[i]>=(qp->opt_up[i]-epsilon_loqo)) {
wolffd@0:       primal[i]=qp->opt_up[i];
wolffd@0:     }
wolffd@0:   }
wolffd@0: 
wolffd@0:   obj_after=0;  /* calculate objective after optimization */
wolffd@0:   for(i=0;i<qp->opt_n;i++) {
wolffd@0:     obj_after+=(qp->opt_g0[i]*primal[i]);
wolffd@0:     obj_after+=(0.5*primal[i]*primal[i]*qp->opt_g[i*qp->opt_n+i]);
wolffd@0:     for(j=0;j<i;j++) {
wolffd@0:       obj_after+=(primal[j]*primal[i]*qp->opt_g[j*qp->opt_n+i]);
wolffd@0:     }
wolffd@0:   }
wolffd@0: 
wolffd@0:   /* if optimizer returned NAN values, reset and retry with smaller */
wolffd@0:   /* working set. */
wolffd@0:   if(isnan(obj_after) || isnan(model_b)) {
wolffd@0:     for(i=0;i<qp->opt_n;i++) {
wolffd@0:       primal[i]=qp->opt_xinit[i];
wolffd@0:     }     
wolffd@0:     model_b=0;
wolffd@0:     if(learn_parm->svm_maxqpsize>2) {
wolffd@0:       learn_parm->svm_maxqpsize--;  /* decrease size of qp-subproblems */
wolffd@0:     }
wolffd@0:   }
wolffd@0: 
wolffd@0:   if(obj_after >= obj_before) { /* check whether there was progress */
wolffd@0:     (opt_precision)/=100.0;
wolffd@0:     precision_violations++;
wolffd@0:     if(verbosity>=2) {
wolffd@0:       printf("NOTICE: Increasing Precision of PR_LOQO.\n");
wolffd@0:     }
wolffd@0:   }
wolffd@0: 
wolffd@0:   if(precision_violations > 500) { 
wolffd@0:     (*epsilon_crit)*=10.0;
wolffd@0:     precision_violations=0;
wolffd@0:     if(verbosity>=1) {
wolffd@0:       printf("\nWARNING: Relaxing epsilon on KT-Conditions.\n");
wolffd@0:     }
wolffd@0:   }	  
wolffd@0: 
wolffd@0:   (*threshold)=model_b;
wolffd@0: 
wolffd@0:   if(result!=OPTIMAL_SOLUTION) {
wolffd@0:     printf("\nERROR: PR_LOQO did not converge. \n");
wolffd@0:     return(qp->opt_xinit);
wolffd@0:   }
wolffd@0:   else {
wolffd@0:     return(primal);
wolffd@0:   }
wolffd@0: }
wolffd@0: