15 #include <shogun/lib/external/libqp.h>
48 libqp_state_T qp_exitflag={0, 0, 0, 0}, qp_exitflag_good={0, 0, 0, 0};
49 float64_t *b, *b2, *beta, *beta_good, *beta_start, *diag_H, *diag_H2;
50 float64_t R, *Rt, **subgrad_t, *A, QPSolverTolRel, *C=NULL;
51 float64_t *prevW, *wt, alpha, alpha_start, alpha_good=0.0, Fd_alpha0=0.0;
53 floatmax_t rsum, sq_norm_W, sq_norm_Wdiff, sq_norm_prevW, eps;
54 uint32_t *I, *I2, *I_start, *I_good;
57 bmrm_ll *CPList_head, *CPList_tail, *cp_ptr, *cp_ptr2, *cp_list=NULL;
59 bool *map=NULL, tuneAlpha=
true, flag=
true;
60 bool alphaChanged=
false, isThereGoodSolution=
false;
65 uint32_t to=0, N=0, cp_i=0;
136 if (
H==NULL || A==NULL || b==NULL || beta==NULL || subgrad_t==NULL ||
137 diag_H==NULL || I==NULL || icp_stats.
ICPcounter==NULL ||
138 icp_stats.
ICPs==NULL || icp_stats.
ACPs==NULL ||
139 cp_list==NULL || prevW==NULL || wt==NULL || Rt==NULL || C==NULL ||
140 S==NULL || info==NULL || icp_stats.
H_buff==NULL)
151 for (uint32_t p=0; p<cp_models; ++p)
158 if (subgrad_t[p]==NULL || info[p]==NULL)
165 to=((p+1)*N > (uint32_t)num_feats) ? (uint32_t)num_feats : (p+1)*N;
177 memset( (
bool*) map,
true, BufSize);
196 if (beta_start==NULL || beta_good==NULL || b2==NULL || diag_H2==NULL ||
197 I_start==NULL || I_good==NULL || I2==NULL ||
H2==NULL)
208 Rt[0] = machine->
risk(subgrad_t[0], W, info[0]);
226 for (uint32_t p=1; p<cp_models; ++p)
228 Rt[p] = machine->
risk(subgrad_t[p], W, info[p]);
236 for (uint32_t p=0; p<cp_models; ++p)
242 for (uint32_t j=0; j<nDim; ++j)
244 sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
249 p3bmrm.
Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff;
254 K = (sq_norm_W == 0.0) ? 0.4 : 0.01*
CMath::sqrt(sq_norm_W);
267 SG_SDEBUG(
"%4d: tim=%.3lf, Fp=%lf, Fd=%lf, R=%lf, K=%lf, CPmodels=%d\n",
268 p3bmrm.
nIter, tstop-tstart, p3bmrm.
Fp, p3bmrm.
Fd, R, K, cp_models);
284 for (cp_i=0; cp_i<cp_models; ++cp_i)
288 for (uint32_t p=0; p<cp_models; ++p)
302 for (cp_i=0; cp_i<p3bmrm.
nCP+cp_models; ++cp_i)
306 for (uint32_t p=0; p<cp_models; ++p)
316 for (uint32_t i=0; i<p3bmrm.
nCP; ++i)
317 for (uint32_t j=0; j<cp_models; ++j)
322 for (uint32_t p=0; p<cp_models; ++p)
325 p3bmrm.
nCP+=cp_models;
330 isThereGoodSolution=
false;
332 for (uint32_t p=0; p<cp_models; ++p)
334 I[p3bmrm.
nCP-cp_models+p]=p+1;
335 beta[p3bmrm.
nCP-cp_models+p]=0.0;
344 alpha_start=alpha; alpha=0.0;
350 for (uint32_t i=0; i<p3bmrm.
nCP; ++i)
357 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
358 diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
360 for (uint32_t j=0; j<p3bmrm.
nCP; ++j)
367 qp_exitflag=libqp_splx_solver(&
get_col, diag_H2, b2, C, I2, S, beta,
371 Fd_alpha0=-qp_exitflag.QP;
377 for (uint32_t j=0; j<p3bmrm.
nCP; ++j)
386 for (uint32_t i=0; i<nDim; ++i)
387 sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]);
393 if (alpha!=alpha_start)
409 for (uint32_t i=0; i<p3bmrm.
nCP; ++i)
416 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
417 diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
419 for (uint32_t j=0; j<p3bmrm.
nCP; ++j)
424 qp_exitflag=libqp_splx_solver(&
get_col, diag_H2, b2, C, I2, S, beta,
433 for (uint32_t j=0; j<p3bmrm.
nCP; ++j)
442 for (uint32_t i=0; i<nDim; ++i)
443 sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]);
449 if (isThereGoodSolution)
454 qp_exitflag=qp_exitflag_good;
479 qp_exitflag_good=qp_exitflag;
480 isThereGoodSolution=
true;
509 for (uint32_t i=0; i<p3bmrm.
nCP; ++i)
516 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
517 diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
519 for (uint32_t j=0; j<p3bmrm.
nCP; ++j)
524 qp_exitflag=libqp_splx_solver(&
get_col, diag_H2, b2, C, I2, S, beta,
534 for (uint32_t aaa=0; aaa<p3bmrm.
nCP; ++aaa)
551 for (uint32_t j=0; j<p3bmrm.
nCP; ++j)
561 for (uint32_t p=0; p<cp_models; ++p)
563 Rt[p] = machine->
risk(subgrad_t[p], W, info[p]);
573 for (uint32_t j=0; j<nDim; ++j)
575 sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
579 p3bmrm.
Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff;
580 p3bmrm.
Fd=-qp_exitflag.QP + ((alpha*_lambda)/(_lambda + 2*alpha))*sq_norm_prevW;
585 eps=1.0-(p3bmrm.
Fd/p3bmrm.
Fp);
586 gamma=(lastFp*(1-eps)-Fd_alpha0)/(Tmax*(1-eps));
589 if ((lastFp-p3bmrm.
Fp) <= gamma)
604 if (p3bmrm.
Fp-p3bmrm.
Fd<=TolAbs)
608 if (p3bmrm.
nCP>=BufSize)
616 for (uint32_t i=0; i<nDim; ++i)
618 sq_norm_Wdiff+=(W[i]-prevW[i])*(W[i]-prevW[i]);
630 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",
631 p3bmrm.
nIter, tstop-tstart, p3bmrm.
Fp, p3bmrm.
Fd, p3bmrm.
Fp-p3bmrm.
Fd,
632 (p3bmrm.
Fp-p3bmrm.
Fd)/p3bmrm.
Fp, R, p3bmrm.
nCP, p3bmrm.
nzA, wdist, alpha,
633 qp_cnt, gamma, tuneAlpha);
636 if (p3bmrm.
nCP>=BufSize)
649 clean_icp(&icp_stats, p3bmrm, &CPList_head,
650 &CPList_tail,
H, diag_H, beta, map,
651 cleanAfter, b, I, cp_models);
655 if (p3bmrm.
nCP+1 >= BufSize)
720 for (uint32_t p=0; p<cp_models; ++p)