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