1 /***********************************************************************/
2 /* */
3 /* svm_loqo.c */
4 /* */
5 /* Interface to the PR_LOQO optimization package for SVM. */
6 /* */
7 /* Author: Thorsten Joachims */
8 /* Date: 19.07.99 */
9 /* */
10 /* Copyright (c) 1999 Universitaet Dortmund - All rights reserved */
11 /* */
12 /* This software is available for non-commercial use only. It must */
13 /* not be modified and distributed without prior permission of the */
14 /* author. The author is not responsible for implications from the */
15 /* use of this software. */
16 /* */
17 /***********************************************************************/
18
19 # include <math.h>
20 # include "pr_loqo/pr_loqo.h"
21 # include "svm_common.h"
22
23 /* Common Block Declarations */
24
25 long verbosity;
26
27 /* /////////////////////////////////////////////////////////////// */
28
29 # define DEF_PRECISION_LINEAR 1E-8
30 # define DEF_PRECISION_NONLINEAR 1E-14
31
32 double *optimize_qp();
33 double *primal=0,*dual=0;
34 double init_margin=0.15;
35 long init_iter=500,precision_violations=0;
36 double model_b;
37 double opt_precision=DEF_PRECISION_LINEAR;
38
39 /* /////////////////////////////////////////////////////////////// */
40
41 void *my_malloc();
42
optimize_qp(qp,epsilon_crit,nx,threshold,learn_parm)43 double *optimize_qp(qp,epsilon_crit,nx,threshold,learn_parm)
44 QP *qp;
45 double *epsilon_crit;
46 long nx; /* Maximum number of variables in QP */
47 double *threshold;
48 LEARN_PARM *learn_parm;
49 /* start the optimizer and return the optimal values */
50 {
51 register long i,j,result;
52 double margin,obj_before,obj_after;
53 double sigdig,dist,epsilon_loqo;
54 int iter;
55
56 if(!primal) { /* allocate memory at first call */
57 primal=(double *)my_malloc(sizeof(double)*nx*3);
58 dual=(double *)my_malloc(sizeof(double)*(nx*2+1));
59 }
60
61 if(verbosity>=4) { /* really verbose */
62 printf("\n\n");
63 for(i=0;i<qp->opt_n;i++) {
64 printf("%f: ",qp->opt_g0[i]);
65 for(j=0;j<qp->opt_n;j++) {
66 printf("%f ",qp->opt_g[i*qp->opt_n+j]);
67 }
68 printf(": a%ld=%.10f < %f",i,qp->opt_xinit[i],qp->opt_up[i]);
69 printf(": y=%f\n",qp->opt_ce[i]);
70 }
71 for(j=0;j<qp->opt_m;j++) {
72 printf("EQ-%ld: %f*a0",j,qp->opt_ce[j]);
73 for(i=1;i<qp->opt_n;i++) {
74 printf(" + %f*a%ld",qp->opt_ce[i],i);
75 }
76 printf(" = %f\n\n",-qp->opt_ce0[0]);
77 }
78 }
79
80 obj_before=0; /* calculate objective before optimization */
81 for(i=0;i<qp->opt_n;i++) {
82 obj_before+=(qp->opt_g0[i]*qp->opt_xinit[i]);
83 obj_before+=(0.5*qp->opt_xinit[i]*qp->opt_xinit[i]*qp->opt_g[i*qp->opt_n+i]);
84 for(j=0;j<i;j++) {
85 obj_before+=(qp->opt_xinit[j]*qp->opt_xinit[i]*qp->opt_g[j*qp->opt_n+i]);
86 }
87 }
88
89 result=STILL_RUNNING;
90 qp->opt_ce0[0]*=(-1.0);
91 /* Run pr_loqo. If a run fails, try again with parameters which lead */
92 /* to a slower, but more robust setting. */
93 for(margin=init_margin,iter=init_iter;
94 (margin<=0.9999999) && (result!=OPTIMAL_SOLUTION);) {
95 sigdig=-log10(opt_precision);
96
97 result=pr_loqo((int)qp->opt_n,(int)qp->opt_m,
98 (double *)qp->opt_g0,(double *)qp->opt_g,
99 (double *)qp->opt_ce,(double *)qp->opt_ce0,
100 (double *)qp->opt_low,(double *)qp->opt_up,
101 (double *)primal,(double *)dual,
102 (int)(verbosity-2),
103 (double)sigdig,(int)iter,
104 (double)margin,(double)(qp->opt_up[0])/4.0,(int)0);
105
106 if(isnan(dual[0])) { /* check for choldc problem */
107 if(verbosity>=2) {
108 printf("NOTICE: Restarting PR_LOQO with more conservative parameters.\n");
109 }
110 if(init_margin<0.80) { /* become more conservative in general */
111 init_margin=(4.0*margin+1.0)/5.0;
112 }
113 margin=(margin+1.0)/2.0;
114 (opt_precision)*=10.0; /* reduce precision */
115 if(verbosity>=2) {
116 printf("NOTICE: Reducing precision of PR_LOQO.\n");
117 }
118 }
119 else if(result!=OPTIMAL_SOLUTION) {
120 iter+=2000;
121 init_iter+=10;
122 (opt_precision)*=10.0; /* reduce precision */
123 if(verbosity>=2) {
124 printf("NOTICE: Reducing precision of PR_LOQO due to (%ld).\n",result);
125 }
126 }
127 }
128
129 if(qp->opt_m) /* Thanks to Alex Smola for this hint */
130 model_b=dual[0];
131 else
132 model_b=0;
133
134 /* Check the precision of the alphas. If results of current optimization */
135 /* violate KT-Conditions, relax the epsilon on the bounds on alphas. */
136 epsilon_loqo=1E-10;
137 for(i=0;i<qp->opt_n;i++) {
138 dist=-model_b*qp->opt_ce[i];
139 dist+=(qp->opt_g0[i]+1.0);
140 for(j=0;j<i;j++) {
141 dist+=(primal[j]*qp->opt_g[j*qp->opt_n+i]);
142 }
143 for(j=i;j<qp->opt_n;j++) {
144 dist+=(primal[j]*qp->opt_g[i*qp->opt_n+j]);
145 }
146 /* printf("LOQO: a[%d]=%f, dist=%f, b=%f\n",i,primal[i],dist,dual[0]); */
147 if((primal[i]<(qp->opt_up[i]-epsilon_loqo)) && (dist < (1.0-(*epsilon_crit)))) {
148 epsilon_loqo=(qp->opt_up[i]-primal[i])*2.0;
149 }
150 else if((primal[i]>(0+epsilon_loqo)) && (dist > (1.0+(*epsilon_crit)))) {
151 epsilon_loqo=primal[i]*2.0;
152 }
153 }
154
155 for(i=0;i<qp->opt_n;i++) { /* clip alphas to bounds */
156 if(primal[i]<=(0+epsilon_loqo)) {
157 primal[i]=0;
158 }
159 else if(primal[i]>=(qp->opt_up[i]-epsilon_loqo)) {
160 primal[i]=qp->opt_up[i];
161 }
162 }
163
164 obj_after=0; /* calculate objective after optimization */
165 for(i=0;i<qp->opt_n;i++) {
166 obj_after+=(qp->opt_g0[i]*primal[i]);
167 obj_after+=(0.5*primal[i]*primal[i]*qp->opt_g[i*qp->opt_n+i]);
168 for(j=0;j<i;j++) {
169 obj_after+=(primal[j]*primal[i]*qp->opt_g[j*qp->opt_n+i]);
170 }
171 }
172
173 /* if optimizer returned NAN values, reset and retry with smaller */
174 /* working set. */
175 if(isnan(obj_after) || isnan(model_b)) {
176 for(i=0;i<qp->opt_n;i++) {
177 primal[i]=qp->opt_xinit[i];
178 }
179 model_b=0;
180 if(learn_parm->svm_maxqpsize>2) {
181 learn_parm->svm_maxqpsize--; /* decrease size of qp-subproblems */
182 }
183 }
184
185 if(obj_after >= obj_before) { /* check whether there was progress */
186 (opt_precision)/=100.0;
187 precision_violations++;
188 if(verbosity>=2) {
189 printf("NOTICE: Increasing Precision of PR_LOQO.\n");
190 }
191 }
192
193 if(precision_violations > 500) {
194 (*epsilon_crit)*=10.0;
195 precision_violations=0;
196 if(verbosity>=1) {
197 printf("\nWARNING: Relaxing epsilon on KT-Conditions.\n");
198 }
199 }
200
201 (*threshold)=model_b;
202
203 if(result!=OPTIMAL_SOLUTION) {
204 printf("\nERROR: PR_LOQO did not converge. \n");
205 return(qp->opt_xinit);
206 }
207 else {
208 return(primal);
209 }
210 }
211
212