SHOGUN v0.9.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * libocas.c: Implementation of the OCAS solver for training 00008 * linear SVM classifiers. 00009 * 00010 * Copyright (C) 2008 Vojtech Franc, xfrancv@cmp.felk.cvut.cz 00011 * Soeren Sonnenburg, soeren.sonnenburg@first.fraunhofer.de 00012 *-------------------------------------------------------------------- */ 00013 00014 #include <stdlib.h> 00015 #include <string.h> 00016 #include <math.h> 00017 #include <sys/time.h> 00018 #include <time.h> 00019 #include <stdio.h> 00020 #include <stdint.h> 00021 00022 #include "classifier/svm/libocas.h" 00023 #include "classifier/svm/libocas_common.h" 00024 #include "classifier/svm/libqp.h" 00025 00026 namespace shogun 00027 { 00028 #define LAMBDA 0.1 /* must be from (0,1> 1..means that OCAS becomes equivalent to CPA */ 00029 00030 static const uint32_t QPSolverMaxIter = 10000000; 00031 00032 static float64_t *H; 00033 static uint32_t BufSize; 00034 00035 /*---------------------------------------------------------------------- 00036 Returns pointer at i-th column of Hessian matrix. 00037 ----------------------------------------------------------------------*/ 00038 static const float64_t *get_col( uint32_t i) 00039 { 00040 return( &H[ BufSize*i ] ); 00041 } 00042 00043 /*---------------------------------------------------------------------- 00044 Returns time of the day in seconds. 00045 ----------------------------------------------------------------------*/ 00046 static float64_t get_time() 00047 { 00048 struct timeval tv; 00049 if (gettimeofday(&tv, NULL)==0) 00050 return tv.tv_sec+((float64_t)(tv.tv_usec))/1e6; 00051 else 00052 return 0.0; 00053 } 00054 00055 /*---------------------------------------------------------------------- 00056 Linear binary Ocas-SVM solver. 00057 ----------------------------------------------------------------------*/ 00058 ocas_return_value_T svm_ocas_solver( 00059 float64_t C, 00060 uint32_t nData, 00061 float64_t TolRel, 00062 float64_t TolAbs, 00063 float64_t QPBound, 00064 float64_t MaxTime, 00065 uint32_t _BufSize, 00066 uint8_t Method, 00067 void (*compute_W)(float64_t*, float64_t*, float64_t*, uint32_t, void*), 00068 float64_t (*update_W)(float64_t, void*), 00069 int (*add_new_cut)(float64_t*, uint32_t*, uint32_t, uint32_t, void*), 00070 int (*compute_output)(float64_t*, void* ), 00071 int (*sort)(float64_t*, float64_t*, uint32_t), 00072 void (*ocas_print)(ocas_return_value_T), 00073 void* user_data) 00074 { 00075 ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0}; 00076 float64_t *b, *alpha, *diag_H; 00077 float64_t *output, *old_output; 00078 float64_t xi, sq_norm_W, QPSolverTolRel, dot_prod_WoldW, sq_norm_oldW; 00079 float64_t A0, B0, GradVal, t, t1, t2, *Ci, *Bi, *hpf, *hpb; 00080 float64_t start_time, ocas_start_time; 00081 uint32_t cut_length; 00082 uint32_t i, *new_cut; 00083 uint32_t *I; 00084 uint8_t S = 1; 00085 libqp_state_T qp_exitflag; 00086 00087 ocas_start_time = get_time(); 00088 ocas.qp_solver_time = 0; 00089 ocas.output_time = 0; 00090 ocas.sort_time = 0; 00091 ocas.add_time = 0; 00092 ocas.w_time = 0; 00093 ocas.print_time = 0; 00094 float64_t gap; 00095 00096 BufSize = _BufSize; 00097 00098 QPSolverTolRel = TolRel*0.5; 00099 00100 H=NULL; 00101 b=NULL; 00102 alpha=NULL; 00103 new_cut=NULL; 00104 I=NULL; 00105 diag_H=NULL; 00106 output=NULL; 00107 old_output=NULL; 00108 hpf=NULL; 00109 hpb = NULL; 00110 Ci=NULL; 00111 Bi=NULL; 00112 00113 /* Hessian matrix contains dot product of normal vectors of selected cutting planes */ 00114 H = (float64_t*)LIBOCAS_CALLOC(BufSize*BufSize,sizeof(float64_t)); 00115 if(H == NULL) 00116 { 00117 ocas.exitflag=-2; 00118 goto cleanup; 00119 } 00120 00121 /* bias of cutting planes */ 00122 b = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 00123 if(b == NULL) 00124 { 00125 ocas.exitflag=-2; 00126 goto cleanup; 00127 } 00128 00129 alpha = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 00130 if(alpha == NULL) 00131 { 00132 ocas.exitflag=-2; 00133 goto cleanup; 00134 } 00135 00136 /* indices of examples which define a new cut */ 00137 new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,sizeof(uint32_t)); 00138 if(new_cut == NULL) 00139 { 00140 ocas.exitflag=-2; 00141 goto cleanup; 00142 } 00143 00144 I = (uint32_t*)LIBOCAS_CALLOC(BufSize,sizeof(uint32_t)); 00145 if(I == NULL) 00146 { 00147 ocas.exitflag=-2; 00148 goto cleanup; 00149 } 00150 00151 for(i=0; i< BufSize; i++) I[i] = 1; 00152 00153 diag_H = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 00154 if(diag_H == NULL) 00155 { 00156 ocas.exitflag=-2; 00157 goto cleanup; 00158 } 00159 00160 output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 00161 if(output == NULL) 00162 { 00163 ocas.exitflag=-2; 00164 goto cleanup; 00165 } 00166 00167 old_output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 00168 if(old_output == NULL) 00169 { 00170 ocas.exitflag=-2; 00171 goto cleanup; 00172 } 00173 00174 /* array of hinge points used in line-serach */ 00175 hpf = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpf[0])); 00176 if(hpf == NULL) 00177 { 00178 ocas.exitflag=-2; 00179 goto cleanup; 00180 } 00181 00182 hpb = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpb[0])); 00183 if(hpb == NULL) 00184 { 00185 ocas.exitflag=-2; 00186 goto cleanup; 00187 } 00188 00189 /* vectors Ci, Bi are used in the line search procedure */ 00190 Ci = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 00191 if(Ci == NULL) 00192 { 00193 ocas.exitflag=-2; 00194 goto cleanup; 00195 } 00196 00197 Bi = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 00198 if(Bi == NULL) 00199 { 00200 ocas.exitflag=-2; 00201 goto cleanup; 00202 } 00203 00204 ocas.nCutPlanes = 0; 00205 ocas.exitflag = 0; 00206 ocas.nIter = 0; 00207 00208 /* Compute initial value of Q_P assuming that W is zero vector.*/ 00209 sq_norm_W = 0; 00210 xi = nData; 00211 ocas.Q_P = 0.5*sq_norm_W + C*xi; 00212 ocas.Q_D = 0; 00213 00214 /* Compute the initial cutting plane */ 00215 cut_length = nData; 00216 for(i=0; i < nData; i++) 00217 new_cut[i] = i; 00218 00219 gap=(ocas.Q_P-ocas.Q_D)/CMath::abs(ocas.Q_P); 00220 SG_SABS_PROGRESS(gap, -CMath::log10(gap), -CMath::log10(1), -CMath::log10(TolRel), 6); 00221 00222 ocas.trn_err = nData; 00223 ocas.ocas_time = get_time() - ocas_start_time; 00224 /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, Q_P-Q_D/abs(Q_P)=%f\n", 00225 ocas.nIter,cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P)); 00226 */ 00227 ocas_print(ocas); 00228 00229 /* main loop */ 00230 while( ocas.exitflag == 0 ) 00231 { 00232 ocas.nIter++; 00233 00234 /* append a new cut to the buffer and update H */ 00235 b[ocas.nCutPlanes] = -(float64_t)cut_length; 00236 00237 start_time = get_time(); 00238 00239 if(add_new_cut( &H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, cut_length, ocas.nCutPlanes, user_data ) != 0) 00240 { 00241 ocas.exitflag=-2; 00242 goto cleanup; 00243 } 00244 00245 ocas.add_time += get_time() - start_time; 00246 00247 /* copy new added row: H(ocas.nCutPlanes,ocas.nCutPlanes,1:ocas.nCutPlanes-1) = H(1:ocas.nCutPlanes-1:ocas.nCutPlanes)' */ 00248 diag_H[ocas.nCutPlanes] = H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)]; 00249 for(i=0; i < ocas.nCutPlanes; i++) { 00250 H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] = H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)]; 00251 } 00252 00253 ocas.nCutPlanes++; 00254 00255 /* call inner QP solver */ 00256 start_time = get_time(); 00257 00258 qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &C, I, &S, alpha, 00259 ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0); 00260 00261 ocas.qp_exitflag = qp_exitflag.exitflag; 00262 00263 ocas.qp_solver_time += get_time() - start_time; 00264 ocas.Q_D = -qp_exitflag.QP; 00265 00266 ocas.nNZAlpha = 0; 00267 for(i=0; i < ocas.nCutPlanes; i++) { 00268 if( alpha[i] != 0) ocas.nNZAlpha++; 00269 } 00270 00271 sq_norm_oldW = sq_norm_W; 00272 start_time = get_time(); 00273 compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data ); 00274 ocas.w_time += get_time() - start_time; 00275 00276 /* select a new cut */ 00277 switch( Method ) 00278 { 00279 /* cutting plane algorithm implemented in SVMperf and BMRM */ 00280 case 0: 00281 00282 start_time = get_time(); 00283 if( compute_output( output, user_data ) != 0) 00284 { 00285 ocas.exitflag=-2; 00286 goto cleanup; 00287 } 00288 ocas.output_time += get_time()-start_time; 00289 00290 xi = 0; 00291 cut_length = 0; 00292 ocas.trn_err = 0; 00293 for(i=0; i < nData; i++) 00294 { 00295 if(output[i] <= 0) ocas.trn_err++; 00296 00297 if(output[i] <= 1) { 00298 xi += 1 - output[i]; 00299 new_cut[cut_length] = i; 00300 cut_length++; 00301 } 00302 } 00303 ocas.Q_P = 0.5*sq_norm_W + C*xi; 00304 00305 gap=(ocas.Q_P-ocas.Q_D)/CMath::abs(ocas.Q_P); 00306 SG_SABS_PROGRESS(gap, -CMath::log10(gap), -CMath::log10(1), -CMath::log10(TolRel), 6); 00307 /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, 1-Q_D/Q_P=%f, nza=%4d, err=%.2f%%, qpf=%d\n", 00308 ocas.nIter,cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P), 00309 ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag ); 00310 */ 00311 00312 start_time = get_time(); 00313 ocas_print(ocas); 00314 ocas.print_time += get_time() - start_time; 00315 00316 break; 00317 00318 00319 /* Ocas strategy */ 00320 case 1: 00321 00322 /* Linesearch */ 00323 A0 = sq_norm_W -2*dot_prod_WoldW + sq_norm_oldW; 00324 B0 = dot_prod_WoldW - sq_norm_oldW; 00325 00326 memcpy( old_output, output, sizeof(float64_t)*nData ); 00327 00328 start_time = get_time(); 00329 if( compute_output( output, user_data ) != 0) 00330 { 00331 ocas.exitflag=-2; 00332 goto cleanup; 00333 } 00334 ocas.output_time += get_time()-start_time; 00335 00336 uint32_t num_hp = 0; 00337 GradVal = B0; 00338 for(i=0; i< nData; i++) { 00339 00340 Ci[i] = C*(1-old_output[i]); 00341 Bi[i] = C*(old_output[i] - output[i]); 00342 00343 float64_t val; 00344 if(Bi[i] != 0) 00345 val = -Ci[i]/Bi[i]; 00346 else 00347 val = -LIBOCAS_PLUS_INF; 00348 00349 if (val>0) 00350 { 00351 /* hpi[num_hp] = i;*/ 00352 hpb[num_hp] = Bi[i]; 00353 hpf[num_hp] = val; 00354 num_hp++; 00355 } 00356 00357 if( (Bi[i] < 0 && val > 0) || (Bi[i] > 0 && val <= 0)) 00358 GradVal += Bi[i]; 00359 00360 } 00361 00362 t = 0; 00363 if( GradVal < 0 ) 00364 { 00365 start_time = get_time(); 00366 /* if( sort(hpf, hpi, num_hp) != 0)*/ 00367 if( sort(hpf, hpb, num_hp) != 0 ) 00368 { 00369 ocas.exitflag=-2; 00370 goto cleanup; 00371 } 00372 ocas.sort_time += get_time() - start_time; 00373 00374 float64_t t_new, GradVal_new; 00375 i = 0; 00376 while( GradVal < 0 && i < num_hp ) 00377 { 00378 t_new = hpf[i]; 00379 GradVal_new = GradVal + LIBOCAS_ABS(hpb[i]) + A0*(t_new-t); 00380 00381 if( GradVal_new >= 0 ) 00382 { 00383 t = t + GradVal*(t-t_new)/(GradVal_new - GradVal); 00384 } 00385 else 00386 { 00387 t = t_new; 00388 i++; 00389 } 00390 00391 GradVal = GradVal_new; 00392 } 00393 } 00394 00395 /* 00396 t = hpf[0] - 1; 00397 i = 0; 00398 GradVal = t*A0 + Bsum; 00399 while( GradVal < 0 && i < num_hp && hpf[i] < LIBOCAS_PLUS_INF ) { 00400 t = hpf[i]; 00401 Bsum = Bsum + LIBOCAS_ABS(Bi[hpi[i]]); 00402 GradVal = t*A0 + Bsum; 00403 i++; 00404 } 00405 */ 00406 t = LIBOCAS_MAX(t,0); /* just sanity check; t < 0 should not ocure */ 00407 00408 t1 = t; /* new (best so far) W */ 00409 t2 = t+LAMBDA*(1.0-t); /* new cutting plane */ 00410 /* t2 = t+(1.0-t)/10.0; */ 00411 00412 /* update W to be the best so far solution */ 00413 sq_norm_W = update_W( t1, user_data ); 00414 00415 /* select a new cut */ 00416 xi = 0; 00417 cut_length = 0; 00418 ocas.trn_err = 0; 00419 for(i=0; i < nData; i++ ) { 00420 00421 if( (old_output[i]*(1-t2) + t2*output[i]) <= 1 ) 00422 { 00423 new_cut[cut_length] = i; 00424 cut_length++; 00425 } 00426 00427 output[i] = old_output[i]*(1-t1) + t1*output[i]; 00428 00429 if( output[i] <= 1) xi += 1-output[i]; 00430 if( output[i] <= 0) ocas.trn_err++; 00431 00432 } 00433 00434 ocas.Q_P = 0.5*sq_norm_W + C*xi; 00435 00436 ocas.ocas_time = get_time() - ocas_start_time; 00437 00438 /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, 1-Q_D/Q_P=%f, nza=%4d, err=%.2f%%, qpf=%d\n", 00439 ocas.nIter, cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P), 00440 ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag ); 00441 */ 00442 00443 start_time = get_time(); 00444 ocas_print(ocas); 00445 ocas.print_time += get_time() - start_time; 00446 00447 break; 00448 } 00449 00450 /* Stopping conditions */ 00451 if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1; 00452 if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2; 00453 if( ocas.Q_P <= QPBound) ocas.exitflag = 3; 00454 if( ocas.ocas_time >= MaxTime) ocas.exitflag = 4; 00455 if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1; 00456 00457 } /* end of the main loop */ 00458 00459 cleanup: 00460 00461 LIBOCAS_FREE(H); 00462 LIBOCAS_FREE(b); 00463 LIBOCAS_FREE(alpha); 00464 LIBOCAS_FREE(new_cut); 00465 LIBOCAS_FREE(I); 00466 LIBOCAS_FREE(diag_H); 00467 LIBOCAS_FREE(output); 00468 LIBOCAS_FREE(old_output); 00469 LIBOCAS_FREE(hpf); 00470 /* LIBOCAS_FREE(hpi);*/ 00471 LIBOCAS_FREE(hpb); 00472 LIBOCAS_FREE(Ci); 00473 LIBOCAS_FREE(Bi); 00474 00475 ocas.ocas_time = get_time() - ocas_start_time; 00476 00477 return(ocas); 00478 } 00479 00480 00481 /*---------------------------------------------------------------------- 00482 Binary linear Ocas-SVM solver which allows using different C for each 00483 training example. 00484 ----------------------------------------------------------------------*/ 00485 ocas_return_value_T svm_ocas_solver_difC( 00486 float64_t *C, 00487 uint32_t nData, 00488 float64_t TolRel, 00489 float64_t TolAbs, 00490 float64_t QPBound, 00491 float64_t MaxTime, 00492 uint32_t _BufSize, 00493 uint8_t Method, 00494 void (*compute_W)(float64_t*, float64_t*, float64_t*, uint32_t, void*), 00495 float64_t (*update_W)(float64_t, void*), 00496 int (*add_new_cut)(float64_t*, uint32_t*, uint32_t, uint32_t, void*), 00497 int (*compute_output)(float64_t*, void* ), 00498 int (*sort)(float64_t*, float64_t*, uint32_t), 00499 void (*ocas_print)(ocas_return_value_T), 00500 void* user_data) 00501 { 00502 ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0}; 00503 float64_t *b, *alpha, *diag_H; 00504 float64_t *output, *old_output; 00505 float64_t xi, sq_norm_W, QPSolverTolRel, dot_prod_WoldW, sq_norm_oldW; 00506 float64_t A0, B0, GradVal, t, t1, t2, *Ci, *Bi, *hpf, *hpb; 00507 float64_t start_time, ocas_start_time; 00508 float64_t qp_b = 1.0; 00509 float64_t new_b; 00510 uint32_t cut_length; 00511 uint32_t i, *new_cut; 00512 uint32_t *I; 00513 uint8_t S = 1; 00514 libqp_state_T qp_exitflag; 00515 00516 ocas_start_time = get_time(); 00517 ocas.qp_solver_time = 0; 00518 ocas.output_time = 0; 00519 ocas.sort_time = 0; 00520 ocas.add_time = 0; 00521 ocas.w_time = 0; 00522 ocas.print_time = 0; 00523 00524 BufSize = _BufSize; 00525 00526 QPSolverTolRel = TolRel*0.5; 00527 00528 H=NULL; 00529 b=NULL; 00530 alpha=NULL; 00531 new_cut=NULL; 00532 I=NULL; 00533 diag_H=NULL; 00534 output=NULL; 00535 old_output=NULL; 00536 hpf=NULL; 00537 hpb = NULL; 00538 Ci=NULL; 00539 Bi=NULL; 00540 00541 /* Hessian matrix contains dot product of normal vectors of selected cutting planes */ 00542 H = (float64_t*)LIBOCAS_CALLOC(BufSize*BufSize,sizeof(float64_t)); 00543 if(H == NULL) 00544 { 00545 ocas.exitflag=-2; 00546 goto cleanup; 00547 } 00548 00549 /* bias of cutting planes */ 00550 b = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 00551 if(b == NULL) 00552 { 00553 ocas.exitflag=-2; 00554 goto cleanup; 00555 } 00556 00557 alpha = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 00558 if(alpha == NULL) 00559 { 00560 ocas.exitflag=-2; 00561 goto cleanup; 00562 } 00563 00564 /* indices of examples which define a new cut */ 00565 new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,sizeof(uint32_t)); 00566 if(new_cut == NULL) 00567 { 00568 ocas.exitflag=-2; 00569 goto cleanup; 00570 } 00571 00572 I = (uint32_t*)LIBOCAS_CALLOC(BufSize,sizeof(uint32_t)); 00573 if(I == NULL) 00574 { 00575 ocas.exitflag=-2; 00576 goto cleanup; 00577 } 00578 00579 for(i=0; i< BufSize; i++) I[i] = 1; 00580 00581 diag_H = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 00582 if(diag_H == NULL) 00583 { 00584 ocas.exitflag=-2; 00585 goto cleanup; 00586 } 00587 00588 output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 00589 if(output == NULL) 00590 { 00591 ocas.exitflag=-2; 00592 goto cleanup; 00593 } 00594 00595 old_output = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 00596 if(old_output == NULL) 00597 { 00598 ocas.exitflag=-2; 00599 goto cleanup; 00600 } 00601 00602 /* array of hinge points used in line-serach */ 00603 hpf = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpf[0])); 00604 if(hpf == NULL) 00605 { 00606 ocas.exitflag=-2; 00607 goto cleanup; 00608 } 00609 00610 hpb = (float64_t*) LIBOCAS_CALLOC(nData, sizeof(hpb[0])); 00611 if(hpb == NULL) 00612 { 00613 ocas.exitflag=-2; 00614 goto cleanup; 00615 } 00616 00617 /* vectors Ci, Bi are used in the line search procedure */ 00618 Ci = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 00619 if(Ci == NULL) 00620 { 00621 ocas.exitflag=-2; 00622 goto cleanup; 00623 } 00624 00625 Bi = (float64_t*)LIBOCAS_CALLOC(nData,sizeof(float64_t)); 00626 if(Bi == NULL) 00627 { 00628 ocas.exitflag=-2; 00629 goto cleanup; 00630 } 00631 00632 ocas.nCutPlanes = 0; 00633 ocas.exitflag = 0; 00634 ocas.nIter = 0; 00635 00636 /* Compute initial value of Q_P assuming that W is zero vector.*/ 00637 sq_norm_W = 0; 00638 xi = nData; 00639 /* ocas.Q_P = 0.5*sq_norm_W + C*xi;*/ 00640 ocas.Q_D = 0; 00641 00642 /* Compute the initial cutting plane */ 00643 cut_length = nData; 00644 new_b = 0; 00645 for(i=0; i < nData; i++) 00646 { 00647 new_cut[i] = i; 00648 new_b += C[i]; 00649 } 00650 00651 ocas.Q_P = 0.5*sq_norm_W + new_b; 00652 00653 00654 ocas.trn_err = nData; 00655 ocas.ocas_time = get_time() - ocas_start_time; 00656 /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, Q_P-Q_D/abs(Q_P)=%f\n", 00657 ocas.nIter,cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P)); 00658 */ 00659 ocas_print(ocas); 00660 00661 /* main loop */ 00662 while( ocas.exitflag == 0 ) 00663 { 00664 ocas.nIter++; 00665 00666 /* append a new cut to the buffer and update H */ 00667 /* b[ocas.nCutPlanes] = -(float64_t)cut_length*C;*/ 00668 b[ocas.nCutPlanes] = -new_b; 00669 00670 start_time = get_time(); 00671 00672 if(add_new_cut( &H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, cut_length, ocas.nCutPlanes, user_data ) != 0) 00673 { 00674 ocas.exitflag=-2; 00675 goto cleanup; 00676 } 00677 00678 ocas.add_time += get_time() - start_time; 00679 00680 /* copy new added row: H(ocas.nCutPlanes,ocas.nCutPlanes,1:ocas.nCutPlanes-1) = H(1:ocas.nCutPlanes-1:ocas.nCutPlanes)' */ 00681 diag_H[ocas.nCutPlanes] = H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)]; 00682 for(i=0; i < ocas.nCutPlanes; i++) { 00683 H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] = H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)]; 00684 } 00685 00686 ocas.nCutPlanes++; 00687 00688 /* call inner QP solver */ 00689 start_time = get_time(); 00690 00691 /* qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &C, I, &S, alpha,*/ 00692 /* ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);*/ 00693 qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &qp_b, I, &S, alpha, 00694 ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0); 00695 00696 ocas.qp_exitflag = qp_exitflag.exitflag; 00697 00698 ocas.qp_solver_time += get_time() - start_time; 00699 ocas.Q_D = -qp_exitflag.QP; 00700 00701 ocas.nNZAlpha = 0; 00702 for(i=0; i < ocas.nCutPlanes; i++) { 00703 if( alpha[i] != 0) ocas.nNZAlpha++; 00704 } 00705 00706 sq_norm_oldW = sq_norm_W; 00707 start_time = get_time(); 00708 compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data ); 00709 ocas.w_time += get_time() - start_time; 00710 00711 /* select a new cut */ 00712 switch( Method ) 00713 { 00714 /* cutting plane algorithm implemented in SVMperf and BMRM */ 00715 case 0: 00716 00717 start_time = get_time(); 00718 if( compute_output( output, user_data ) != 0) 00719 { 00720 ocas.exitflag=-2; 00721 goto cleanup; 00722 } 00723 ocas.output_time += get_time()-start_time; 00724 00725 xi = 0; 00726 cut_length = 0; 00727 ocas.trn_err = 0; 00728 new_b = 0; 00729 for(i=0; i < nData; i++) 00730 { 00731 if(output[i] <= 0) ocas.trn_err++; 00732 00733 /* if(output[i] <= 1) {*/ 00734 /* xi += 1 - output[i];*/ 00735 if(output[i] <= C[i]) { 00736 xi += C[i] - output[i]; 00737 new_cut[cut_length] = i; 00738 cut_length++; 00739 new_b += C[i]; 00740 } 00741 } 00742 /* ocas.Q_P = 0.5*sq_norm_W + C*xi;*/ 00743 ocas.Q_P = 0.5*sq_norm_W + xi; 00744 00745 ocas.ocas_time = get_time() - ocas_start_time; 00746 00747 /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, 1-Q_D/Q_P=%f, nza=%4d, err=%.2f%%, qpf=%d\n", 00748 ocas.nIter,cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P), 00749 ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag ); 00750 */ 00751 00752 start_time = get_time(); 00753 ocas_print(ocas); 00754 ocas.print_time += get_time() - start_time; 00755 00756 break; 00757 00758 00759 /* Ocas strategy */ 00760 case 1: 00761 00762 /* Linesearch */ 00763 A0 = sq_norm_W -2*dot_prod_WoldW + sq_norm_oldW; 00764 B0 = dot_prod_WoldW - sq_norm_oldW; 00765 00766 memcpy( old_output, output, sizeof(float64_t)*nData ); 00767 00768 start_time = get_time(); 00769 if( compute_output( output, user_data ) != 0) 00770 { 00771 ocas.exitflag=-2; 00772 goto cleanup; 00773 } 00774 ocas.output_time += get_time()-start_time; 00775 00776 uint32_t num_hp = 0; 00777 GradVal = B0; 00778 for(i=0; i< nData; i++) { 00779 00780 /* Ci[i] = C*(1-old_output[i]);*/ 00781 /* Bi[i] = C*(old_output[i] - output[i]);*/ 00782 Ci[i] = (C[i]-old_output[i]); 00783 Bi[i] = old_output[i] - output[i]; 00784 00785 float64_t val; 00786 if(Bi[i] != 0) 00787 val = -Ci[i]/Bi[i]; 00788 else 00789 val = -LIBOCAS_PLUS_INF; 00790 00791 if (val>0) 00792 { 00793 /* hpi[num_hp] = i;*/ 00794 hpb[num_hp] = Bi[i]; 00795 hpf[num_hp] = val; 00796 num_hp++; 00797 } 00798 00799 if( (Bi[i] < 0 && val > 0) || (Bi[i] > 0 && val <= 0)) 00800 GradVal += Bi[i]; 00801 00802 } 00803 00804 t = 0; 00805 if( GradVal < 0 ) 00806 { 00807 start_time = get_time(); 00808 /* if( sort(hpf, hpi, num_hp) != 0)*/ 00809 if( sort(hpf, hpb, num_hp) != 0 ) 00810 { 00811 ocas.exitflag=-2; 00812 goto cleanup; 00813 } 00814 ocas.sort_time += get_time() - start_time; 00815 00816 float64_t t_new, GradVal_new; 00817 i = 0; 00818 while( GradVal < 0 && i < num_hp ) 00819 { 00820 t_new = hpf[i]; 00821 GradVal_new = GradVal + LIBOCAS_ABS(hpb[i]) + A0*(t_new-t); 00822 00823 if( GradVal_new >= 0 ) 00824 { 00825 t = t + GradVal*(t-t_new)/(GradVal_new - GradVal); 00826 } 00827 else 00828 { 00829 t = t_new; 00830 i++; 00831 } 00832 00833 GradVal = GradVal_new; 00834 } 00835 } 00836 00837 /* 00838 t = hpf[0] - 1; 00839 i = 0; 00840 GradVal = t*A0 + Bsum; 00841 while( GradVal < 0 && i < num_hp && hpf[i] < LIBOCAS_PLUS_INF ) { 00842 t = hpf[i]; 00843 Bsum = Bsum + LIBOCAS_ABS(Bi[hpi[i]]); 00844 GradVal = t*A0 + Bsum; 00845 i++; 00846 } 00847 */ 00848 t = LIBOCAS_MAX(t,0); /* just sanity check; t < 0 should not ocure */ 00849 00850 t1 = t; /* new (best so far) W */ 00851 t2 = t+(1.0-t)*LAMBDA; /* new cutting plane */ 00852 /* t2 = t+(1.0-t)/10.0; new cutting plane */ 00853 00854 /* update W to be the best so far solution */ 00855 sq_norm_W = update_W( t1, user_data ); 00856 00857 /* select a new cut */ 00858 xi = 0; 00859 cut_length = 0; 00860 ocas.trn_err = 0; 00861 new_b = 0; 00862 for(i=0; i < nData; i++ ) { 00863 00864 /* if( (old_output[i]*(1-t2) + t2*output[i]) <= 1 ) */ 00865 if( (old_output[i]*(1-t2) + t2*output[i]) <= C[i] ) 00866 { 00867 new_cut[cut_length] = i; 00868 cut_length++; 00869 new_b += C[i]; 00870 } 00871 00872 output[i] = old_output[i]*(1-t1) + t1*output[i]; 00873 00874 /* if( output[i] <= 1) xi += 1-output[i];*/ 00875 if( output[i] <= C[i]) xi += C[i]-output[i]; 00876 if( output[i] <= 0) ocas.trn_err++; 00877 00878 } 00879 00880 /* ocas.Q_P = 0.5*sq_norm_W + C*xi;*/ 00881 ocas.Q_P = 0.5*sq_norm_W + xi; 00882 00883 ocas.ocas_time = get_time() - ocas_start_time; 00884 00885 /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, 1-Q_D/Q_P=%f, nza=%4d, err=%.2f%%, qpf=%d\n", 00886 ocas.nIter, cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P), 00887 ocas.nNZAlpha, 100*(float64_t)ocas.trn_err/(float64_t)nData, ocas.qp_exitflag ); 00888 */ 00889 00890 start_time = get_time(); 00891 ocas_print(ocas); 00892 ocas.print_time += get_time() - start_time; 00893 00894 break; 00895 } 00896 00897 /* Stopping conditions */ 00898 if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1; 00899 if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2; 00900 if( ocas.Q_P <= QPBound) ocas.exitflag = 3; 00901 if( ocas.ocas_time >= MaxTime) ocas.exitflag = 4; 00902 if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1; 00903 00904 } /* end of the main loop */ 00905 00906 cleanup: 00907 00908 LIBOCAS_FREE(H); 00909 LIBOCAS_FREE(b); 00910 LIBOCAS_FREE(alpha); 00911 LIBOCAS_FREE(new_cut); 00912 LIBOCAS_FREE(I); 00913 LIBOCAS_FREE(diag_H); 00914 LIBOCAS_FREE(output); 00915 LIBOCAS_FREE(old_output); 00916 LIBOCAS_FREE(hpf); 00917 /* LIBOCAS_FREE(hpi);*/ 00918 LIBOCAS_FREE(hpb); 00919 LIBOCAS_FREE(Ci); 00920 LIBOCAS_FREE(Bi); 00921 00922 ocas.ocas_time = get_time() - ocas_start_time; 00923 00924 return(ocas); 00925 } 00926 00927 00928 00929 /*---------------------------------------------------------------------- 00930 Multiclass SVM-Ocas solver 00931 ----------------------------------------------------------------------*/ 00932 00933 /* Helper function needed by the multi-class SVM linesearch. 00934 00935 - This function finds a simplified representation of a piece-wise linear function 00936 by splitting the domain into intervals and fining active terms for these intevals */ 00937 static void findactive(float64_t *Theta, float64_t *SortedA, uint32_t *nSortedA, float64_t *A, float64_t *B, int n, 00938 int (*sort)(float64_t*, float64_t*, uint32_t)) 00939 { 00940 float64_t tmp, theta; 00941 int32_t i, j, idx, idx2 = 0, start; 00942 00943 sort(A,B,n); 00944 00945 tmp = B[0]; 00946 idx = 0; 00947 i = 0; 00948 while( i < n-1 && A[i] == A[i+1]) 00949 { 00950 if( B[i+1] > B[idx] ) 00951 { 00952 idx = i+1; 00953 tmp = B[i+1]; 00954 } 00955 i++; 00956 } 00957 00958 (*nSortedA) = 1; 00959 SortedA[0] = A[idx]; 00960 00961 while(1) 00962 { 00963 start = idx + 1; 00964 while( start < n && A[idx] == A[start]) 00965 start++; 00966 00967 theta = LIBOCAS_PLUS_INF; 00968 for(j=start; j < n; j++) 00969 { 00970 tmp = (B[j] - B[idx])/(A[idx]-A[j]); 00971 if( tmp < theta) 00972 { 00973 theta = tmp; 00974 idx2 = j; 00975 } 00976 } 00977 00978 if( theta < LIBOCAS_PLUS_INF) 00979 { 00980 Theta[(*nSortedA) - 1] = theta; 00981 SortedA[(*nSortedA)] = A[idx2]; 00982 (*nSortedA)++; 00983 idx = idx2; 00984 } 00985 else 00986 return; 00987 } 00988 } 00989 00990 00991 /*---------------------------------------------------------------------- 00992 Multiclass linear OCAS-SVM solver. 00993 ----------------------------------------------------------------------*/ 00994 ocas_return_value_T msvm_ocas_solver( 00995 float64_t C, 00996 float64_t *data_y, 00997 uint32_t nY, 00998 uint32_t nData, 00999 float64_t TolRel, 01000 float64_t TolAbs, 01001 float64_t QPBound, 01002 float64_t MaxTime, 01003 uint32_t _BufSize, 01004 uint8_t Method, 01005 void (*compute_W)(float64_t*, float64_t*, float64_t*, uint32_t, void*), 01006 float64_t (*update_W)(float64_t, void*), 01007 int (*add_new_cut)(float64_t*, uint32_t*, uint32_t, void*), 01008 int (*compute_output)(float64_t*, void* ), 01009 int (*sort)(float64_t*, float64_t*, uint32_t), 01010 void (*ocas_print)(ocas_return_value_T), 01011 void* user_data) 01012 { 01013 ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0}; 01014 float64_t *b, *alpha, *diag_H; 01015 float64_t *output, *old_output; 01016 float64_t xi, sq_norm_W, QPSolverTolRel, QPSolverTolAbs, dot_prod_WoldW, sq_norm_oldW; 01017 float64_t A0, B0, t, t1, t2, R, tmp, element_b, x; 01018 float64_t *A, *B, *theta, *Theta, *sortedA, *Add; 01019 float64_t start_time, ocas_start_time, grad_sum, grad, min_x = 0, old_x, old_grad; 01020 uint32_t i, y, y2, ypred = 0, *new_cut, cnt1, cnt2, j, nSortedA, idx; 01021 uint32_t *I; 01022 uint8_t S = 1; 01023 libqp_state_T qp_exitflag; 01024 01025 ocas_start_time = get_time(); 01026 ocas.qp_solver_time = 0; 01027 ocas.output_time = 0; 01028 ocas.sort_time = 0; 01029 ocas.add_time = 0; 01030 ocas.w_time = 0; 01031 ocas.print_time = 0; 01032 01033 BufSize = _BufSize; 01034 01035 QPSolverTolRel = TolRel*0.5; 01036 QPSolverTolAbs = TolAbs*0.5; 01037 01038 H=NULL; 01039 b=NULL; 01040 alpha=NULL; 01041 new_cut=NULL; 01042 I=NULL; 01043 diag_H=NULL; 01044 output=NULL; 01045 old_output=NULL; 01046 A = NULL; 01047 B = NULL; 01048 theta = NULL; 01049 Theta = NULL; 01050 sortedA = NULL; 01051 Add = NULL; 01052 01053 /* Hessian matrix contains dot product of normal vectors of selected cutting planes */ 01054 H = (float64_t*)LIBOCAS_CALLOC(BufSize*BufSize,sizeof(float64_t)); 01055 if(H == NULL) 01056 { 01057 ocas.exitflag=-2; 01058 goto cleanup; 01059 } 01060 01061 /* bias of cutting planes */ 01062 b = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 01063 if(b == NULL) 01064 { 01065 ocas.exitflag=-2; 01066 goto cleanup; 01067 } 01068 01069 alpha = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 01070 if(alpha == NULL) 01071 { 01072 ocas.exitflag=-2; 01073 goto cleanup; 01074 } 01075 01076 /* indices of examples which define a new cut */ 01077 new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,sizeof(uint32_t)); 01078 if(new_cut == NULL) 01079 { 01080 ocas.exitflag=-2; 01081 goto cleanup; 01082 } 01083 01084 I = (uint32_t*)LIBOCAS_CALLOC(BufSize,sizeof(uint32_t)); 01085 if(I == NULL) 01086 { 01087 ocas.exitflag=-2; 01088 goto cleanup; 01089 } 01090 01091 for(i=0; i< BufSize; i++) 01092 I[i] = 1; 01093 01094 diag_H = (float64_t*)LIBOCAS_CALLOC(BufSize,sizeof(float64_t)); 01095 if(diag_H == NULL) 01096 { 01097 ocas.exitflag=-2; 01098 goto cleanup; 01099 } 01100 01101 output = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t)); 01102 if(output == NULL) 01103 { 01104 ocas.exitflag=-2; 01105 goto cleanup; 01106 } 01107 01108 old_output = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t)); 01109 if(old_output == NULL) 01110 { 01111 ocas.exitflag=-2; 01112 goto cleanup; 01113 } 01114 01115 /* auxciliary variables used in the linesearch */ 01116 A = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t)); 01117 if(A == NULL) 01118 { 01119 ocas.exitflag=-2; 01120 goto cleanup; 01121 } 01122 01123 B = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t)); 01124 if(B == NULL) 01125 { 01126 ocas.exitflag=-2; 01127 goto cleanup; 01128 } 01129 01130 theta = (float64_t*)LIBOCAS_CALLOC(nY,sizeof(float64_t)); 01131 if(theta == NULL) 01132 { 01133 ocas.exitflag=-2; 01134 goto cleanup; 01135 } 01136 01137 sortedA = (float64_t*)LIBOCAS_CALLOC(nY,sizeof(float64_t)); 01138 if(sortedA == NULL) 01139 { 01140 ocas.exitflag=-2; 01141 goto cleanup; 01142 } 01143 01144 Theta = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t)); 01145 if(Theta == NULL) 01146 { 01147 ocas.exitflag=-2; 01148 goto cleanup; 01149 } 01150 01151 Add = (float64_t*)LIBOCAS_CALLOC(nData*nY,sizeof(float64_t)); 01152 if(Add == NULL) 01153 { 01154 ocas.exitflag=-2; 01155 goto cleanup; 01156 } 01157 01158 /* Set initial values*/ 01159 ocas.nCutPlanes = 0; 01160 ocas.exitflag = 0; 01161 ocas.nIter = 0; 01162 ocas.Q_D = 0; 01163 ocas.trn_err = nData; 01164 R = (float64_t)nData; 01165 sq_norm_W = 0; 01166 element_b = (float64_t)nData; 01167 ocas.Q_P = 0.5*sq_norm_W + C*R; 01168 01169 /* initial cutting plane */ 01170 for(i=0; i < nData; i++) 01171 { 01172 y2 = (uint32_t)data_y[i]-1; 01173 01174 if(y2 > 0) 01175 new_cut[i] = 0; 01176 else 01177 new_cut[i] = 1; 01178 01179 } 01180 01181 ocas.ocas_time = get_time() - ocas_start_time; 01182 01183 start_time = get_time(); 01184 ocas_print(ocas); 01185 ocas.print_time += get_time() - start_time; 01186 01187 /* main loop of the OCAS */ 01188 while( ocas.exitflag == 0 ) 01189 { 01190 ocas.nIter++; 01191 01192 /* append a new cut to the buffer and update H */ 01193 b[ocas.nCutPlanes] = -(float64_t)element_b; 01194 01195 start_time = get_time(); 01196 01197 if(add_new_cut( &H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, ocas.nCutPlanes, user_data ) != 0) 01198 { 01199 ocas.exitflag=-2; 01200 goto cleanup; 01201 } 01202 01203 ocas.add_time += get_time() - start_time; 01204 01205 /* copy newly appended row: H(ocas.nCutPlanes,ocas.nCutPlanes,1:ocas.nCutPlanes-1) = H(1:ocas.nCutPlanes-1:ocas.nCutPlanes)' */ 01206 diag_H[ocas.nCutPlanes] = H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)]; 01207 for(i=0; i < ocas.nCutPlanes; i++) 01208 { 01209 H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] = H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)]; 01210 } 01211 01212 ocas.nCutPlanes++; 01213 01214 /* call inner QP solver */ 01215 start_time = get_time(); 01216 01217 qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &C, I, &S, alpha, 01218 ocas.nCutPlanes, QPSolverMaxIter, QPSolverTolAbs, QPSolverTolRel, -LIBOCAS_PLUS_INF,0); 01219 01220 ocas.qp_exitflag = qp_exitflag.exitflag; 01221 01222 ocas.qp_solver_time += get_time() - start_time; 01223 ocas.Q_D = -qp_exitflag.QP; 01224 01225 ocas.nNZAlpha = 0; 01226 for(i=0; i < ocas.nCutPlanes; i++) 01227 if( alpha[i] != 0) ocas.nNZAlpha++; 01228 01229 sq_norm_oldW = sq_norm_W; 01230 start_time = get_time(); 01231 compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data ); 01232 ocas.w_time += get_time() - start_time; 01233 01234 /* select a new cut */ 01235 switch( Method ) 01236 { 01237 /* cutting plane algorithm implemented in SVMperf and BMRM */ 01238 case 0: 01239 01240 start_time = get_time(); 01241 if( compute_output( output, user_data ) != 0) 01242 { 01243 ocas.exitflag=-2; 01244 goto cleanup; 01245 } 01246 ocas.output_time += get_time()-start_time; 01247 01248 /* the following loop computes: */ 01249 element_b = 0.0; /* element_b = R(old_W) - g'*old_W */ 01250 R = 0; /* R(W) = sum_i max_y ( [[y != y_i]] + (w_y- w_y_i)'*x_i ) */ 01251 ocas.trn_err = 0; /* trn_err = sum_i [[y != y_i ]] */ 01252 /* new_cut[i] = argmax_i ( [[y != y_i]] + (w_y- w_y_i)'*x_i ) */ 01253 for(i=0; i < nData; i++) 01254 { 01255 y2 = (uint32_t)data_y[i]-1; 01256 01257 for(xi=-LIBOCAS_PLUS_INF, y=0; y < nY; y++) 01258 { 01259 if(y2 != y && xi < output[LIBOCAS_INDEX(y,i,nY)]) 01260 { 01261 xi = output[LIBOCAS_INDEX(y,i,nY)]; 01262 ypred = y; 01263 } 01264 } 01265 01266 if(xi >= output[LIBOCAS_INDEX(y2,i,nY)]) 01267 ocas.trn_err ++; 01268 01269 xi = LIBOCAS_MAX(0,xi+1-output[LIBOCAS_INDEX(y2,i,nY)]); 01270 R += xi; 01271 if(xi > 0) 01272 { 01273 element_b++; 01274 new_cut[i] = ypred; 01275 } 01276 else 01277 new_cut[i] = y2; 01278 } 01279 01280 ocas.Q_P = 0.5*sq_norm_W + C*R; 01281 01282 ocas.ocas_time = get_time() - ocas_start_time; 01283 01284 start_time = get_time(); 01285 ocas_print(ocas); 01286 ocas.print_time += get_time() - start_time; 01287 01288 break; 01289 01290 /* The OCAS solver */ 01291 case 1: 01292 memcpy( old_output, output, sizeof(float64_t)*nData*nY ); 01293 01294 start_time = get_time(); 01295 if( compute_output( output, user_data ) != 0) 01296 { 01297 ocas.exitflag=-2; 01298 goto cleanup; 01299 } 01300 ocas.output_time += get_time()-start_time; 01301 01302 A0 = sq_norm_W - 2*dot_prod_WoldW + sq_norm_oldW; 01303 B0 = dot_prod_WoldW - sq_norm_oldW; 01304 01305 for(i=0; i < nData; i++) 01306 { 01307 y2 = (uint32_t)data_y[i]-1; 01308 01309 for(y=0; y < nY; y++) 01310 { 01311 A[LIBOCAS_INDEX(y,i,nY)] = C*(output[LIBOCAS_INDEX(y,i,nY)] - old_output[LIBOCAS_INDEX(y,i,nY)] 01312 + old_output[LIBOCAS_INDEX(y2,i,nY)] - output[LIBOCAS_INDEX(y2,i,nY)]); 01313 B[LIBOCAS_INDEX(y,i,nY)] = C*(old_output[LIBOCAS_INDEX(y,i,nY)] - old_output[LIBOCAS_INDEX(y2,i,nY)] 01314 + (float64_t)(y != y2)); 01315 } 01316 } 01317 01318 /* linesearch */ 01319 /* new_x = msvm_linesearch_mex(A0,B0,AA*C,BB*C);*/ 01320 01321 grad_sum = B0; 01322 cnt1 = 0; 01323 cnt2 = 0; 01324 for(i=0; i < nData; i++) 01325 { 01326 findactive(theta,sortedA,&nSortedA,&A[i*nY],&B[i*nY],nY,sort); 01327 01328 idx = 0; 01329 while( idx < nSortedA-1 && theta[idx] < 0 ) 01330 idx++; 01331 01332 grad_sum += sortedA[idx]; 01333 01334 for(j=idx; j < nSortedA-1; j++) 01335 { 01336 Theta[cnt1] = theta[j]; 01337 cnt1++; 01338 } 01339 01340 for(j=idx+1; j < nSortedA; j++) 01341 { 01342 Add[cnt2] = -sortedA[j-1]+sortedA[j]; 01343 cnt2++; 01344 } 01345 } 01346 01347 start_time = get_time(); 01348 sort(Theta,Add,cnt1); 01349 ocas.sort_time += get_time() - start_time; 01350 01351 grad = grad_sum; 01352 if(grad >= 0) 01353 { 01354 min_x = 0; 01355 } 01356 else 01357 { 01358 old_x = 0; 01359 old_grad = grad; 01360 01361 for(i=0; i < cnt1; i++) 01362 { 01363 x = Theta[i]; 01364 01365 grad = x*A0 + grad_sum; 01366 01367 if(grad >=0) 01368 { 01369 01370 min_x = (grad*old_x - old_grad*x)/(grad - old_grad); 01371 01372 break; 01373 } 01374 else 01375 { 01376 grad_sum = grad_sum + Add[i]; 01377 01378 grad = x*A0 + grad_sum; 01379 if( grad >= 0) 01380 { 01381 min_x = x; 01382 break; 01383 } 01384 } 01385 01386 old_grad = grad; 01387 old_x = x; 01388 } 01389 } 01390 /* end of the linesearch which outputs min_x */ 01391 01392 t = min_x; 01393 t1 = t; /* new (best so far) W */ 01394 t2 = t+(1.0-t)*LAMBDA; /* new cutting plane */ 01395 /* t2 = t+(1.0-t)/10.0; */ 01396 01397 /* update W to be the best so far solution */ 01398 sq_norm_W = update_W( t1, user_data ); 01399 01400 /* the following code computes a new cutting plane: */ 01401 element_b = 0.0; /* element_b = R(old_W) - g'*old_W */ 01402 /* new_cut[i] = argmax_i ( [[y != y_i]] + (w_y- w_y_i)'*x_i ) */ 01403 for(i=0; i < nData; i++) 01404 { 01405 y2 = (uint32_t)data_y[i]-1; 01406 01407 for(xi=-LIBOCAS_PLUS_INF, y=0; y < nY; y++) 01408 { 01409 tmp = old_output[LIBOCAS_INDEX(y,i,nY)]*(1-t2) + t2*output[LIBOCAS_INDEX(y,i,nY)]; 01410 if(y2 != y && xi < tmp) 01411 { 01412 xi = tmp; 01413 ypred = y; 01414 } 01415 } 01416 01417 tmp = old_output[LIBOCAS_INDEX(y2,i,nY)]*(1-t2) + t2*output[LIBOCAS_INDEX(y2,i,nY)]; 01418 xi = LIBOCAS_MAX(0,xi+1-tmp); 01419 if(xi > 0) 01420 { 01421 element_b++; 01422 new_cut[i] = ypred; 01423 } 01424 else 01425 new_cut[i] = y2; 01426 } 01427 01428 /* compute Risk, class. error and update outputs to correspond to the new W */ 01429 ocas.trn_err = 0; /* trn_err = sum_i [[y != y_i ]] */ 01430 R = 0; 01431 for(i=0; i < nData; i++) 01432 { 01433 y2 = (uint32_t)data_y[i]-1; 01434 01435 for(tmp=-LIBOCAS_PLUS_INF, y=0; y < nY; y++) 01436 { 01437 output[LIBOCAS_INDEX(y,i,nY)] = old_output[LIBOCAS_INDEX(y,i,nY)]*(1-t1) + t1*output[LIBOCAS_INDEX(y,i,nY)]; 01438 01439 if(y2 != y && tmp < output[LIBOCAS_INDEX(y,i,nY)]) 01440 { 01441 ypred = y; 01442 tmp = output[LIBOCAS_INDEX(y,i,nY)]; 01443 } 01444 } 01445 01446 R += LIBOCAS_MAX(0,1+tmp - output[LIBOCAS_INDEX(y2,i,nY)]); 01447 if( tmp >= output[LIBOCAS_INDEX(y2,i,nY)]) 01448 ocas.trn_err ++; 01449 } 01450 01451 ocas.Q_P = 0.5*sq_norm_W + C*R; 01452 01453 01454 /* get time and print status */ 01455 ocas.ocas_time = get_time() - ocas_start_time; 01456 01457 start_time = get_time(); 01458 ocas_print(ocas); 01459 ocas.print_time += get_time() - start_time; 01460 01461 break; 01462 01463 } 01464 01465 /* Stopping conditions */ 01466 if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1; 01467 if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2; 01468 if( ocas.Q_P <= QPBound) ocas.exitflag = 3; 01469 if( ocas.ocas_time >= MaxTime) ocas.exitflag = 4; 01470 if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1; 01471 01472 } /* end of the main loop */ 01473 01474 cleanup: 01475 01476 LIBOCAS_FREE(H); 01477 LIBOCAS_FREE(b); 01478 LIBOCAS_FREE(alpha); 01479 LIBOCAS_FREE(new_cut); 01480 LIBOCAS_FREE(I); 01481 LIBOCAS_FREE(diag_H); 01482 LIBOCAS_FREE(output); 01483 LIBOCAS_FREE(old_output); 01484 LIBOCAS_FREE(A); 01485 LIBOCAS_FREE(B); 01486 LIBOCAS_FREE(theta); 01487 LIBOCAS_FREE(Theta); 01488 LIBOCAS_FREE(sortedA); 01489 LIBOCAS_FREE(Add); 01490 01491 ocas.ocas_time = get_time() - ocas_start_time; 01492 01493 return(ocas); 01494 } 01495 }