16 #include <shogun/lib/external/libqp.h>
25 static const uint32_t QPSolverMaxIter=0xFFFFFFFF;
31 void add_cutting_plane(
40 "add_cutting_plane: CP index %u is not free\n", free_idx)
42 LIBBMRM_MEMCPY(A+free_idx*dim, cp_data, dim*sizeof(
float64_t));
45 bmrm_ll *cp=(bmrm_ll*)LIBBMRM_CALLOC(1, bmrm_ll);
53 cp->address=A+(free_idx*dim);
61 void remove_cutting_plane(
67 bmrm_ll *cp_list_ptr=*head;
69 while(cp_list_ptr->address != icp)
71 cp_list_ptr=cp_list_ptr->next;
74 if (cp_list_ptr==*head)
77 cp_list_ptr->next->prev=NULL;
79 else if (cp_list_ptr==*tail)
82 cp_list_ptr->prev->next=NULL;
86 cp_list_ptr->prev->next=cp_list_ptr->next;
87 cp_list_ptr->next->prev=cp_list_ptr->prev;
90 map[cp_list_ptr->idx]=
true;
91 LIBBMRM_FREE(cp_list_ptr);
94 void clean_icp(ICP_stats* icp_stats,
111 bmrm_ll* cp_ptr=*head;
114 while (cp_ptr != *tail)
116 if (icp_stats->ICPcounter[tmp_idx++]>=cleanAfter)
118 icp_stats->ICPs[cntICP++]=cp_ptr->address;
122 icp_stats->ACPs[cntACP++]=tmp_idx-1;
131 uint32_t nCP_new=bmrm.nCP-cntICP;
133 for (uint32_t i=0; i<cntICP; ++i)
138 while(cp_ptr->address != icp_stats->ICPs[i])
144 remove_cutting_plane(head, tail, map, icp_stats->ICPs[i]);
146 LIBBMRM_MEMMOVE(b+tmp_idx, b+tmp_idx+1,
147 (bmrm.nCP+cp_models-tmp_idx)*
sizeof(float64_t));
148 LIBBMRM_MEMMOVE(beta+tmp_idx, beta+tmp_idx+1,
149 (bmrm.nCP-tmp_idx)*
sizeof(float64_t));
150 LIBBMRM_MEMMOVE(diag_H+tmp_idx, diag_H+tmp_idx+1,
151 (bmrm.nCP-tmp_idx)*
sizeof(float64_t));
152 LIBBMRM_MEMMOVE(I+tmp_idx, I+tmp_idx+1,
153 (bmrm.nCP-tmp_idx)*
sizeof(uint32_t));
154 LIBBMRM_MEMMOVE(icp_stats->ICPcounter+tmp_idx, icp_stats->ICPcounter+tmp_idx+1,
155 (bmrm.nCP-tmp_idx)*
sizeof(uint32_t));
159 for (uint32_t i=0; i < nCP_new; ++i)
161 for (uint32_t j=0; j < nCP_new; ++j)
163 icp_stats->H_buff[LIBBMRM_INDEX(i, j, icp_stats->maxCPs)]=
164 Hmat[LIBBMRM_INDEX(icp_stats->ACPs[i], icp_stats->ACPs[j], icp_stats->maxCPs)];
168 for (uint32_t i=0; i<nCP_new; ++i)
169 for (uint32_t j=0; j<nCP_new; ++j)
170 Hmat[LIBBMRM_INDEX(i, j, icp_stats->maxCPs)]=
171 icp_stats->H_buff[LIBBMRM_INDEX(i, j, icp_stats->maxCPs)];
181 static const float64_t *get_col( uint32_t i)
183 return( &H[ BufSize*i ] );
186 BmrmStatistics svm_bmrm_solver(
187 CDualLibQPBMSOSVM *machine,
197 bool store_train_info)
200 libqp_state_T qp_exitflag={0, 0, 0, 0};
201 float64_t *b, *beta, *diag_H, *prevW;
202 float64_t R, *subgrad, *A, QPSolverTolRel, C=1.0, wdist=0.0;
203 floatmax_t rsum, sq_norm_W, sq_norm_Wdiff=0.0;
206 CStructuredModel* model=machine->get_model();
207 uint32_t nDim=model->get_dim();
208 CSOSVMHelper* helper = NULL;
211 float64_t tstart, tstop;
213 bmrm_ll *CPList_head, *CPList_tail, *cp_ptr, *cp_ptr2, *cp_list=NULL;
214 float64_t *A_1=NULL, *A_2=NULL;
218 tstart=ttime.cur_time_diff(
false);
223 uint32_t histSize = BufSize;
234 H= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, float64_t);
245 "overflow: %u * %u > %u -- biggest possible BufSize=%u or nDim=%u\n",
250 A= (float64_t*) LIBBMRM_CALLOC(
size_t(nDim)*size_t(BufSize),
float64_t);
258 b= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
266 beta= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
274 subgrad= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t);
282 diag_H= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
290 I= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
299 icp_stats.maxCPs = BufSize;
301 icp_stats.ICPcounter= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
302 if (icp_stats.ICPcounter==NULL)
308 icp_stats.ICPs= (float64_t**) LIBBMRM_CALLOC(BufSize, float64_t*);
309 if (icp_stats.ICPs==NULL)
315 icp_stats.ACPs= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
316 if (icp_stats.ACPs==NULL)
323 icp_stats.H_buff= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, float64_t);
324 if (icp_stats.H_buff==NULL)
330 map= (
bool*) LIBBMRM_CALLOC(BufSize,
bool);
338 memset( (
bool*) map,
true, BufSize);
340 cp_list= (bmrm_ll*) LIBBMRM_CALLOC(1, bmrm_ll);
348 prevW= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t);
361 R=machine->risk(subgrad, W);
371 LIBBMRM_MEMCPY(A, subgrad, nDim*
sizeof(float64_t));
373 cp_list->address=&A[0];
383 bmrm.Fp=R+0.5*_lambda*sq_norm_W;
384 bmrm.Fd=-LIBBMRM_PLUS_INF;
386 tstop=ttime.cur_time_diff(
false);
389 SG_SINFO(
"%4d: tim=%.3lf, Fp=%lf, Fd=%lf, R=%lf\n",
390 bmrm.nIter, tstop-tstart, bmrm.Fp, bmrm.Fd, R);
393 bmrm.hist_Fp[0]=bmrm.Fp;
394 bmrm.hist_Fd[0]=bmrm.Fd;
395 bmrm.hist_wdist[0]=0.0;
397 if (store_train_info)
398 helper = machine->get_helper();
402 while (bmrm.exitflag==0)
404 tstart=ttime.cur_time_diff(
false);
411 A_2=get_cutting_plane(CPList_tail);
414 for (uint32_t i=0; i<bmrm.nCP; ++i)
416 A_1=get_cutting_plane(cp_ptr);
420 H[LIBBMRM_INDEX(bmrm.nCP, i, BufSize)]
421 = H[LIBBMRM_INDEX(i, bmrm.nCP, BufSize)]
426 A_2=get_cutting_plane(CPList_tail);
429 H[LIBBMRM_INDEX(bmrm.nCP, bmrm.nCP, BufSize)]=rsum/_lambda;
431 diag_H[bmrm.nCP]=H[LIBBMRM_INDEX(bmrm.nCP, bmrm.nCP, BufSize)];
441 SGVector<float64_t> sb(bmrm.nCP);
443 sb.vec1_plus_scalar_times_vec2(sb.vector, 1/scale, b, bmrm.nCP);
445 SGVector<float64_t> sh(bmrm.nCP);
447 sb.vec1_plus_scalar_times_vec2(sh.vector, 1/scale, diag_H, bmrm.nCP);
450 libqp_splx_solver(&get_col, sh.vector, sb.vector, &C, I, &S, beta,
451 bmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
454 qp_exitflag=libqp_splx_solver(&get_col, diag_H, b, &C, I, &S, beta,
455 bmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
458 bmrm.qp_exitflag=qp_exitflag.exitflag;
464 for (uint32_t aaa=0; aaa<bmrm.nCP; ++aaa)
466 if (beta[aaa]>epsilon)
469 icp_stats.ICPcounter[aaa]=0;
473 icp_stats.ICPcounter[aaa]+=1;
478 memset(W, 0,
sizeof(float64_t)*nDim);
480 for (uint32_t j=0; j<bmrm.nCP; ++j)
482 A_1=get_cutting_plane(cp_ptr);
488 R = machine->risk(subgrad, W);
489 add_cutting_plane(&CPList_tail, map, A,
490 find_free_idx(map, BufSize), subgrad, nDim);
496 for (uint32_t j=0; j<nDim; ++j)
498 sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
501 bmrm.Fp=R+0.5*_lambda*sq_norm_W;
502 bmrm.Fd=-qp_exitflag.QP;
506 if (bmrm.Fp - bmrm.Fd <= TolRel*LIBBMRM_ABS(bmrm.Fp))
509 if (bmrm.Fp - bmrm.Fd <= TolAbs)
512 tstop=ttime.cur_time_diff(
false);
515 SG_SINFO(
"%4d: tim=%.3lf, Fp=%lf, Fd=%lf, (Fp-Fd)=%lf, (Fp-Fd)/Fp=%lf, R=%lf, nCP=%d, nzA=%d, QPexitflag=%d\n",
516 bmrm.nIter, tstop-tstart, bmrm.Fp, bmrm.Fd, bmrm.Fp-bmrm.Fd,
517 (bmrm.Fp-bmrm.Fd)/bmrm.Fp, R, bmrm.nCP, bmrm.nzA, qp_exitflag.exitflag);
520 if (bmrm.nIter >= histSize)
523 bmrm.hist_Fp.resize_vector(histSize);
524 bmrm.hist_Fd.resize_vector(histSize);
525 bmrm.hist_wdist.resize_vector(histSize);
529 ASSERT(bmrm.nIter < histSize);
530 bmrm.hist_Fp[bmrm.nIter]=bmrm.Fp;
531 bmrm.hist_Fd[bmrm.nIter]=bmrm.Fd;
532 bmrm.hist_wdist[bmrm.nIter]=wdist;
535 LIBBMRM_MEMCPY(prevW, W, nDim*
sizeof(float64_t));
540 clean_icp(&icp_stats, bmrm, &CPList_head, &CPList_tail, H, diag_H, beta, map, cleanAfter, b, I);
545 if (bmrm.nCP+1 >= BufSize)
549 if (store_train_info)
551 float64_t info_tstart=ttime.cur_time_diff(
false);
553 SGVector<float64_t> w_info(W, nDim,
false);
556 helper->add_debug_info(primal, bmrm.nIter, train_error);
558 float64_t info_tstop=ttime.cur_time_diff(
false);
560 SG_SINFO(
"On iteration %4d, tim=%.3lf, primal=%.3lf, train_error=%lf\n", bmrm.nIter, info_tstop-info_tstart, primal, train_error);
565 if (store_train_info)
571 ASSERT(bmrm.nIter+1 <= histSize);
572 bmrm.hist_Fp.resize_vector(bmrm.nIter+1);
573 bmrm.hist_Fd.resize_vector(bmrm.nIter+1);
574 bmrm.hist_wdist.resize_vector(bmrm.nIter+1);
582 LIBBMRM_FREE(cp_ptr2);
594 LIBBMRM_FREE(subgrad);
595 LIBBMRM_FREE(diag_H);
597 LIBBMRM_FREE(icp_stats.ICPcounter);
598 LIBBMRM_FREE(icp_stats.ICPs);
599 LIBBMRM_FREE(icp_stats.ACPs);
600 LIBBMRM_FREE(icp_stats.H_buff);
605 LIBBMRM_FREE(cp_list);
612 #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)
static void vec1_plus_scalar_times_vec2(T *vec1, const T scalar, const T *vec2, int32_t n)
x=x+alpha*y
all of classes and functions are contained in the shogun namespace
void scale(Matrix A, Matrix B, typename Matrix::Scalar alpha)
Matrix::Scalar max(Matrix m)
static float32_t sqrt(float32_t x)