SHOGUN
v3.2.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 * libp3bm.h: Implementation of the Proximal Point P-BMRM solver for SO training 00008 * 00009 * Copyright (C) 2012 Michal Uricar, uricamic@cmp.felk.cvut.cz 00010 * 00011 * Implementation of the Proximal Point P-BMRM 00012 *--------------------------------------------------------------------- */ 00013 00014 #include <shogun/structure/libp3bm.h> 00015 #include <shogun/lib/external/libqp.h> 00016 #include <shogun/lib/Time.h> 00017 00018 namespace shogun 00019 { 00020 static const uint32_t QPSolverMaxIter=0xFFFFFFFF; 00021 static const float64_t epsilon=0.0; 00022 00023 static float64_t *H, *H2; 00024 static uint32_t BufSize; 00025 00026 /*---------------------------------------------------------------------- 00027 Returns pointer at i-th column of Hessian matrix. 00028 ----------------------------------------------------------------------*/ 00029 static const float64_t *get_col( uint32_t i) 00030 { 00031 return( &H2[ BufSize*i ] ); 00032 } 00033 00034 BmrmStatistics svm_ppbm_solver( 00035 CDualLibQPBMSOSVM *machine, 00036 float64_t* W, 00037 float64_t TolRel, 00038 float64_t TolAbs, 00039 float64_t _lambda, 00040 uint32_t _BufSize, 00041 bool cleanICP, 00042 uint32_t cleanAfter, 00043 float64_t K, 00044 uint32_t Tmax, 00045 bool verbose) 00046 { 00047 BmrmStatistics ppbmrm; 00048 libqp_state_T qp_exitflag={0, 0, 0, 0}, qp_exitflag_good={0, 0, 0, 0}; 00049 float64_t *b, *b2, *beta, *beta_good, *beta_start, *diag_H, *diag_H2; 00050 float64_t R, *subgrad, *A, QPSolverTolRel, C=1.0; 00051 float64_t *prevW, *wt, alpha, alpha_start, alpha_good=0.0, Fd_alpha0=0.0; 00052 float64_t lastFp, wdist, gamma=0.0; 00053 floatmax_t rsum, sq_norm_W, sq_norm_Wdiff, sq_norm_prevW, eps; 00054 uint32_t *I, *I2, *I_start, *I_good; 00055 uint8_t S=1; 00056 CStructuredModel* model=machine->get_model(); 00057 uint32_t nDim=model->get_dim(); 00058 CSOSVMHelper* helper = NULL; 00059 uint32_t qp_cnt=0; 00060 bmrm_ll *CPList_head, *CPList_tail, *cp_ptr, *cp_ptr2, *cp_list=NULL; 00061 float64_t *A_1=NULL, *A_2=NULL; 00062 bool *map=NULL, tuneAlpha=true, flag=true, alphaChanged=false, isThereGoodSolution=false; 00063 00064 CTime ttime; 00065 float64_t tstart, tstop; 00066 00067 00068 tstart=ttime.cur_time_diff(false); 00069 00070 BufSize=_BufSize; 00071 QPSolverTolRel=1e-9; 00072 00073 H=NULL; 00074 b=NULL; 00075 beta=NULL; 00076 A=NULL; 00077 subgrad=NULL; 00078 diag_H=NULL; 00079 I=NULL; 00080 prevW=NULL; 00081 wt=NULL; 00082 diag_H2=NULL; 00083 b2=NULL; 00084 I2=NULL; 00085 H2=NULL; 00086 I_good=NULL; 00087 I_start=NULL; 00088 beta_start=NULL; 00089 beta_good=NULL; 00090 00091 alpha=0.0; 00092 00093 H= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, float64_t); 00094 00095 A= (float64_t*) LIBBMRM_CALLOC(nDim*BufSize, float64_t); 00096 00097 b= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t); 00098 00099 beta= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t); 00100 00101 subgrad= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t); 00102 00103 diag_H= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t); 00104 00105 I= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t); 00106 00107 /* structure for maintaining inactive CPs info */ 00108 ICP_stats icp_stats; 00109 icp_stats.maxCPs = BufSize; 00110 icp_stats.ICPcounter= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t); 00111 icp_stats.ICPs= (float64_t**) LIBBMRM_CALLOC(BufSize, float64_t*); 00112 icp_stats.ACPs= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t); 00113 00114 cp_list= (bmrm_ll*) LIBBMRM_CALLOC(1, bmrm_ll); 00115 00116 prevW= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t); 00117 00118 wt= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t); 00119 00120 if (H==NULL || A==NULL || b==NULL || beta==NULL || subgrad==NULL || 00121 diag_H==NULL || I==NULL || icp_stats.ICPcounter==NULL || 00122 icp_stats.ICPs==NULL || icp_stats.ACPs==NULL || 00123 cp_list==NULL || prevW==NULL || wt==NULL) 00124 { 00125 ppbmrm.exitflag=-2; 00126 goto cleanup; 00127 } 00128 00129 map= (bool*) LIBBMRM_CALLOC(BufSize, bool); 00130 00131 if (map==NULL) 00132 { 00133 ppbmrm.exitflag=-2; 00134 goto cleanup; 00135 } 00136 00137 memset( (bool*) map, true, BufSize); 00138 00139 /* Temporary buffers for ICP removal */ 00140 icp_stats.H_buff= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, float64_t); 00141 00142 if (icp_stats.H_buff==NULL) 00143 { 00144 ppbmrm.exitflag=-2; 00145 goto cleanup; 00146 } 00147 00148 /* Temporary buffers */ 00149 beta_start= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t); 00150 00151 beta_good= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t); 00152 00153 b2= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t); 00154 00155 diag_H2= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t); 00156 00157 H2= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, float64_t); 00158 00159 I_start= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t); 00160 00161 I_good= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t); 00162 00163 I2= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t); 00164 00165 if (beta_start==NULL || beta_good==NULL || b2==NULL || diag_H2==NULL || 00166 I_start==NULL || I_good==NULL || I2==NULL || H2==NULL) 00167 { 00168 ppbmrm.exitflag=-2; 00169 goto cleanup; 00170 } 00171 00172 ppbmrm.hist_Fp.resize_vector(BufSize); 00173 ppbmrm.hist_Fd.resize_vector(BufSize); 00174 ppbmrm.hist_wdist.resize_vector(BufSize); 00175 00176 /* Iinitial solution */ 00177 R = machine->risk(subgrad, W); 00178 00179 ppbmrm.nCP=0; 00180 ppbmrm.nIter=0; 00181 ppbmrm.exitflag=0; 00182 00183 b[0]=-R; 00184 00185 /* Cutting plane auxiliary double linked list */ 00186 LIBBMRM_MEMCPY(A, subgrad, nDim*sizeof(float64_t)); 00187 map[0]=false; 00188 cp_list->address=&A[0]; 00189 cp_list->idx=0; 00190 cp_list->prev=NULL; 00191 cp_list->next=NULL; 00192 CPList_head=cp_list; 00193 CPList_tail=cp_list; 00194 00195 /* Compute initial value of Fp, Fd, assuming that W is zero vector */ 00196 sq_norm_Wdiff=0.0; 00197 00198 b[0] = SGVector<float64_t>::dot(subgrad, W, nDim); 00199 sq_norm_W = SGVector<float64_t>::dot(W, W, nDim); 00200 for (uint32_t j=0; j<nDim; ++j) 00201 { 00202 sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]); 00203 } 00204 00205 ppbmrm.Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff; 00206 ppbmrm.Fd=-LIBBMRM_PLUS_INF; 00207 lastFp=ppbmrm.Fp; 00208 wdist=CMath::sqrt(sq_norm_Wdiff); 00209 00210 K = (sq_norm_W == 0.0) ? 0.4 : 0.01*CMath::sqrt(sq_norm_W); 00211 00212 LIBBMRM_MEMCPY(prevW, W, nDim*sizeof(float64_t)); 00213 00214 tstop=ttime.cur_time_diff(false); 00215 00216 /* Keep history of Fp, Fd, wdist */ 00217 ppbmrm.hist_Fp[0]=ppbmrm.Fp; 00218 ppbmrm.hist_Fd[0]=ppbmrm.Fd; 00219 ppbmrm.hist_wdist[0]=wdist; 00220 00221 /* Verbose output */ 00222 00223 if (verbose) 00224 SG_SDEBUG("%4d: tim=%.3lf, Fp=%lf, Fd=%lf, R=%lf, K=%lf\n", 00225 ppbmrm.nIter, tstop-tstart, ppbmrm.Fp, ppbmrm.Fd, R, K); 00226 00227 if (verbose) 00228 helper = machine->get_helper(); 00229 00230 /* main loop */ 00231 00232 while (ppbmrm.exitflag==0) 00233 { 00234 tstart=ttime.cur_time_diff(false); 00235 ppbmrm.nIter++; 00236 00237 /* Update H */ 00238 00239 if (ppbmrm.nCP>0) 00240 { 00241 A_2=get_cutting_plane(CPList_tail); 00242 cp_ptr=CPList_head; 00243 00244 for (uint32_t i=0; i<ppbmrm.nCP; ++i) 00245 { 00246 A_1=get_cutting_plane(cp_ptr); 00247 cp_ptr=cp_ptr->next; 00248 rsum=SGVector<float64_t>::dot(A_1, A_2, nDim); 00249 00250 H[LIBBMRM_INDEX(ppbmrm.nCP, i, BufSize)] 00251 = H[LIBBMRM_INDEX(i, ppbmrm.nCP, BufSize)] 00252 = rsum; 00253 } 00254 } 00255 00256 A_2=get_cutting_plane(CPList_tail); 00257 rsum = SGVector<float64_t>::dot(A_2, A_2, nDim); 00258 00259 H[LIBBMRM_INDEX(ppbmrm.nCP, ppbmrm.nCP, BufSize)]=rsum; 00260 00261 diag_H[ppbmrm.nCP]=H[LIBBMRM_INDEX(ppbmrm.nCP, ppbmrm.nCP, BufSize)]; 00262 I[ppbmrm.nCP]=1; 00263 00264 ppbmrm.nCP++; 00265 beta[ppbmrm.nCP]=0.0; // [beta; 0] 00266 00267 /* tune alpha cycle */ 00268 /* ---------------------------------------------------------------------- */ 00269 00270 flag=true; 00271 isThereGoodSolution=false; 00272 LIBBMRM_MEMCPY(beta_start, beta, ppbmrm.nCP*sizeof(float64_t)); 00273 LIBBMRM_MEMCPY(I_start, I, ppbmrm.nCP*sizeof(uint32_t)); 00274 qp_cnt=0; 00275 alpha_good=alpha; 00276 00277 if (tuneAlpha) 00278 { 00279 alpha_start=alpha; alpha=0.0; 00280 beta[ppbmrm.nCP]=0.0; 00281 LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*sizeof(uint32_t)); 00282 I2[ppbmrm.nCP]=1; 00283 00284 /* add alpha-dependent terms to H, diag_h and b */ 00285 cp_ptr=CPList_head; 00286 00287 for (uint32_t i=0; i<ppbmrm.nCP; ++i) 00288 { 00289 A_1=get_cutting_plane(cp_ptr); 00290 cp_ptr=cp_ptr->next; 00291 00292 rsum = SGVector<float64_t>::dot(A_1, prevW, nDim); 00293 00294 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum; 00295 diag_H2[i]=diag_H[i]/(_lambda+2*alpha); 00296 00297 for (uint32_t j=0; j<ppbmrm.nCP; ++j) 00298 H2[LIBBMRM_INDEX(i, j, BufSize)]= 00299 H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha); 00300 } 00301 00302 /* solve QP with current alpha */ 00303 qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta, 00304 ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0); 00305 ppbmrm.qp_exitflag=qp_exitflag.exitflag; 00306 qp_cnt++; 00307 Fd_alpha0=-qp_exitflag.QP; 00308 00309 /* obtain w_t and check if norm(w_{t+1} -w_t) <= K */ 00310 for (uint32_t i=0; i<nDim; ++i) 00311 { 00312 rsum=0.0; 00313 cp_ptr=CPList_head; 00314 00315 for (uint32_t j=0; j<ppbmrm.nCP; ++j) 00316 { 00317 A_1=get_cutting_plane(cp_ptr); 00318 cp_ptr=cp_ptr->next; 00319 rsum+=A_1[i]*beta[j]; 00320 } 00321 00322 wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha); 00323 } 00324 00325 sq_norm_Wdiff=0.0; 00326 00327 for (uint32_t i=0; i<nDim; ++i) 00328 sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]); 00329 00330 if (CMath::sqrt(sq_norm_Wdiff) <= K) 00331 { 00332 flag=false; 00333 00334 if (alpha!=alpha_start) 00335 alphaChanged=true; 00336 } 00337 else 00338 { 00339 alpha=alpha_start; 00340 } 00341 00342 while(flag) 00343 { 00344 LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*sizeof(uint32_t)); 00345 LIBBMRM_MEMCPY(beta, beta_start, ppbmrm.nCP*sizeof(float64_t)); 00346 I2[ppbmrm.nCP]=1; 00347 beta[ppbmrm.nCP]=0.0; 00348 00349 /* add alpha-dependent terms to H, diag_h and b */ 00350 cp_ptr=CPList_head; 00351 00352 for (uint32_t i=0; i<ppbmrm.nCP; ++i) 00353 { 00354 A_1=get_cutting_plane(cp_ptr); 00355 cp_ptr=cp_ptr->next; 00356 00357 rsum = SGVector<float64_t>::dot(A_1, prevW, nDim); 00358 00359 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum; 00360 diag_H2[i]=diag_H[i]/(_lambda+2*alpha); 00361 00362 for (uint32_t j=0; j<ppbmrm.nCP; ++j) 00363 H2[LIBBMRM_INDEX(i, j, BufSize)]=H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha); 00364 } 00365 00366 /* solve QP with current alpha */ 00367 qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta, 00368 ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0); 00369 ppbmrm.qp_exitflag=qp_exitflag.exitflag; 00370 qp_cnt++; 00371 00372 /* obtain w_t and check if norm(w_{t+1}-w_t) <= K */ 00373 for (uint32_t i=0; i<nDim; ++i) 00374 { 00375 rsum=0.0; 00376 cp_ptr=CPList_head; 00377 00378 for (uint32_t j=0; j<ppbmrm.nCP; ++j) 00379 { 00380 A_1=get_cutting_plane(cp_ptr); 00381 cp_ptr=cp_ptr->next; 00382 rsum+=A_1[i]*beta[j]; 00383 } 00384 00385 wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha); 00386 } 00387 00388 sq_norm_Wdiff=0.0; 00389 for (uint32_t i=0; i<nDim; ++i) 00390 sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]); 00391 00392 if (CMath::sqrt(sq_norm_Wdiff) > K) 00393 { 00394 /* if there is a record of some good solution 00395 * (i.e. adjust alpha by division by 2) */ 00396 00397 if (isThereGoodSolution) 00398 { 00399 LIBBMRM_MEMCPY(beta, beta_good, ppbmrm.nCP*sizeof(float64_t)); 00400 LIBBMRM_MEMCPY(I2, I_good, ppbmrm.nCP*sizeof(uint32_t)); 00401 alpha=alpha_good; 00402 qp_exitflag=qp_exitflag_good; 00403 flag=false; 00404 } 00405 else 00406 { 00407 if (alpha == 0) 00408 { 00409 alpha=1.0; 00410 alphaChanged=true; 00411 } 00412 else 00413 { 00414 alpha*=2; 00415 alphaChanged=true; 00416 } 00417 } 00418 } 00419 else 00420 { 00421 if (alpha > 0) 00422 { 00423 /* keep good solution and try for alpha /= 2 if previous alpha was 1 */ 00424 LIBBMRM_MEMCPY(beta_good, beta, ppbmrm.nCP*sizeof(float64_t)); 00425 LIBBMRM_MEMCPY(I_good, I2, ppbmrm.nCP*sizeof(uint32_t)); 00426 alpha_good=alpha; 00427 qp_exitflag_good=qp_exitflag; 00428 isThereGoodSolution=true; 00429 00430 if (alpha!=1.0) 00431 { 00432 alpha/=2.0; 00433 alphaChanged=true; 00434 } 00435 else 00436 { 00437 alpha=0.0; 00438 alphaChanged=true; 00439 } 00440 } 00441 else 00442 { 00443 flag=false; 00444 } 00445 } 00446 } 00447 } 00448 else 00449 { 00450 alphaChanged=false; 00451 LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*sizeof(uint32_t)); 00452 LIBBMRM_MEMCPY(beta, beta_start, ppbmrm.nCP*sizeof(float64_t)); 00453 00454 /* add alpha-dependent terms to H, diag_h and b */ 00455 cp_ptr=CPList_head; 00456 00457 for (uint32_t i=0; i<ppbmrm.nCP; ++i) 00458 { 00459 A_1=get_cutting_plane(cp_ptr); 00460 cp_ptr=cp_ptr->next; 00461 00462 rsum = SGVector<float64_t>::dot(A_1, prevW, nDim); 00463 00464 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum; 00465 diag_H2[i]=diag_H[i]/(_lambda+2*alpha); 00466 00467 for (uint32_t j=0; j<ppbmrm.nCP; ++j) 00468 H2[LIBBMRM_INDEX(i, j, BufSize)]= 00469 H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha); 00470 } 00471 /* solve QP with current alpha */ 00472 qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta, 00473 ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0); 00474 ppbmrm.qp_exitflag=qp_exitflag.exitflag; 00475 qp_cnt++; 00476 } 00477 00478 /* ----------------------------------------------------------------------------------------------- */ 00479 00480 /* Update ICPcounter (add one to unused and reset used) + compute number of active CPs */ 00481 ppbmrm.nzA=0; 00482 00483 for (uint32_t aaa=0; aaa<ppbmrm.nCP; ++aaa) 00484 { 00485 if (beta[aaa]>epsilon) 00486 { 00487 ++ppbmrm.nzA; 00488 icp_stats.ICPcounter[aaa]=0; 00489 } 00490 else 00491 { 00492 icp_stats.ICPcounter[aaa]+=1; 00493 } 00494 } 00495 00496 /* W update */ 00497 for (uint32_t i=0; i<nDim; ++i) 00498 { 00499 rsum=0.0; 00500 cp_ptr=CPList_head; 00501 00502 for (uint32_t j=0; j<ppbmrm.nCP; ++j) 00503 { 00504 A_1=get_cutting_plane(cp_ptr); 00505 cp_ptr=cp_ptr->next; 00506 rsum+=A_1[i]*beta[j]; 00507 } 00508 00509 W[i]=(2*alpha*prevW[i]-rsum)/(_lambda+2*alpha); 00510 } 00511 00512 /* risk and subgradient computation */ 00513 R = machine->risk(subgrad, W); 00514 add_cutting_plane(&CPList_tail, map, A, 00515 find_free_idx(map, BufSize), subgrad, nDim); 00516 00517 sq_norm_W=SGVector<float64_t>::dot(W, W, nDim); 00518 sq_norm_prevW=SGVector<float64_t>::dot(prevW, prevW, nDim); 00519 b[ppbmrm.nCP]=SGVector<float64_t>::dot(subgrad, W, nDim) - R; 00520 00521 sq_norm_Wdiff=0.0; 00522 for (uint32_t j=0; j<nDim; ++j) 00523 { 00524 sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]); 00525 } 00526 00527 /* compute Fp and Fd */ 00528 ppbmrm.Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff; 00529 ppbmrm.Fd=-qp_exitflag.QP+((alpha*_lambda)/(_lambda + 2*alpha))*sq_norm_prevW; 00530 00531 /* gamma + tuneAlpha flag */ 00532 if (alphaChanged) 00533 { 00534 eps=1.0-(ppbmrm.Fd/ppbmrm.Fp); 00535 gamma=(lastFp*(1-eps)-Fd_alpha0)/(Tmax*(1-eps)); 00536 } 00537 00538 if ((lastFp-ppbmrm.Fp) <= gamma) 00539 { 00540 tuneAlpha=true; 00541 } 00542 else 00543 { 00544 tuneAlpha=false; 00545 } 00546 00547 /* Stopping conditions - set only with nonzero alpha */ 00548 if (alpha==0.0) 00549 { 00550 if (ppbmrm.Fp-ppbmrm.Fd<=TolRel*LIBBMRM_ABS(ppbmrm.Fp)) 00551 ppbmrm.exitflag=1; 00552 00553 if (ppbmrm.Fp-ppbmrm.Fd<=TolAbs) 00554 ppbmrm.exitflag=2; 00555 } 00556 00557 if (ppbmrm.nCP>=BufSize) 00558 ppbmrm.exitflag=-1; 00559 00560 tstop=ttime.cur_time_diff(false); 00561 00562 /* compute wdist (= || W_{t+1} - W_{t} || ) */ 00563 sq_norm_Wdiff=0.0; 00564 00565 for (uint32_t i=0; i<nDim; ++i) 00566 { 00567 sq_norm_Wdiff+=(W[i]-prevW[i])*(W[i]-prevW[i]); 00568 } 00569 00570 wdist=CMath::sqrt(sq_norm_Wdiff); 00571 00572 /* Keep history of Fp, Fd, wdist */ 00573 ppbmrm.hist_Fp[ppbmrm.nIter]=ppbmrm.Fp; 00574 ppbmrm.hist_Fd[ppbmrm.nIter]=ppbmrm.Fd; 00575 ppbmrm.hist_wdist[ppbmrm.nIter]=wdist; 00576 00577 /* Verbose output */ 00578 if (verbose) 00579 SG_SDEBUG("%4d: tim=%.3lf, Fp=%lf, Fd=%lf, (Fp-Fd)=%lf, (Fp-Fd)/Fp=%lf, R=%lf, nCP=%d, nzA=%d, wdist=%lf, alpha=%lf, qp_cnt=%d, gamma=%lf, tuneAlpha=%d\n", 00580 ppbmrm.nIter, tstop-tstart, ppbmrm.Fp, ppbmrm.Fd, ppbmrm.Fp-ppbmrm.Fd, 00581 (ppbmrm.Fp-ppbmrm.Fd)/ppbmrm.Fp, R, ppbmrm.nCP, ppbmrm.nzA, wdist, alpha, 00582 qp_cnt, gamma, tuneAlpha); 00583 00584 /* Check size of Buffer */ 00585 if (ppbmrm.nCP>=BufSize) 00586 { 00587 ppbmrm.exitflag=-2; 00588 SG_SERROR("Buffer exceeded.\n") 00589 } 00590 00591 /* keep w_t + Fp */ 00592 LIBBMRM_MEMCPY(prevW, W, nDim*sizeof(float64_t)); 00593 lastFp=ppbmrm.Fp; 00594 00595 /* Inactive Cutting Planes (ICP) removal */ 00596 if (cleanICP) 00597 { 00598 clean_icp(&icp_stats, ppbmrm, &CPList_head, &CPList_tail, H, diag_H, beta, map, cleanAfter, b, I); 00599 } 00600 00601 /* Debug: compute objective and training error */ 00602 if (verbose) 00603 { 00604 SGVector<float64_t> w_debug(W, nDim, false); 00605 float64_t primal = CSOSVMHelper::primal_objective(w_debug, model, _lambda); 00606 float64_t train_error = CSOSVMHelper::average_loss(w_debug, model); 00607 helper->add_debug_info(primal, ppbmrm.nIter, train_error); 00608 } 00609 } /* end of main loop */ 00610 00611 if (verbose) 00612 { 00613 helper->terminate(); 00614 SG_UNREF(helper); 00615 } 00616 00617 ppbmrm.hist_Fp.resize_vector(ppbmrm.nIter); 00618 ppbmrm.hist_Fd.resize_vector(ppbmrm.nIter); 00619 ppbmrm.hist_wdist.resize_vector(ppbmrm.nIter); 00620 00621 cp_ptr=CPList_head; 00622 00623 while(cp_ptr!=NULL) 00624 { 00625 cp_ptr2=cp_ptr; 00626 cp_ptr=cp_ptr->next; 00627 LIBBMRM_FREE(cp_ptr2); 00628 cp_ptr2=NULL; 00629 } 00630 00631 cp_list=NULL; 00632 00633 cleanup: 00634 00635 LIBBMRM_FREE(H); 00636 LIBBMRM_FREE(b); 00637 LIBBMRM_FREE(beta); 00638 LIBBMRM_FREE(A); 00639 LIBBMRM_FREE(subgrad); 00640 LIBBMRM_FREE(diag_H); 00641 LIBBMRM_FREE(I); 00642 LIBBMRM_FREE(icp_stats.ICPcounter); 00643 LIBBMRM_FREE(icp_stats.ICPs); 00644 LIBBMRM_FREE(icp_stats.ACPs); 00645 LIBBMRM_FREE(icp_stats.H_buff); 00646 LIBBMRM_FREE(map); 00647 LIBBMRM_FREE(prevW); 00648 LIBBMRM_FREE(wt); 00649 LIBBMRM_FREE(beta_start); 00650 LIBBMRM_FREE(beta_good); 00651 LIBBMRM_FREE(I_start); 00652 LIBBMRM_FREE(I_good); 00653 LIBBMRM_FREE(I2); 00654 LIBBMRM_FREE(b2); 00655 LIBBMRM_FREE(diag_H2); 00656 LIBBMRM_FREE(H2); 00657 00658 if (cp_list) 00659 LIBBMRM_FREE(cp_list); 00660 00661 SG_UNREF(model); 00662 00663 return(ppbmrm); 00664 } 00665 }