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)
536 if (beta[aaa]>epsilon)
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)
#define LIBBMRM_CALLOC(x, y)
Class Time that implements a stopwatch based on either cpu time or wall clock time.
static const double * get_col(uint32_t j)
Class DualLibQPBMSOSVM that uses Bundle Methods for Regularized Risk Minimization algorithms for stru...
float64_t * get_cutting_plane(bmrm_ll *ptr)
CSOSVMHelper * get_helper() const
SGVector< float64_t > hist_wdist
virtual int32_t get_num_vectors() const =0
SGVector< float64_t > hist_Fd
static float64_t primal_objective(SGVector< float64_t > w, CStructuredModel *model, float64_t lbda)
virtual int32_t get_dim() const =0
uint32_t find_free_idx(bool *map, uint32_t size)
static const float64_t epsilon
float64_t cur_time_diff(bool verbose=false)
static const uint32_t QPSolverMaxIter
class CSOSVMHelper contains helper functions to compute primal objectives, dual objectives, average training losses, duality gaps etc. These values will be recorded to check convergence. This class is inspired by the matlab implementation of the block coordinate Frank-Wolfe SOSVM solver [1].
virtual float64_t risk(float64_t *subgrad, float64_t *W, TMultipleCPinfo *info=0, EStructRiskType rtype=N_SLACK_MARGIN_RESCALING)
SGVector< float64_t > hist_Fp
#define LIBBMRM_MEMCPY(x, y, z)
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
static float64_t dot(const bool *v1, const bool *v2, int32_t n)
Compute dot product between v1 and v2 (blas optimized)
Class CStructuredModel that represents the application specific model and contains most of the applic...
#define LIBBMRM_INDEX(ROW, COL, NUM_ROWS)
all of classes and functions are contained in the shogun namespace
virtual void add_debug_info(float64_t primal, float64_t eff_pass, float64_t train_error, float64_t dual=-1, float64_t dgap=-1)
void add_cutting_plane(bmrm_ll **tail, bool *map, float64_t *A, uint32_t free_idx, float64_t *cp_data, uint32_t dim)
The class Features is the base class of all feature objects.
BmrmStatistics svm_p3bm_solver(CDualLibQPBMSOSVM *machine, float64_t *W, float64_t TolRel, float64_t TolAbs, float64_t _lambda, uint32_t _BufSize, bool cleanICP, uint32_t cleanAfter, float64_t K, uint32_t Tmax, uint32_t cp_models, bool verbose)
CStructuredModel * get_model() const
void resize_vector(int32_t n)
static float32_t sqrt(float32_t x)
x^0.5
CFeatures * get_features()
void clean_icp(ICP_stats *icp_stats, BmrmStatistics &bmrm, bmrm_ll **head, bmrm_ll **tail, float64_t *&Hmat, float64_t *&diag_H, float64_t *&beta, bool *&map, uint32_t cleanAfter, float64_t *&b, uint32_t *&I, uint32_t cp_models)