15 #include <shogun/lib/external/libqp.h>
47 libqp_state_T qp_exitflag={0, 0, 0, 0}, qp_exitflag_good={0, 0, 0, 0};
48 float64_t *b, *b2, *beta, *beta_good, *beta_start, *diag_H, *diag_H2;
49 float64_t R, *subgrad, *A, QPSolverTolRel, C=1.0;
50 float64_t *prevW, *wt, alpha, alpha_start, alpha_good=0.0, Fd_alpha0=0.0;
52 floatmax_t rsum, sq_norm_W, sq_norm_Wdiff, sq_norm_prevW, eps;
53 uint32_t *I, *I2, *I_start, *I_good;
59 bmrm_ll *CPList_head, *CPList_tail, *cp_ptr, *cp_ptr2, *cp_list=NULL;
61 bool *map=NULL, tuneAlpha=
true, flag=
true, alphaChanged=
false, isThereGoodSolution=
false;
96 REQUIRE(BufSize < (std::numeric_limits<size_t>::max() / nDim),
97 "overflow: %u * %u > %u -- biggest possible BufSize=%u or nDim=%u\n",
98 BufSize, nDim, std::numeric_limits<size_t>::max(),
99 (std::numeric_limits<size_t>::max() / nDim),
100 (std::numeric_limits<size_t>::max() / BufSize));
127 if (
H==NULL || A==NULL || b==NULL || beta==NULL || subgrad==NULL ||
128 diag_H==NULL || I==NULL || icp_stats.
ICPcounter==NULL ||
129 icp_stats.
ICPs==NULL || icp_stats.
ACPs==NULL ||
130 cp_list==NULL || prevW==NULL || wt==NULL)
144 memset( (
bool*) map,
true, BufSize);
149 if (icp_stats.
H_buff==NULL)
172 if (beta_start==NULL || beta_good==NULL || b2==NULL || diag_H2==NULL ||
173 I_start==NULL || I_good==NULL || I2==NULL ||
H2==NULL)
184 R = machine->
risk(subgrad, W);
207 for (uint32_t j=0; j<nDim; ++j)
209 sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
212 ppbmrm.
Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff;
217 K = (sq_norm_W == 0.0) ? 0.4 : 0.01*
CMath::sqrt(sq_norm_W);
231 SG_SDEBUG(
"%4d: tim=%.3lf, Fp=%lf, Fd=%lf, R=%lf, K=%lf\n",
232 ppbmrm.
nIter, tstop-tstart, ppbmrm.
Fp, ppbmrm.
Fd, R, K);
251 for (uint32_t i=0; i<ppbmrm.
nCP; ++i)
271 beta[ppbmrm.
nCP]=0.0;
278 isThereGoodSolution=
false;
286 alpha_start=alpha; alpha=0.0;
287 beta[ppbmrm.
nCP]=0.0;
294 for (uint32_t i=0; i<ppbmrm.
nCP; ++i)
301 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
302 diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
304 for (uint32_t j=0; j<ppbmrm.
nCP; ++j)
310 qp_exitflag=libqp_splx_solver(&
get_col, diag_H2, b2, &C, I2, &S, beta,
314 Fd_alpha0=-qp_exitflag.QP;
317 for (uint32_t i=0; i<nDim; ++i)
322 for (uint32_t j=0; j<ppbmrm.
nCP; ++j)
326 rsum+=A_1[i]*beta[j];
329 wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha);
334 for (uint32_t i=0; i<nDim; ++i)
335 sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]);
341 if (alpha!=alpha_start)
354 beta[ppbmrm.
nCP]=0.0;
359 for (uint32_t i=0; i<ppbmrm.
nCP; ++i)
366 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
367 diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
369 for (uint32_t j=0; j<ppbmrm.
nCP; ++j)
374 qp_exitflag=libqp_splx_solver(&
get_col, diag_H2, b2, &C, I2, &S, beta,
380 for (uint32_t i=0; i<nDim; ++i)
385 for (uint32_t j=0; j<ppbmrm.
nCP; ++j)
389 rsum+=A_1[i]*beta[j];
392 wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha);
396 for (uint32_t i=0; i<nDim; ++i)
397 sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]);
404 if (isThereGoodSolution)
409 qp_exitflag=qp_exitflag_good;
434 qp_exitflag_good=qp_exitflag;
435 isThereGoodSolution=
true;
464 for (uint32_t i=0; i<ppbmrm.
nCP; ++i)
471 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
472 diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
474 for (uint32_t j=0; j<ppbmrm.
nCP; ++j)
479 qp_exitflag=libqp_splx_solver(&
get_col, diag_H2, b2, &C, I2, &S, beta,
490 for (uint32_t aaa=0; aaa<ppbmrm.
nCP; ++aaa)
504 for (uint32_t i=0; i<nDim; ++i)
509 for (uint32_t j=0; j<ppbmrm.
nCP; ++j)
513 rsum+=A_1[i]*beta[j];
516 W[i]=(2*alpha*prevW[i]-rsum)/(_lambda+2*alpha);
520 R = machine->
risk(subgrad, W);
529 for (uint32_t j=0; j<nDim; ++j)
531 sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
535 ppbmrm.
Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff;
536 ppbmrm.
Fd=-qp_exitflag.QP+((alpha*_lambda)/(_lambda + 2*alpha))*sq_norm_prevW;
541 eps=1.0-(ppbmrm.
Fd/ppbmrm.
Fp);
542 gamma=(lastFp*(1-eps)-Fd_alpha0)/(Tmax*(1-eps));
545 if ((lastFp-ppbmrm.
Fp) <= gamma)
560 if (ppbmrm.
Fp-ppbmrm.
Fd<=TolAbs)
564 if (ppbmrm.
nCP>=BufSize)
572 for (uint32_t i=0; i<nDim; ++i)
574 sq_norm_Wdiff+=(W[i]-prevW[i])*(W[i]-prevW[i]);
586 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",
587 ppbmrm.
nIter, tstop-tstart, ppbmrm.
Fp, ppbmrm.
Fd, ppbmrm.
Fp-ppbmrm.
Fd,
588 (ppbmrm.
Fp-ppbmrm.
Fd)/ppbmrm.
Fp, R, ppbmrm.
nCP, ppbmrm.
nzA, wdist, alpha,
589 qp_cnt, gamma, tuneAlpha);
592 if (ppbmrm.
nCP>=BufSize)
605 clean_icp(&icp_stats, ppbmrm, &CPList_head, &CPList_tail,
H, diag_H, beta, map, cleanAfter, b, I);
609 if (ppbmrm.
nCP+1 >= BufSize)