15 #include <shogun/lib/external/libqp.h>
31 return( &H2[ BufSize*i ] );
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, *subgrad, *A, QPSolverTolRel, C=1.0;
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;
60 bmrm_ll *CPList_head, *CPList_tail, *cp_ptr, *cp_ptr2, *cp_list=NULL;
62 bool *map=NULL, tuneAlpha=
true, flag=
true, alphaChanged=
false, isThereGoodSolution=
false;
120 if (H==NULL || A==NULL || b==NULL || beta==NULL || subgrad==NULL ||
121 diag_H==NULL || I==NULL || icp_stats.
ICPcounter==NULL ||
122 icp_stats.
ICPs==NULL || icp_stats.
ACPs==NULL ||
123 cp_list==NULL || prevW==NULL || wt==NULL)
137 memset( (
bool*) map,
true, BufSize);
142 if (icp_stats.
H_buff==NULL)
165 if (beta_start==NULL || beta_good==NULL || b2==NULL || diag_H2==NULL ||
166 I_start==NULL || I_good==NULL || I2==NULL || H2==NULL)
177 R = machine->
risk(subgrad, W);
200 for (uint32_t j=0; j<nDim; ++j)
202 sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
205 ppbmrm.
Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff;
210 K = (sq_norm_W == 0.0) ? 0.4 : 0.01*
CMath::sqrt(sq_norm_W);
224 SG_SDEBUG(
"%4d: tim=%.3lf, Fp=%lf, Fd=%lf, R=%lf, K=%lf\n",
225 ppbmrm.
nIter, tstop-tstart, ppbmrm.
Fp, ppbmrm.
Fd, R, K);
244 for (uint32_t i=0; i<ppbmrm.
nCP; ++i)
265 beta[ppbmrm.
nCP]=0.0;
271 isThereGoodSolution=
false;
279 alpha_start=alpha; alpha=0.0;
280 beta[ppbmrm.
nCP]=0.0;
287 for (uint32_t i=0; i<ppbmrm.
nCP; ++i)
294 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
295 diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
297 for (uint32_t j=0; j<ppbmrm.
nCP; ++j)
303 qp_exitflag=libqp_splx_solver(&
get_col, diag_H2, b2, &C, I2, &S, beta,
307 Fd_alpha0=-qp_exitflag.QP;
310 for (uint32_t i=0; i<nDim; ++i)
315 for (uint32_t j=0; j<ppbmrm.
nCP; ++j)
319 rsum+=A_1[i]*beta[j];
322 wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha);
327 for (uint32_t i=0; i<nDim; ++i)
328 sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]);
334 if (alpha!=alpha_start)
347 beta[ppbmrm.
nCP]=0.0;
352 for (uint32_t i=0; i<ppbmrm.
nCP; ++i)
359 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
360 diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
362 for (uint32_t j=0; j<ppbmrm.
nCP; ++j)
367 qp_exitflag=libqp_splx_solver(&
get_col, diag_H2, b2, &C, I2, &S, beta,
373 for (uint32_t i=0; i<nDim; ++i)
378 for (uint32_t j=0; j<ppbmrm.
nCP; ++j)
382 rsum+=A_1[i]*beta[j];
385 wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha);
389 for (uint32_t i=0; i<nDim; ++i)
390 sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]);
397 if (isThereGoodSolution)
402 qp_exitflag=qp_exitflag_good;
427 qp_exitflag_good=qp_exitflag;
428 isThereGoodSolution=
true;
457 for (uint32_t i=0; i<ppbmrm.
nCP; ++i)
464 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
465 diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
467 for (uint32_t j=0; j<ppbmrm.
nCP; ++j)
472 qp_exitflag=libqp_splx_solver(&
get_col, diag_H2, b2, &C, I2, &S, beta,
483 for (uint32_t aaa=0; aaa<ppbmrm.
nCP; ++aaa)
485 if (beta[aaa]>epsilon)
497 for (uint32_t i=0; i<nDim; ++i)
502 for (uint32_t j=0; j<ppbmrm.
nCP; ++j)
506 rsum+=A_1[i]*beta[j];
509 W[i]=(2*alpha*prevW[i]-rsum)/(_lambda+2*alpha);
513 R = machine->
risk(subgrad, W);
522 for (uint32_t j=0; j<nDim; ++j)
524 sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
528 ppbmrm.
Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff;
529 ppbmrm.
Fd=-qp_exitflag.QP+((alpha*_lambda)/(_lambda + 2*alpha))*sq_norm_prevW;
534 eps=1.0-(ppbmrm.
Fd/ppbmrm.
Fp);
535 gamma=(lastFp*(1-eps)-Fd_alpha0)/(Tmax*(1-eps));
538 if ((lastFp-ppbmrm.
Fp) <= gamma)
553 if (ppbmrm.
Fp-ppbmrm.
Fd<=TolAbs)
557 if (ppbmrm.
nCP>=BufSize)
565 for (uint32_t i=0; i<nDim; ++i)
567 sq_norm_Wdiff+=(W[i]-prevW[i])*(W[i]-prevW[i]);
579 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",
580 ppbmrm.
nIter, tstop-tstart, ppbmrm.
Fp, ppbmrm.
Fd, ppbmrm.
Fp-ppbmrm.
Fd,
581 (ppbmrm.
Fp-ppbmrm.
Fd)/ppbmrm.
Fp, R, ppbmrm.
nCP, ppbmrm.
nzA, wdist, alpha,
582 qp_cnt, gamma, tuneAlpha);
585 if (ppbmrm.
nCP>=BufSize)
598 clean_icp(&icp_stats, ppbmrm, &CPList_head, &CPList_tail, H, diag_H, beta, map, cleanAfter, b, I);
#define LIBBMRM_CALLOC(x, y)
Class Time that implements a stopwatch based on either cpu time or wall clock time.
static float64_t dot(const bool *v1, const bool *v2, int32_t n)
compute dot product between v1 and v2 (blas optimized)
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 float64_t average_loss(SGVector< float64_t > w, CStructuredModel *model)
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)
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)
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)