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