SHOGUN
v2.0.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 * libppbm.h: Implementation of the Proximal Point BM 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 (p3bm) 00012 *--------------------------------------------------------------------- */ 00013 00014 #include <stdlib.h> 00015 #include <string.h> 00016 #include <math.h> 00017 00018 #include <shogun/structure/libppbm.h> 00019 #include <shogun/lib/external/libqp.h> 00020 #include <shogun/lib/Time.h> 00021 00022 namespace shogun 00023 { 00024 static const uint32_t QPSolverMaxIter=0xFFFFFFFF; 00025 static const float64_t epsilon=0.0; 00026 00027 static float64_t *H, *H2; 00028 static uint32_t BufSize; 00029 00030 /*---------------------------------------------------------------------- 00031 Returns pointer at i-th column of Hessian matrix. 00032 ----------------------------------------------------------------------*/ 00033 static const float64_t *get_col( uint32_t i) 00034 { 00035 return( &H2[ BufSize*i ] ); 00036 } 00037 00038 bmrm_return_value_T svm_p3bm_solver( 00039 CStructuredModel* model, 00040 float64_t* W, 00041 float64_t TolRel, 00042 float64_t TolAbs, 00043 float64_t _lambda, 00044 uint32_t _BufSize, 00045 bool cleanICP, 00046 uint32_t cleanAfter, 00047 float64_t K, 00048 uint32_t Tmax, 00049 uint32_t cp_models, 00050 bool verbose) 00051 { 00052 bmrm_return_value_T p3bmrm; 00053 libqp_state_T qp_exitflag={0, 0, 0, 0}, qp_exitflag_good={0, 0, 0, 0}; 00054 float64_t *b, *b2, *beta, *beta_good, *beta_start, *diag_H, *diag_H2; 00055 float64_t R, *Rt, **subgrad_t, *A, QPSolverTolRel, *C=NULL; 00056 float64_t *prevW, *wt, alpha, alpha_start, alpha_good=0.0, Fd_alpha0=0.0; 00057 float64_t lastFp, wdist, gamma=0.0; 00058 floatmax_t rsum, sq_norm_W, sq_norm_Wdiff, sq_norm_prevW, eps; 00059 uint32_t *I, *I2, *I_start, *I_good, *ICPcounter, *ACPs, cntICP=0, cntACP=0; 00060 uint8_t *S=NULL; 00061 uint32_t nCP_new=0, qp_cnt=0; 00062 bmrm_ll *CPList_head, *CPList_tail, *cp_ptr, *cp_ptr2, *cp_list=NULL; 00063 float64_t *A_1=NULL, *H_buff, **ICPs; 00064 bool *map=NULL, tuneAlpha=true, flag=true; 00065 bool alphaChanged=false, isThereGoodSolution=false; 00066 TMultipleCPinfo **info=NULL; 00067 uint32_t nDim=model->get_dim(); 00068 uint32_t to=0, N=0, cp_i=0; 00069 00070 CTime ttime; 00071 float64_t tstart, tstop; 00072 00073 00074 tstart=ttime.cur_time_diff(false); 00075 00076 BufSize=_BufSize*cp_models; 00077 QPSolverTolRel=1e-9; 00078 00079 H=NULL; 00080 b=NULL; 00081 beta=NULL; 00082 A=NULL; 00083 subgrad_t=NULL; 00084 diag_H=NULL; 00085 I=NULL; 00086 ICPcounter=NULL; 00087 ICPs=NULL; 00088 ACPs=NULL; 00089 prevW=NULL; 00090 wt=NULL; 00091 H_buff=NULL; 00092 diag_H2=NULL; 00093 b2=NULL; 00094 I2=NULL; 00095 H2=NULL; 00096 I_good=NULL; 00097 I_start=NULL; 00098 beta_start=NULL; 00099 beta_good=NULL; 00100 00101 alpha=0.0; 00102 00103 H= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, sizeof(float64_t)); 00104 00105 A= (float64_t*) LIBBMRM_CALLOC(nDim*BufSize, sizeof(float64_t)); 00106 00107 b= (float64_t*) LIBBMRM_CALLOC(BufSize, sizeof(float64_t)); 00108 00109 beta= (float64_t*) LIBBMRM_CALLOC(BufSize, sizeof(float64_t)); 00110 00111 subgrad_t= (float64_t**) LIBBMRM_CALLOC(cp_models, sizeof(float64_t*)); 00112 00113 Rt= (float64_t*) LIBBMRM_CALLOC(cp_models, sizeof(float64_t)); 00114 00115 diag_H= (float64_t*) LIBBMRM_CALLOC(BufSize, sizeof(float64_t)); 00116 00117 I= (uint32_t*) LIBBMRM_CALLOC(BufSize, sizeof(uint32_t)); 00118 00119 ICPcounter= (uint32_t*) LIBBMRM_CALLOC(BufSize, sizeof(uint32_t)); 00120 00121 ICPs= (float64_t**) LIBBMRM_CALLOC(BufSize, sizeof(float64_t*)); 00122 00123 ACPs= (uint32_t*) LIBBMRM_CALLOC(BufSize, sizeof(uint32_t)); 00124 00125 cp_list= (bmrm_ll*) LIBBMRM_CALLOC(1, sizeof(bmrm_ll)); 00126 00127 prevW= (float64_t*) LIBBMRM_CALLOC(nDim, sizeof(float64_t)); 00128 00129 wt= (float64_t*) LIBBMRM_CALLOC(nDim, sizeof(float64_t)); 00130 00131 C= (float64_t*) LIBBMRM_CALLOC(cp_models, sizeof(float64_t)); 00132 00133 S= (uint8_t*) LIBBMRM_CALLOC(cp_models, sizeof(uint8_t)); 00134 00135 info= (TMultipleCPinfo**) LIBBMRM_CALLOC(cp_models, sizeof(TMultipleCPinfo*)); 00136 00137 int32_t num_feats = model->get_features()->get_num_vectors(); 00138 00139 if (H==NULL || A==NULL || b==NULL || beta==NULL || subgrad_t==NULL || 00140 diag_H==NULL || I==NULL || ICPcounter==NULL || ICPs==NULL || ACPs==NULL || 00141 cp_list==NULL || prevW==NULL || wt==NULL || Rt==NULL || C==NULL || 00142 S==NULL || info==NULL) 00143 { 00144 p3bmrm.exitflag=-2; 00145 goto cleanup; 00146 } 00147 00148 /* multiple cutting plane model init */ 00149 00150 to=0; 00151 N= (uint32_t) round( (float64_t) ((float64_t)num_feats / (float64_t) cp_models)); 00152 00153 for (uint32_t p=0; p<cp_models; ++p) 00154 { 00155 S[p]=1; 00156 C[p]=1.0; 00157 info[p]=(TMultipleCPinfo*)LIBBMRM_CALLOC(1, sizeof(TMultipleCPinfo)); 00158 subgrad_t[p]=(float64_t*)LIBBMRM_CALLOC(nDim, sizeof(float64_t)); 00159 00160 if (subgrad_t[p]==NULL || info[p]==NULL) 00161 { 00162 p3bmrm.exitflag=-2; 00163 goto cleanup; 00164 } 00165 00166 info[p]->_from=to; 00167 to=((p+1)*N > (uint32_t)num_feats) ? (uint32_t)num_feats : (p+1)*N; 00168 info[p]->N=to-info[p]->_from; 00169 } 00170 00171 map= (bool*) LIBBMRM_CALLOC(BufSize, sizeof(bool)); 00172 00173 if (map==NULL) 00174 { 00175 p3bmrm.exitflag=-2; 00176 goto cleanup; 00177 } 00178 00179 memset( (bool*) map, true, BufSize); 00180 00181 /* Temporary buffers for ICP removal */ 00182 H_buff= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, sizeof(float64_t)); 00183 00184 if (H_buff==NULL) 00185 { 00186 p3bmrm.exitflag=-2; 00187 goto cleanup; 00188 } 00189 00190 /* Temporary buffers */ 00191 beta_start= (float64_t*) LIBBMRM_CALLOC(BufSize, sizeof(float64_t)); 00192 00193 beta_good= (float64_t*) LIBBMRM_CALLOC(BufSize, sizeof(float64_t)); 00194 00195 b2= (float64_t*) LIBBMRM_CALLOC(BufSize, sizeof(float64_t)); 00196 00197 diag_H2= (float64_t*) LIBBMRM_CALLOC(BufSize, sizeof(float64_t)); 00198 00199 H2= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, sizeof(float64_t)); 00200 00201 I_start= (uint32_t*) LIBBMRM_CALLOC(BufSize, sizeof(uint32_t)); 00202 00203 I_good= (uint32_t*) LIBBMRM_CALLOC(BufSize, sizeof(uint32_t)); 00204 00205 I2= (uint32_t*) LIBBMRM_CALLOC(BufSize, sizeof(uint32_t)); 00206 00207 if (beta_start==NULL || beta_good==NULL || b2==NULL || diag_H2==NULL || 00208 I_start==NULL || I_good==NULL || I2==NULL || H2==NULL) 00209 { 00210 p3bmrm.exitflag=-2; 00211 goto cleanup; 00212 } 00213 00214 p3bmrm.hist_Fp.resize_vector(BufSize); 00215 p3bmrm.hist_Fd.resize_vector(BufSize); 00216 p3bmrm.hist_wdist.resize_vector(BufSize); 00217 00218 /* Iinitial solution */ 00219 Rt[0] = model->risk(subgrad_t[0], W, info[0]); 00220 00221 p3bmrm.nCP=0; 00222 p3bmrm.nIter=0; 00223 p3bmrm.exitflag=0; 00224 00225 b[0]=-Rt[0]; 00226 00227 /* Cutting plane auxiliary double linked list */ 00228 LIBBMRM_MEMCPY(A, subgrad_t[0], nDim*sizeof(float64_t)); 00229 map[0]=false; 00230 cp_list->address=&A[0]; 00231 cp_list->idx=0; 00232 cp_list->prev=NULL; 00233 cp_list->next=NULL; 00234 CPList_head=cp_list; 00235 CPList_tail=cp_list; 00236 00237 for (uint32_t p=1; p<cp_models; ++p) 00238 { 00239 Rt[p] = model->risk(subgrad_t[p], W, info[p]); 00240 b[p]=-Rt[p]; 00241 add_cutting_plane(&CPList_tail, map, A, find_free_idx(map, BufSize), subgrad_t[p], nDim); 00242 } 00243 00244 /* Compute initial value of Fp, Fd, assuming that W is zero vector */ 00245 R=0.0; 00246 00247 for (uint32_t p=0; p<cp_models; ++p) 00248 R+=Rt[p]; 00249 00250 sq_norm_W=0.0; 00251 sq_norm_Wdiff=0.0; 00252 00253 for (uint32_t j=0; j<nDim; ++j) 00254 { 00255 for (uint32_t p=0; p<cp_models; ++p) 00256 b[p]+=subgrad_t[p][j]*W[j]; 00257 00258 sq_norm_W+=W[j]*W[j]; 00259 sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]); 00260 } 00261 00262 wdist=::sqrt(sq_norm_Wdiff); 00263 00264 p3bmrm.Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff; 00265 p3bmrm.Fd=-LIBBMRM_PLUS_INF; 00266 lastFp=p3bmrm.Fp; 00267 00268 /* if there is initial W, then set K to be 0.01 times its norm */ 00269 K = (sq_norm_W == 0.0) ? 0.4 : 0.01*::sqrt(sq_norm_W); 00270 00271 LIBBMRM_MEMCPY(prevW, W, nDim*sizeof(float64_t)); 00272 00273 tstop=ttime.cur_time_diff(false); 00274 00275 /* Keep history of Fp, Fd, and wdist */ 00276 p3bmrm.hist_Fp[0]=p3bmrm.Fp; 00277 p3bmrm.hist_Fd[0]=p3bmrm.Fd; 00278 p3bmrm.hist_wdist[0]=wdist; 00279 00280 /* Verbose output */ 00281 if (verbose) 00282 SG_SPRINT("%4d: tim=%.3lf, Fp=%lf, Fd=%lf, R=%lf, K=%lf, CPmodels=%d\n", 00283 p3bmrm.nIter, tstop-tstart, p3bmrm.Fp, p3bmrm.Fd, R, K, cp_models); 00284 00285 /* main loop */ 00286 while (p3bmrm.exitflag==0) 00287 { 00288 tstart=ttime.cur_time_diff(false); 00289 p3bmrm.nIter++; 00290 00291 /* Update H */ 00292 if (p3bmrm.nIter==1) 00293 { 00294 cp_ptr=CPList_head; 00295 00296 for (cp_i=0; cp_i<cp_models; ++cp_i) /* for all cutting planes */ 00297 { 00298 A_1=get_cutting_plane(cp_ptr); 00299 00300 for (uint32_t p=0; p<cp_models; ++p) 00301 { 00302 rsum=0.0; 00303 00304 for (uint32_t j=0; j<nDim; ++j) 00305 { 00306 rsum+=subgrad_t[p][j]*A_1[j]; 00307 } 00308 00309 H[LIBBMRM_INDEX(p, cp_i, BufSize)]=rsum; 00310 } 00311 00312 cp_ptr=cp_ptr->next; 00313 } 00314 } 00315 else 00316 { 00317 cp_ptr=CPList_head; 00318 00319 for (cp_i=0; cp_i<p3bmrm.nCP+cp_models; ++cp_i) /* for all cutting planes */ 00320 { 00321 A_1=get_cutting_plane(cp_ptr); 00322 00323 for (uint32_t p=0; p<cp_models; ++p) 00324 { 00325 rsum=0.0; 00326 00327 for (uint32_t j=0; j<nDim; ++j) 00328 rsum+=subgrad_t[p][j]*A_1[j]; 00329 00330 H[LIBBMRM_INDEX(p3bmrm.nCP+p, cp_i, BufSize)]=rsum; 00331 } 00332 00333 cp_ptr=cp_ptr->next; 00334 } 00335 00336 for (uint32_t i=0; i<p3bmrm.nCP; ++i) 00337 for (uint32_t j=0; j<cp_models; ++j) 00338 H[LIBBMRM_INDEX(i, p3bmrm.nCP+j, BufSize)]= 00339 H[LIBBMRM_INDEX(p3bmrm.nCP+j, i, BufSize)]; 00340 } 00341 00342 for (uint32_t p=0; p<cp_models; ++p) 00343 diag_H[p3bmrm.nCP+p]=H[LIBBMRM_INDEX(p3bmrm.nCP+p, p3bmrm.nCP+p, BufSize)]; 00344 00345 p3bmrm.nCP+=cp_models; 00346 00347 /* tune alpha cycle */ 00348 /* ------------------------------------------------------------------------ */ 00349 flag=true; 00350 isThereGoodSolution=false; 00351 00352 for (uint32_t p=0; p<cp_models; ++p) 00353 { 00354 I[p3bmrm.nCP-cp_models+p]=p+1; 00355 beta[p3bmrm.nCP-cp_models+p]=0.0; 00356 } 00357 00358 LIBBMRM_MEMCPY(beta_start, beta, p3bmrm.nCP*sizeof(float64_t)); 00359 LIBBMRM_MEMCPY(I_start, I, p3bmrm.nCP*sizeof(uint32_t)); 00360 qp_cnt=0; 00361 00362 if (tuneAlpha) 00363 { 00364 alpha_start=alpha; alpha=0.0; 00365 LIBBMRM_MEMCPY(I2, I_start, p3bmrm.nCP*sizeof(uint32_t)); 00366 00367 /* add alpha-dependent terms to H, diag_h and b */ 00368 cp_ptr=CPList_head; 00369 00370 for (uint32_t i=0; i<p3bmrm.nCP; ++i) 00371 { 00372 rsum=0.0; 00373 A_1=get_cutting_plane(cp_ptr); 00374 cp_ptr=cp_ptr->next; 00375 00376 for (uint32_t j=0; j<nDim; ++j) 00377 rsum+=A_1[j]*prevW[j]; 00378 00379 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum; 00380 diag_H2[i]=diag_H[i]/(_lambda+2*alpha); 00381 00382 for (uint32_t j=0; j<p3bmrm.nCP; ++j) 00383 H2[LIBBMRM_INDEX(i, j, BufSize)]= 00384 H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha); 00385 00386 } 00387 00388 /* solve QP with current alpha */ 00389 qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, C, I2, S, beta, 00390 p3bmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0); 00391 p3bmrm.qp_exitflag=qp_exitflag.exitflag; 00392 qp_cnt++; 00393 Fd_alpha0=-qp_exitflag.QP; 00394 00395 /* obtain w_t and check if norm(w_{t+1} -w_t) <= K */ 00396 for (uint32_t i=0; i<nDim; ++i) 00397 { 00398 rsum=0.0; 00399 cp_ptr=CPList_head; 00400 00401 for (uint32_t j=0; j<p3bmrm.nCP; ++j) 00402 { 00403 A_1=get_cutting_plane(cp_ptr); 00404 cp_ptr=cp_ptr->next; 00405 rsum+=A_1[i]*beta[j]; 00406 } 00407 00408 wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha); 00409 } 00410 00411 sq_norm_Wdiff=0.0; 00412 00413 for (uint32_t i=0; i<nDim; ++i) 00414 sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]); 00415 00416 if (::sqrt(sq_norm_Wdiff) <= K) 00417 { 00418 flag=false; 00419 00420 if (alpha!=alpha_start) 00421 alphaChanged=true; 00422 } 00423 else 00424 { 00425 alpha=alpha_start; 00426 } 00427 00428 while(flag) 00429 { 00430 LIBBMRM_MEMCPY(I2, I_start, p3bmrm.nCP*sizeof(uint32_t)); 00431 LIBBMRM_MEMCPY(beta, beta_start, p3bmrm.nCP*sizeof(float64_t)); 00432 00433 /* add alpha-dependent terms to H, diag_h and b */ 00434 cp_ptr=CPList_head; 00435 00436 for (uint32_t i=0; i<p3bmrm.nCP; ++i) 00437 { 00438 rsum=0.0; 00439 A_1=get_cutting_plane(cp_ptr); 00440 cp_ptr=cp_ptr->next; 00441 00442 for (uint32_t j=0; j<nDim; ++j) 00443 rsum+=A_1[j]*prevW[j]; 00444 00445 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum; 00446 diag_H2[i]=diag_H[i]/(_lambda+2*alpha); 00447 00448 for (uint32_t j=0; j<p3bmrm.nCP; ++j) 00449 H2[LIBBMRM_INDEX(i, j, BufSize)]=H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha); 00450 } 00451 00452 /* solve QP with current alpha */ 00453 qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, C, I2, S, beta, 00454 p3bmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0); 00455 p3bmrm.qp_exitflag=qp_exitflag.exitflag; 00456 qp_cnt++; 00457 00458 /* obtain w_t and check if norm(w_{t+1}-w_t) <= K */ 00459 for (uint32_t i=0; i<nDim; ++i) 00460 { 00461 rsum=0.0; 00462 cp_ptr=CPList_head; 00463 00464 for (uint32_t j=0; j<p3bmrm.nCP; ++j) 00465 { 00466 A_1=get_cutting_plane(cp_ptr); 00467 cp_ptr=cp_ptr->next; 00468 rsum+=A_1[i]*beta[j]; 00469 } 00470 00471 wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha); 00472 } 00473 00474 sq_norm_Wdiff=0.0; 00475 00476 for (uint32_t i=0; i<nDim; ++i) 00477 sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]); 00478 00479 if (::sqrt(sq_norm_Wdiff) > K) 00480 { 00481 /* if there is a record of some good solution (i.e. adjust alpha by division by 2) */ 00482 00483 if (isThereGoodSolution) 00484 { 00485 LIBBMRM_MEMCPY(beta, beta_good, p3bmrm.nCP*sizeof(float64_t)); 00486 LIBBMRM_MEMCPY(I2, I_good, p3bmrm.nCP*sizeof(uint32_t)); 00487 alpha=alpha_good; 00488 qp_exitflag=qp_exitflag_good; 00489 flag=false; 00490 } 00491 else 00492 { 00493 if (alpha == 0) 00494 { 00495 alpha=1.0; 00496 alphaChanged=true; 00497 } 00498 else 00499 { 00500 alpha*=2; 00501 alphaChanged=true; 00502 } 00503 } 00504 } 00505 else 00506 { 00507 if (alpha > 0) 00508 { 00509 /* keep good solution and try for alpha /= 2 if previous alpha was 1 */ 00510 LIBBMRM_MEMCPY(beta_good, beta, p3bmrm.nCP*sizeof(float64_t)); 00511 LIBBMRM_MEMCPY(I_good, I2, p3bmrm.nCP*sizeof(uint32_t)); 00512 alpha_good=alpha; 00513 qp_exitflag_good=qp_exitflag; 00514 isThereGoodSolution=true; 00515 00516 if (alpha!=1.0) 00517 { 00518 alpha/=2.0; 00519 alphaChanged=true; 00520 } 00521 else 00522 { 00523 alpha=0.0; 00524 alphaChanged=true; 00525 } 00526 } 00527 else 00528 { 00529 flag=false; 00530 } 00531 } 00532 } 00533 } 00534 else 00535 { 00536 alphaChanged=false; 00537 LIBBMRM_MEMCPY(I2, I_start, p3bmrm.nCP*sizeof(uint32_t)); 00538 LIBBMRM_MEMCPY(beta, beta_start, p3bmrm.nCP*sizeof(float64_t)); 00539 00540 /* add alpha-dependent terms to H, diag_h and b */ 00541 cp_ptr=CPList_head; 00542 00543 for (uint32_t i=0; i<p3bmrm.nCP; ++i) 00544 { 00545 rsum=0.0; 00546 A_1=get_cutting_plane(cp_ptr); 00547 cp_ptr=cp_ptr->next; 00548 00549 for (uint32_t j=0; j<nDim; ++j) 00550 rsum+=A_1[j]*prevW[j]; 00551 00552 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum; 00553 diag_H2[i]=diag_H[i]/(_lambda+2*alpha); 00554 00555 for (uint32_t j=0; j<p3bmrm.nCP; ++j) 00556 H2[LIBBMRM_INDEX(i, j, BufSize)]=H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha); 00557 } 00558 00559 /* solve QP with current alpha */ 00560 qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, C, I2, S, beta, 00561 p3bmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0); 00562 p3bmrm.qp_exitflag=qp_exitflag.exitflag; 00563 qp_cnt++; 00564 } 00565 /* ----------------------------------------------------------------------------------------------- */ 00566 00567 /* Update ICPcounter (add one to unused and reset used) + compute number of active CPs */ 00568 p3bmrm.nzA=0; 00569 00570 for (uint32_t aaa=0; aaa<p3bmrm.nCP; ++aaa) 00571 { 00572 if (beta[aaa]>epsilon) 00573 { 00574 ++p3bmrm.nzA; 00575 ICPcounter[aaa]=0; 00576 } 00577 else 00578 { 00579 ICPcounter[aaa]+=1; 00580 } 00581 } 00582 00583 /* W update */ 00584 for (uint32_t i=0; i<nDim; ++i) 00585 { 00586 rsum=0.0; 00587 cp_ptr=CPList_head; 00588 00589 for (uint32_t j=0; j<p3bmrm.nCP; ++j) 00590 { 00591 A_1=get_cutting_plane(cp_ptr); 00592 cp_ptr=cp_ptr->next; 00593 rsum+=A_1[i]*beta[j]; 00594 } 00595 00596 W[i]=(2*alpha*prevW[i]-rsum)/(_lambda+2*alpha); 00597 } 00598 00599 /* risk and subgradient computation */ 00600 R=0.0; 00601 00602 for (uint32_t p=0; p<cp_models; ++p) 00603 { 00604 Rt[p] = model->risk(subgrad_t[p], W, info[p]); 00605 b[p3bmrm.nCP+p]=-Rt[p]; 00606 add_cutting_plane(&CPList_tail, map, A, find_free_idx(map, BufSize), subgrad_t[p], nDim); 00607 R+=Rt[p]; 00608 } 00609 00610 sq_norm_W=0.0; 00611 sq_norm_Wdiff=0.0; 00612 sq_norm_prevW=0.0; 00613 00614 for (uint32_t j=0; j<nDim; ++j) 00615 { 00616 for (uint32_t p=0; p<cp_models; ++p) 00617 b[p3bmrm.nCP+p]+=subgrad_t[p][j]*W[j]; 00618 00619 sq_norm_W+=W[j]*W[j]; 00620 sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]); 00621 sq_norm_prevW+=prevW[j]*prevW[j]; 00622 } 00623 00624 /* compute Fp and Fd */ 00625 p3bmrm.Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff; 00626 p3bmrm.Fd=-qp_exitflag.QP + ((alpha*_lambda)/(_lambda + 2*alpha))*sq_norm_prevW; 00627 00628 /* gamma + tuneAlpha flag */ 00629 if (alphaChanged) 00630 { 00631 eps=1.0-(p3bmrm.Fd/p3bmrm.Fp); 00632 gamma=(lastFp*(1-eps)-Fd_alpha0)/(Tmax*(1-eps)); 00633 } 00634 00635 if ((lastFp-p3bmrm.Fp) <= gamma) 00636 { 00637 tuneAlpha=true; 00638 } 00639 else 00640 { 00641 tuneAlpha=false; 00642 } 00643 00644 /* Stopping conditions - set only with nonzero alpha */ 00645 if (alpha==0.0) 00646 { 00647 if (p3bmrm.Fp-p3bmrm.Fd<=TolRel*LIBBMRM_ABS(p3bmrm.Fp)) 00648 p3bmrm.exitflag=1; 00649 00650 if (p3bmrm.Fp-p3bmrm.Fd<=TolAbs) 00651 p3bmrm.exitflag=2; 00652 } 00653 00654 if (p3bmrm.nCP>=BufSize) 00655 p3bmrm.exitflag=-1; 00656 00657 tstop=ttime.cur_time_diff(false); 00658 00659 /* compute wdist (= || W_{t+1} - W_{t} || ) */ 00660 sq_norm_Wdiff=0.0; 00661 00662 for (uint32_t i=0; i<nDim; ++i) 00663 { 00664 sq_norm_Wdiff+=(W[i]-prevW[i])*(W[i]-prevW[i]); 00665 } 00666 00667 wdist=::sqrt(sq_norm_Wdiff); 00668 00669 /* Keep history of Fp, Fd and wdist */ 00670 p3bmrm.hist_Fp[p3bmrm.nIter]=p3bmrm.Fp; 00671 p3bmrm.hist_Fd[p3bmrm.nIter]=p3bmrm.Fd; 00672 p3bmrm.hist_wdist[p3bmrm.nIter]=wdist; 00673 00674 /* Verbose output */ 00675 if (verbose) 00676 SG_SPRINT("%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", 00677 p3bmrm.nIter, tstop-tstart, p3bmrm.Fp, p3bmrm.Fd, p3bmrm.Fp-p3bmrm.Fd, 00678 (p3bmrm.Fp-p3bmrm.Fd)/p3bmrm.Fp, R, p3bmrm.nCP, p3bmrm.nzA, wdist, alpha, 00679 qp_cnt, gamma, tuneAlpha); 00680 00681 /* Check size of Buffer */ 00682 if (p3bmrm.nCP>=BufSize) 00683 { 00684 p3bmrm.exitflag=-2; 00685 SG_SERROR("Buffer exceeded.\n"); 00686 } 00687 00688 /* keep w_t + Fp */ 00689 LIBBMRM_MEMCPY(prevW, W, nDim*sizeof(float64_t)); 00690 lastFp=p3bmrm.Fp; 00691 00692 /* Inactive Cutting Planes (ICP) removal */ 00693 if (cleanICP) 00694 { 00695 /* find ICP */ 00696 cntICP=0; 00697 cntACP=0; 00698 cp_ptr=CPList_head; 00699 uint32_t tmp_idx=0; 00700 00701 while (cp_ptr != CPList_tail) 00702 { 00703 if (ICPcounter[tmp_idx++]>=cleanAfter) 00704 { 00705 ICPs[cntICP++]=cp_ptr->address; 00706 } 00707 else 00708 { 00709 ACPs[cntACP++]=tmp_idx-1; 00710 } 00711 00712 cp_ptr=cp_ptr->next; 00713 } 00714 00715 /* do ICP removal */ 00716 if (cntICP > 0) 00717 { 00718 nCP_new=p3bmrm.nCP-cntICP; 00719 00720 for (uint32_t i=0; i<cntICP; ++i) 00721 { 00722 tmp_idx=0; 00723 cp_ptr=CPList_head; 00724 00725 while(cp_ptr->address != ICPs[i]) 00726 { 00727 cp_ptr=cp_ptr->next; 00728 tmp_idx++; 00729 } 00730 00731 remove_cutting_plane(&CPList_head, &CPList_tail, map, ICPs[i]); 00732 00733 LIBBMRM_MEMMOVE(b+tmp_idx, b+tmp_idx+1, (p3bmrm.nCP+cp_models-tmp_idx)*sizeof(float64_t)); 00734 LIBBMRM_MEMMOVE(beta+tmp_idx, beta+tmp_idx+1, (p3bmrm.nCP-tmp_idx)*sizeof(float64_t)); 00735 LIBBMRM_MEMMOVE(diag_H+tmp_idx, diag_H+tmp_idx+1, (p3bmrm.nCP-tmp_idx)*sizeof(float64_t)); 00736 LIBBMRM_MEMMOVE(I+tmp_idx, I+tmp_idx+1, (p3bmrm.nCP-tmp_idx)*sizeof(uint32_t)); 00737 LIBBMRM_MEMMOVE(ICPcounter+tmp_idx, ICPcounter+tmp_idx+1, (p3bmrm.nCP-tmp_idx)*sizeof(uint32_t)); 00738 } 00739 00740 /* H */ 00741 for (uint32_t i=0; i < nCP_new; ++i) 00742 { 00743 for (uint32_t j=0; j < nCP_new; ++j) 00744 { 00745 H_buff[LIBBMRM_INDEX(i, j, BufSize)]=H[LIBBMRM_INDEX(ACPs[i], ACPs[j], BufSize)]; 00746 } 00747 } 00748 00749 for (uint32_t i=0; i<nCP_new; ++i) 00750 for (uint32_t j=0; j<nCP_new; ++j) 00751 H[LIBBMRM_INDEX(i, j, BufSize)]=H_buff[LIBBMRM_INDEX(i, j, BufSize)]; 00752 00753 p3bmrm.nCP=nCP_new; 00754 } 00755 } 00756 } /* end of main loop */ 00757 00758 p3bmrm.hist_Fp.resize_vector(p3bmrm.nIter); 00759 p3bmrm.hist_Fd.resize_vector(p3bmrm.nIter); 00760 p3bmrm.hist_wdist.resize_vector(p3bmrm.nIter); 00761 00762 cp_ptr=CPList_head; 00763 00764 while(cp_ptr!=NULL) 00765 { 00766 cp_ptr2=cp_ptr; 00767 cp_ptr=cp_ptr->next; 00768 LIBBMRM_FREE(cp_ptr2); 00769 cp_ptr2=NULL; 00770 } 00771 00772 cp_list=NULL; 00773 00774 cleanup: 00775 00776 LIBBMRM_FREE(H); 00777 LIBBMRM_FREE(b); 00778 LIBBMRM_FREE(beta); 00779 LIBBMRM_FREE(A); 00780 LIBBMRM_FREE(diag_H); 00781 LIBBMRM_FREE(I); 00782 LIBBMRM_FREE(ICPcounter); 00783 LIBBMRM_FREE(ICPs); 00784 LIBBMRM_FREE(ACPs); 00785 LIBBMRM_FREE(H_buff); 00786 LIBBMRM_FREE(map); 00787 LIBBMRM_FREE(prevW); 00788 LIBBMRM_FREE(wt); 00789 LIBBMRM_FREE(beta_start); 00790 LIBBMRM_FREE(beta_good); 00791 LIBBMRM_FREE(I_start); 00792 LIBBMRM_FREE(I_good); 00793 LIBBMRM_FREE(I2); 00794 LIBBMRM_FREE(b2); 00795 LIBBMRM_FREE(diag_H2); 00796 LIBBMRM_FREE(H2); 00797 LIBBMRM_FREE(C); 00798 LIBBMRM_FREE(S); 00799 LIBBMRM_FREE(Rt); 00800 00801 if (cp_list) 00802 LIBBMRM_FREE(cp_list); 00803 00804 for (uint32_t p=0; p<cp_models; ++p) 00805 { 00806 LIBBMRM_FREE(subgrad_t[p]); 00807 LIBBMRM_FREE(info[p]); 00808 } 00809 00810 LIBBMRM_FREE(subgrad_t); 00811 LIBBMRM_FREE(info); 00812 00813 return(p3bmrm); 00814 } 00815 }