16 #include <shogun/lib/external/libqp.h>
22 static const uint32_t QPSolverMaxIter=0xFFFFFFFF;
30 static const float64_t *get_col( uint32_t i)
32 return( &H2[ BufSize*i ] );
35 BmrmStatistics svm_ppbm_solver(
36 CDualLibQPBMSOSVM *machine,
48 BmrmStatistics ppbmrm;
49 libqp_state_T qp_exitflag={0, 0, 0, 0}, qp_exitflag_good={0, 0, 0, 0};
50 float64_t *b, *b2, *beta, *beta_good, *beta_start, *diag_H, *diag_H2;
51 float64_t R, *subgrad, *A, QPSolverTolRel, C=1.0;
52 float64_t *prevW, *wt, alpha, alpha_start, alpha_good=0.0, Fd_alpha0=0.0;
54 floatmax_t rsum, sq_norm_W, sq_norm_Wdiff, sq_norm_prevW, eps;
55 uint32_t *I, *I2, *I_start, *I_good;
57 CStructuredModel* model=machine->get_model();
58 uint32_t nDim=model->get_dim();
59 CSOSVMHelper* helper = NULL;
61 bmrm_ll *CPList_head, *CPList_tail, *cp_ptr, *cp_ptr2, *cp_list=NULL;
63 bool *map=NULL, tuneAlpha=
true, flag=
true, alphaChanged=
false, isThereGoodSolution=
false;
69 tstart=ttime.cur_time_diff(
false);
99 "overflow: %u * %u > %u -- biggest possible BufSize=%u or nDim=%u\n",
114 I= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
118 icp_stats.maxCPs = BufSize;
119 icp_stats.ICPcounter= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
121 icp_stats.ACPs= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
123 cp_list= (bmrm_ll*) LIBBMRM_CALLOC(1, bmrm_ll);
129 if (H==NULL || A==NULL || b==NULL || beta==NULL || subgrad==NULL ||
130 diag_H==NULL || I==NULL || icp_stats.ICPcounter==NULL ||
131 icp_stats.ICPs==NULL || icp_stats.ACPs==NULL ||
132 cp_list==NULL || prevW==NULL || wt==NULL)
138 map= (
bool*) LIBBMRM_CALLOC(BufSize,
bool);
146 memset( (
bool*) map,
true, BufSize);
151 if (icp_stats.H_buff==NULL)
168 I_start= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
170 I_good= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
172 I2= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
174 if (beta_start==NULL || beta_good==NULL || b2==NULL || diag_H2==NULL ||
175 I_start==NULL || I_good==NULL || I2==NULL || H2==NULL)
181 ppbmrm.hist_Fp.resize_vector(BufSize);
182 ppbmrm.hist_Fd.resize_vector(BufSize);
183 ppbmrm.hist_wdist.resize_vector(BufSize);
186 R = machine->risk(subgrad, W);
195 LIBBMRM_MEMCPY(A, subgrad, nDim*
sizeof(
float64_t));
197 cp_list->address=&A[0];
209 for (uint32_t j=0; j<nDim; ++j)
211 sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
214 ppbmrm.Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff;
215 ppbmrm.Fd=-LIBBMRM_PLUS_INF;
219 K = (sq_norm_W == 0.0) ? 0.4 : 0.01*
CMath::sqrt(sq_norm_W);
221 LIBBMRM_MEMCPY(prevW, W, nDim*
sizeof(
float64_t));
223 tstop=ttime.cur_time_diff(
false);
226 ppbmrm.hist_Fp[0]=ppbmrm.Fp;
227 ppbmrm.hist_Fd[0]=ppbmrm.Fd;
228 ppbmrm.hist_wdist[0]=wdist;
233 SG_SDEBUG(
"%4d: tim=%.3lf, Fp=%lf, Fd=%lf, R=%lf, K=%lf\n",
234 ppbmrm.nIter, tstop-tstart, ppbmrm.Fp, ppbmrm.Fd, R, K);
237 helper = machine->get_helper();
241 while (ppbmrm.exitflag==0)
243 tstart=ttime.cur_time_diff(
false);
250 A_2=get_cutting_plane(CPList_tail);
253 for (uint32_t i=0; i<ppbmrm.nCP; ++i)
255 A_1=get_cutting_plane(cp_ptr);
259 H[LIBBMRM_INDEX(ppbmrm.nCP, i, BufSize)]
260 = H[LIBBMRM_INDEX(i, ppbmrm.nCP, BufSize)]
265 A_2=get_cutting_plane(CPList_tail);
268 H[LIBBMRM_INDEX(ppbmrm.nCP, ppbmrm.nCP, BufSize)]=rsum;
270 diag_H[ppbmrm.nCP]=H[LIBBMRM_INDEX(ppbmrm.nCP, ppbmrm.nCP, BufSize)];
273 beta[ppbmrm.nCP]=0.0;
280 isThereGoodSolution=
false;
281 LIBBMRM_MEMCPY(beta_start, beta, ppbmrm.nCP*
sizeof(
float64_t));
282 LIBBMRM_MEMCPY(I_start, I, ppbmrm.nCP*
sizeof(uint32_t));
288 alpha_start=alpha; alpha=0.0;
289 beta[ppbmrm.nCP]=0.0;
290 LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*
sizeof(uint32_t));
296 for (uint32_t i=0; i<ppbmrm.nCP; ++i)
298 A_1=get_cutting_plane(cp_ptr);
303 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
304 diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
306 for (uint32_t j=0; j<ppbmrm.nCP; ++j)
307 H2[LIBBMRM_INDEX(i, j, BufSize)]=
308 H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha);
312 qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta,
313 ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
314 ppbmrm.qp_exitflag=qp_exitflag.exitflag;
316 Fd_alpha0=-qp_exitflag.QP;
319 for (uint32_t i=0; i<nDim; ++i)
324 for (uint32_t j=0; j<ppbmrm.nCP; ++j)
326 A_1=get_cutting_plane(cp_ptr);
328 rsum+=A_1[i]*beta[j];
331 wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha);
336 for (uint32_t i=0; i<nDim; ++i)
337 sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]);
343 if (alpha!=alpha_start)
353 LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*
sizeof(uint32_t));
354 LIBBMRM_MEMCPY(beta, beta_start, ppbmrm.nCP*
sizeof(
float64_t));
356 beta[ppbmrm.nCP]=0.0;
361 for (uint32_t i=0; i<ppbmrm.nCP; ++i)
363 A_1=get_cutting_plane(cp_ptr);
368 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
369 diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
371 for (uint32_t j=0; j<ppbmrm.nCP; ++j)
372 H2[LIBBMRM_INDEX(i, j, BufSize)]=H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha);
376 qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta,
377 ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
378 ppbmrm.qp_exitflag=qp_exitflag.exitflag;
382 for (uint32_t i=0; i<nDim; ++i)
387 for (uint32_t j=0; j<ppbmrm.nCP; ++j)
389 A_1=get_cutting_plane(cp_ptr);
391 rsum+=A_1[i]*beta[j];
394 wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha);
398 for (uint32_t i=0; i<nDim; ++i)
399 sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]);
406 if (isThereGoodSolution)
408 LIBBMRM_MEMCPY(beta, beta_good, ppbmrm.nCP*
sizeof(
float64_t));
409 LIBBMRM_MEMCPY(I2, I_good, ppbmrm.nCP*
sizeof(uint32_t));
411 qp_exitflag=qp_exitflag_good;
433 LIBBMRM_MEMCPY(beta_good, beta, ppbmrm.nCP*
sizeof(
float64_t));
434 LIBBMRM_MEMCPY(I_good, I2, ppbmrm.nCP*
sizeof(uint32_t));
436 qp_exitflag_good=qp_exitflag;
437 isThereGoodSolution=
true;
460 LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*
sizeof(uint32_t));
461 LIBBMRM_MEMCPY(beta, beta_start, ppbmrm.nCP*
sizeof(
float64_t));
466 for (uint32_t i=0; i<ppbmrm.nCP; ++i)
468 A_1=get_cutting_plane(cp_ptr);
473 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
474 diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
476 for (uint32_t j=0; j<ppbmrm.nCP; ++j)
477 H2[LIBBMRM_INDEX(i, j, BufSize)]=
478 H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha);
481 qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta,
482 ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
483 ppbmrm.qp_exitflag=qp_exitflag.exitflag;
492 for (uint32_t aaa=0; aaa<ppbmrm.nCP; ++aaa)
494 if (beta[aaa]>epsilon)
497 icp_stats.ICPcounter[aaa]=0;
501 icp_stats.ICPcounter[aaa]+=1;
506 for (uint32_t i=0; i<nDim; ++i)
511 for (uint32_t j=0; j<ppbmrm.nCP; ++j)
513 A_1=get_cutting_plane(cp_ptr);
515 rsum+=A_1[i]*beta[j];
518 W[i]=(2*alpha*prevW[i]-rsum)/(_lambda+2*alpha);
522 R = machine->risk(subgrad, W);
523 add_cutting_plane(&CPList_tail, map, A,
524 find_free_idx(map, BufSize), subgrad, nDim);
528 b[ppbmrm.nCP]=
CMath::dot(subgrad, W, nDim) - R;
531 for (uint32_t j=0; j<nDim; ++j)
533 sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
537 ppbmrm.Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff;
538 ppbmrm.Fd=-qp_exitflag.QP+((alpha*_lambda)/(_lambda + 2*alpha))*sq_norm_prevW;
543 eps=1.0-(ppbmrm.Fd/ppbmrm.Fp);
544 gamma=(lastFp*(1-eps)-Fd_alpha0)/(Tmax*(1-eps));
547 if ((lastFp-ppbmrm.Fp) <= gamma)
559 if (ppbmrm.Fp-ppbmrm.Fd<=TolRel*LIBBMRM_ABS(ppbmrm.Fp))
562 if (ppbmrm.Fp-ppbmrm.Fd<=TolAbs)
566 if (ppbmrm.nCP>=BufSize)
569 tstop=ttime.cur_time_diff(
false);
574 for (uint32_t i=0; i<nDim; ++i)
576 sq_norm_Wdiff+=(W[i]-prevW[i])*(W[i]-prevW[i]);
582 ppbmrm.hist_Fp[ppbmrm.nIter]=ppbmrm.Fp;
583 ppbmrm.hist_Fd[ppbmrm.nIter]=ppbmrm.Fd;
584 ppbmrm.hist_wdist[ppbmrm.nIter]=wdist;
588 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",
589 ppbmrm.nIter, tstop-tstart, ppbmrm.Fp, ppbmrm.Fd, ppbmrm.Fp-ppbmrm.Fd,
590 (ppbmrm.Fp-ppbmrm.Fd)/ppbmrm.Fp, R, ppbmrm.nCP, ppbmrm.nzA, wdist, alpha,
591 qp_cnt, gamma, tuneAlpha);
594 if (ppbmrm.nCP>=BufSize)
601 LIBBMRM_MEMCPY(prevW, W, nDim*
sizeof(
float64_t));
607 clean_icp(&icp_stats, ppbmrm, &CPList_head, &CPList_tail, H, diag_H, beta, map, cleanAfter, b, I);
611 if (ppbmrm.nCP+1 >= BufSize)
617 SGVector<float64_t> w_debug(W, nDim,
false);
620 helper->add_debug_info(primal, ppbmrm.nIter, train_error);
630 ppbmrm.hist_Fp.resize_vector(ppbmrm.nIter);
631 ppbmrm.hist_Fd.resize_vector(ppbmrm.nIter);
632 ppbmrm.hist_wdist.resize_vector(ppbmrm.nIter);
640 LIBBMRM_FREE(cp_ptr2);
652 LIBBMRM_FREE(subgrad);
653 LIBBMRM_FREE(diag_H);
655 LIBBMRM_FREE(icp_stats.ICPcounter);
656 LIBBMRM_FREE(icp_stats.ICPs);
657 LIBBMRM_FREE(icp_stats.ACPs);
658 LIBBMRM_FREE(icp_stats.H_buff);
662 LIBBMRM_FREE(beta_start);
663 LIBBMRM_FREE(beta_good);
664 LIBBMRM_FREE(I_start);
665 LIBBMRM_FREE(I_good);
668 LIBBMRM_FREE(diag_H2);
672 LIBBMRM_FREE(cp_list);
679 #endif //USE_GPL_SHOGUN
static float64_t primal_objective(SGVector< float64_t > w, CStructuredModel *model, float64_t lbda)
static float64_t dot(const bool *v1, const bool *v2, int32_t n)
Compute dot product between v1 and v2 (blas optimized)
static float64_t average_loss(SGVector< float64_t > w, CStructuredModel *model, bool is_ub=false)
all of classes and functions are contained in the shogun namespace
Matrix::Scalar max(Matrix m)
static float32_t sqrt(float32_t x)