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;
97 "overflow: %u * %u > %u -- biggest possible BufSize=%u or nDim=%u\n",
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)
492 if (beta[aaa]>epsilon)
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)
#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
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
BmrmStatistics svm_ppbm_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, bool verbose)
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 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)
CStructuredModel * get_model() const
void resize_vector(int32_t n)
Matrix::Scalar max(Matrix m)
static float32_t sqrt(float32_t x)
x^0.5
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)