15 using namespace shogun;
42 MSK_deleteenv(&m_msk_env);
62 MSKidxt *qsubi,*qsubj;
68 aptrb = (MSKlidxt*) SG_MALLOC(MSKlidxt, k);
69 aptre = (MSKlidxt*) SG_MALLOC(MSKlidxt, k);
70 asub = (MSKidxt*) SG_MALLOC(MSKidxt, k);
71 bkx = (MSKboundkeye*) SG_MALLOC(MSKboundkeye, k);
72 qsubi = (MSKidxt*) SG_MALLOC(MSKidxt, Q_size);
73 qsubj = (MSKidxt*) SG_MALLOC(MSKidxt, Q_size);
92 for (int32_t i=0; i < k;i++)
101 bux[i] = MSK_INFINITY;
113 r = MSK_maketask(m_msk_env, 1, k, &task);
116 SG_ERROR(
"Could not create MOSEK task: %d\n", r)
118 r = MSK_inputdata(task,
128 SG_ERROR("Error setting input data: %d\n", r)
132 for (int32_t i=0;i<k;i++)
134 for (int32_t j=0;j<=i;j++)
138 qval[t] = G[i][j]/(1+rho);
143 r = MSK_putqobj(task, k*(k+1)/2, qsubi, qsubj, qval);
145 SG_ERROR(
"Error MSK_putqobj: %d\n", r)
158 MSK_putdouparam(task, MSK_DPAR_INTPNT_TOL_REL_GAP, 1E-8);
160 r = MSK_optimize(task);
163 SG_ERROR("Error MSK_optimize: %d\n", r)
165 MSK_getsolutionslice(task,
173 MSK_getprimalobj(task, MSK_SOL_ITR, dual_obj);
175 MSK_deletetask(&task);
203 int32_t iter, size_active;
214 float64_t expected_descent, primal_obj_b=-1, reg_master_obj;
218 float64_t proximal_term, primal_lower_bound;
230 int32_t suff_decrease_cond=0;
243 new_constraint = find_cutting_plane(&margin);
247 primal_lower_bound = 0;
248 expected_descent = -primal_obj_b;
249 initial_primal_obj = primal_obj_b;
251 SG_INFO(
"Running CCCP inner loop solver: ")
253 while ((!suff_decrease_cond) && (expected_descent<-m_eps) && (iter<m_max_iter))
263 dXc[size_active - 1] = new_constraint;
273 delta[size_active-1] = margin;
275 alpha[size_active-1] = 0.0;
277 idle[size_active-1] = 0;
283 cut_error[size_active-1] += (primal_obj_b - 0.5*w_b.
dot(w_b.
vector, w_b.
vector, w_b.
vlen));
289 G = SG_REALLOC(
float64_t*, G, size_active-1, size_active);
290 G[size_active - 1] = NULL;
291 for (
index_t j=0; j < size_active;j++)
293 G[j] = SG_REALLOC(
float64_t, G[j], size_active-1, size_active);
295 for (
index_t j=0; j < size_active-1; j++)
297 G[size_active-1][j] = dXc[size_active-1].sparse_dot(dXc[j]);
298 G[j][size_active-1] = G[size_active-1][j];
300 G[size_active-1][size_active-1] = dXc[size_active-1].sparse_dot(dXc[size_active-1]);
305 gammaG0[size_active-1] = dXc[size_active-1].dense_dot(1.0, w_b.
vector, w_b.
vlen, 0);
309 for (
index_t i = 0; i < size_active; i++)
310 gammaG0[i] = dXc[i].dense_dot(1.0, w_b.
vector, w_b.
vlen, 0);
314 for (
index_t i = 0; i < size_active; i++)
319 proximal_rhs[i] = delta[i] - rho/(1+rho)*gammaG0[i];
322 proximal_rhs[i] = (1+rho)*delta[i] - rho*gammaG0[i];
325 SG_ERROR(
"Invalid QPType: %d\n", m_qp_type)
334 mosek_qp_optimize(G, proximal_rhs.
vector, alpha.
vector, size_active, &dual_obj, rho);
355 SG_ERROR(
"Invalid QPType: %d\n", m_qp_type)
364 for (
index_t j = 0; j < size_active; j++)
366 if (alpha[j]>m_C*m_alpha_thrld)
370 for (
index_t k = 0; k < dXc[j].num_feat_entries; k++)
372 index_t idx = dXc[j].features[k].feat_index;
373 m_w[idx] += alpha[j]/(1+rho)*dXc[j].features[k].entry;
382 for (int32_t j=0;j<size_active;j++)
383 dual_obj -= proximal_rhs[j]/(1+rho)*alpha[j];
392 for (
index_t j = 0; j < size_active; j++)
394 sigma_k += alpha[j]*cut_error[j];
401 for (
index_t j = 0; j < size_active; j++)
402 SG_DEBUG(
"alpha[%d]: %.8g, cut_error[%d]: %.8g\n", j, alpha[j], j, cut_error[j])
403 SG_DEBUG(
"sigma_k: %.8g\n", sigma_k)
404 SG_DEBUG(
"alphasum: %.8g\n", alphasum)
408 for (
index_t j = 0; j < size_active; j++)
410 if (alpha[j]<m_alpha_thrld*m_C)
416 new_constraint = find_cutting_plane(&margin);
422 SG_DEBUG(
"ITER PRIMAL_OBJ %.4f\n", primal_obj)
427 proximal_term += (
m_w[i]-w_b[i])*(
m_w[i]-w_b[i]);
429 reg_master_obj = -dual_obj+0.5*rho*temp_var/(1+rho);
430 expected_descent = reg_master_obj - primal_obj_b;
432 v_k = (reg_master_obj - proximal_term*rho/2) - primal_obj_b;
434 primal_lower_bound =
CMath::max(primal_lower_bound, reg_master_obj - 0.5*rho*(1+rho)*proximal_term);
436 SG_DEBUG(
"ITER REG_MASTER_OBJ: %.4f\n", reg_master_obj)
437 SG_DEBUG(
"ITER EXPECTED_DESCENT: %.4f\n", expected_descent)
438 SG_DEBUG(
"ITER PRIMLA_OBJ_B: %.4f\n", primal_obj_b)
440 SG_DEBUG(
"ITER ||w-w_b||^2: %.4f\n", proximal_term)
441 SG_DEBUG(
"ITER PRIMAL_LOWER_BOUND: %.4f\n", primal_lower_bound)
443 SG_DEBUG(
"ITER margin: %.4f\n", margin)
444 SG_DEBUG(
"ITER psi*-psi: %.4f\n", value-margin)
446 obj_difference = primal_obj - primal_obj_b;
448 if (primal_obj<primal_obj_b+kappa*expected_descent)
451 if ((gTd>m2*v_k)||(rho<min_rho+1E-8))
456 for (
index_t i = 0; i < size_active; i++)
459 cut_error[i] -= m_C*dXc[i].dense_dot(1.0, w_b.
vector, w_b.
vlen, 0);
463 primal_obj_b = primal_obj;
475 SG_DEBUG(
"NULL STEP: SS(ii) FAILS.\n")
485 if ((cut_error[size_active-1]>m3*last_sigma_k)&&(
CMath::abs(obj_difference)>last_z_k_norm+last_sigma_k))
487 SG_DEBUG(
"NULL STEP: NS(ii) FAILS.\n")
494 last_sigma_k = sigma_k;
495 last_z_k_norm = z_k_norm;
499 if (primal_obj_b/initial_primal_obj<1-decrease_proportion)
500 suff_decrease_cond = 1;
503 if (iter % m_cleanup_check == 0)
505 size_active = resize_cleanup(size_active, idle, alpha, delta, gammaG0, proximal_rhs, &G, dXc, cut_error);
510 SG_INFO(
" Inner loop optimization finished.\n")
512 for (
index_t j = 0; j < size_active; j++)
520 m_primal_obj = primal_obj_b;
533 new_constraint.zero();
534 for (
index_t i = 0; i < num_samples; i++)
541 new_constraint.add(result->
psi_pred);
546 new_constraint.vlen);
548 new_constraint.vlen);
552 SG_ERROR(
"model(%s) should have either of psi_computed or psi_computed_sparse"
560 *margin += result->
delta;
565 new_constraint.scale(scale);
570 for (
index_t i=0; i < psi_size; i++)
581 for (
index_t i = 0; i < psi_size; i++)
585 cut_plane.features[k].feat_index = i;
586 cut_plane.features[k].entry = new_constraint[i];
599 int32_t i, j, new_size_active;
605 while ((i<size_active)&&(idle[i]<m_idle_iter))
609 while((j<size_active)&&(idle[j]>=m_idle_iter))
612 while (j<size_active)
617 gammaG0[i] = gammaG0[j];
618 cut_error[i] = cut_error[j];
629 while((j<size_active)&&(idle[j]>=m_idle_iter))
632 for (k=i;k<size_active;k++)
634 if (G[k]!=NULL) SG_FREE(G[k]);
642 G = SG_REALLOC(
float64_t*, G, size_active, new_size_active);
643 dXc.resize_array(new_size_active);
648 while ((i<size_active)&&(idle[i]<m_idle_iter))
652 while((j<size_active)&&(idle[j]>=m_idle_iter))
655 while (j<size_active)
658 for (k=0;k<new_size_active;k++)
663 while((j<size_active)&&(idle[j]>=m_idle_iter))
667 for (k=0;k<new_size_active;k++)
668 G[k] = SG_REALLOC(
float64_t, G[k], size_active, new_size_active);
672 return new_size_active;
675 void CCCSOSVM::init()
679 m_alpha_thrld = 1E-14;
680 m_cleanup_check = 100;
688 MSKrescodee r = MSK_RES_OK;
691 #if (MSK_VERSION_MAJOR == 6)
692 r = MSK_makeenv(&m_msk_env, NULL, NULL, NULL, NULL);
693 #elif (MSK_VERSION_MAJOR == 7)
694 r = MSK_makeenv(&m_msk_env, NULL);
696 #error "Unsupported Mosek version"
701 SG_ERROR(
"Error while creating mosek env: %d\n", r)
704 r = MSK_initenv(m_msk_env);
706 SG_ERROR("Error while initializing mosek env: %d\n", r)
712 SG_ADD(&m_cleanup_check,
"m_cleanup_check",
"Cleanup after given number of iterations",
MS_NOT_AVAILABLE);
SGVector< float64_t > psi_truth
static float64_t dot(const bool *v1, const bool *v2, int32_t n)
compute dot product between v1 and v2 (blas optimized)
static const float64_t INFTY
infinity
void add_to_dense(T alpha, T *vec, int32_t dim, bool abs_val=false)
virtual int32_t get_num_vectors() const =0
CStructuredModel * m_model
virtual int32_t get_dim() const =0
bool train_machine(CFeatures *data=NULL)
void scale(T alpha)
scale vector inplace
void set_features(CFeatures *f)
Template Dynamic array class that creates an array that can be used like a list or an array...
static T max(T a, T b)
return the maximum of two integers
SGVector< T > clone() const
bool resize_array(int32_t n, bool exact_resize=false)
Class CStructuredModel that represents the application specific model and contains most of the applic...
void set_w(SGVector< float64_t > W)
T dense_dot(T alpha, T *vec, int32_t dim, T b)
virtual CResultSet * argmax(SGVector< float64_t > w, int32_t feat_idx, bool const training=true)=0
The class Features is the base class of all feature objects.
static T min(T a, T b)
return the minimum of two integers
virtual EMachineType get_classifier_type()
SGVector< float64_t > m_w
SGVector< float64_t > psi_pred
SGSparseVector< float64_t > psi_truth_sparse
void resize_vector(int32_t n)
static float32_t sqrt(float32_t x)
x^0.5
virtual const char * get_name() const
SGSparseVector< float64_t > psi_pred_sparse
CFeatures * get_features()
static void vec1_plus_scalar_times_vec2(T *vec1, const T scalar, const T *vec2, int32_t n)
x=x+alpha*y
static T abs(T a)
return the absolute value of a number