14 using namespace shogun;
41 MSK_deleteenv(&m_msk_env);
61 MSKidxt *qsubi,*qsubj;
67 aptrb = (MSKlidxt*) SG_MALLOC(MSKlidxt, k);
68 aptre = (MSKlidxt*) SG_MALLOC(MSKlidxt, k);
69 asub = (MSKidxt*) SG_MALLOC(MSKidxt, k);
70 bkx = (MSKboundkeye*) SG_MALLOC(MSKboundkeye, k);
71 qsubi = (MSKidxt*) SG_MALLOC(MSKidxt, Q_size);
72 qsubj = (MSKidxt*) SG_MALLOC(MSKidxt, Q_size);
91 for (int32_t i=0; i < k;i++)
100 bux[i] = MSK_INFINITY;
112 r = MSK_maketask(m_msk_env, 1, k, &task);
115 SG_ERROR(
"Could not create MOSEK task: %d\n", r)
117 r = MSK_inputdata(task,
127 SG_ERROR("Error setting input data: %d\n", r)
131 for (int32_t i=0;i<k;i++)
133 for (int32_t j=0;j<=i;j++)
137 qval[t] = G[i][j]/(1+rho);
142 r = MSK_putqobj(task, k*(k+1)/2, qsubi, qsubj, qval);
144 SG_ERROR(
"Error MSK_putqobj: %d\n", r)
157 MSK_putdouparam(task, MSK_DPAR_INTPNT_TOL_REL_GAP, 1E-8);
159 r = MSK_optimize(task);
162 SG_ERROR("Error MSK_optimize: %d\n", r)
164 MSK_getsolutionslice(task,
172 MSK_getprimalobj(task, MSK_SOL_ITR, dual_obj);
174 MSK_deletetask(&task);
202 int32_t iter, size_active;
213 float64_t expected_descent, primal_obj_b=-1, reg_master_obj;
217 float64_t proximal_term, primal_lower_bound;
229 int32_t suff_decrease_cond=0;
242 new_constraint = find_cutting_plane(&margin);
246 primal_lower_bound = 0;
247 expected_descent = -primal_obj_b;
248 initial_primal_obj = primal_obj_b;
250 SG_INFO(
"Running CCCP inner loop solver: ")
252 while ((!suff_decrease_cond) && (expected_descent<-m_eps) && (iter<m_max_iter))
262 dXc[size_active - 1] = new_constraint;
272 delta[size_active-1] = margin;
274 alpha[size_active-1] = 0.0;
276 idle[size_active-1] = 0;
282 cut_error[size_active-1] += (primal_obj_b - 0.5*w_b.
dot(w_b.
vector, w_b.
vector, w_b.
vlen));
288 G = SG_REALLOC(
float64_t*, G, size_active-1, size_active);
289 G[size_active - 1] = NULL;
290 for (
index_t j=0; j < size_active;j++)
292 G[j] = SG_REALLOC(
float64_t, G[j], size_active-1, size_active);
294 for (
index_t j=0; j < size_active-1; j++)
296 G[size_active-1][j] = dXc[size_active-1].sparse_dot(dXc[j]);
297 G[j][size_active-1] = G[size_active-1][j];
299 G[size_active-1][size_active-1] = dXc[size_active-1].sparse_dot(dXc[size_active-1]);
304 gammaG0[size_active-1] = dXc[size_active-1].dense_dot(1.0, w_b.
vector, w_b.
vlen, 0);
308 for (
index_t i = 0; i < size_active; i++)
309 gammaG0[i] = dXc[i].dense_dot(1.0, w_b.
vector, w_b.
vlen, 0);
313 for (
index_t i = 0; i < size_active; i++)
318 proximal_rhs[i] = delta[i] - rho/(1+rho)*gammaG0[i];
321 proximal_rhs[i] = (1+rho)*delta[i] - rho*gammaG0[i];
324 SG_ERROR(
"Invalid QPType: %d\n", m_qp_type)
333 mosek_qp_optimize(G, proximal_rhs.
vector, alpha.
vector, size_active, &dual_obj, rho);
354 SG_ERROR(
"Invalid QPType: %d\n", m_qp_type)
363 for (
index_t j = 0; j < size_active; j++)
365 if (alpha[j]>m_C*m_alpha_thrld)
369 for (
index_t k = 0; k < dXc[j].num_feat_entries; k++)
371 index_t idx = dXc[j].features[k].feat_index;
372 m_w[idx] += alpha[j]/(1+rho)*dXc[j].features[k].entry;
381 for (int32_t j=0;j<size_active;j++)
382 dual_obj -= proximal_rhs[j]/(1+rho)*alpha[j];
391 for (
index_t j = 0; j < size_active; j++)
393 sigma_k += alpha[j]*cut_error[j];
400 for (
index_t j = 0; j < size_active; j++)
401 SG_DEBUG(
"alpha[%d]: %.8g, cut_error[%d]: %.8g\n", j, alpha[j], j, cut_error[j])
402 SG_DEBUG(
"sigma_k: %.8g\n", sigma_k)
403 SG_DEBUG(
"alphasum: %.8g\n", alphasum)
407 for (
index_t j = 0; j < size_active; j++)
409 if (alpha[j]<m_alpha_thrld*m_C)
415 new_constraint = find_cutting_plane(&margin);
421 SG_DEBUG(
"ITER PRIMAL_OBJ %.4f\n", primal_obj)
426 proximal_term += (
m_w[i]-w_b[i])*(
m_w[i]-w_b[i]);
428 reg_master_obj = -dual_obj+0.5*rho*temp_var/(1+rho);
429 expected_descent = reg_master_obj - primal_obj_b;
431 v_k = (reg_master_obj - proximal_term*rho/2) - primal_obj_b;
433 primal_lower_bound =
CMath::max(primal_lower_bound, reg_master_obj - 0.5*rho*(1+rho)*proximal_term);
435 SG_DEBUG(
"ITER REG_MASTER_OBJ: %.4f\n", reg_master_obj)
436 SG_DEBUG(
"ITER EXPECTED_DESCENT: %.4f\n", expected_descent)
437 SG_DEBUG(
"ITER PRIMLA_OBJ_B: %.4f\n", primal_obj_b)
439 SG_DEBUG(
"ITER ||w-w_b||^2: %.4f\n", proximal_term)
440 SG_DEBUG(
"ITER PRIMAL_LOWER_BOUND: %.4f\n", primal_lower_bound)
442 SG_DEBUG(
"ITER margin: %.4f\n", margin)
443 SG_DEBUG(
"ITER psi*-psi: %.4f\n", value-margin)
445 obj_difference = primal_obj - primal_obj_b;
447 if (primal_obj<primal_obj_b+kappa*expected_descent)
450 if ((gTd>m2*v_k)||(rho<min_rho+1E-8))
455 for (
index_t i = 0; i < size_active; i++)
458 cut_error[i] -= m_C*dXc[i].dense_dot(1.0, w_b.
vector, w_b.
vlen, 0);
462 primal_obj_b = primal_obj;
474 SG_DEBUG(
"NULL STEP: SS(ii) FAILS.\n")
484 if ((cut_error[size_active-1]>m3*last_sigma_k)&&(
CMath::abs(obj_difference)>last_z_k_norm+last_sigma_k))
486 SG_DEBUG(
"NULL STEP: NS(ii) FAILS.\n")
493 last_sigma_k = sigma_k;
494 last_z_k_norm = z_k_norm;
498 if (primal_obj_b/initial_primal_obj<1-decrease_proportion)
499 suff_decrease_cond = 1;
502 if (iter % m_cleanup_check == 0)
504 size_active = resize_cleanup(size_active, idle, alpha, delta, gammaG0, proximal_rhs, &G, dXc, cut_error);
509 SG_INFO(
" Inner loop optimization finished.\n")
511 for (
index_t j = 0; j < size_active; j++)
519 m_primal_obj = primal_obj_b;
532 new_constraint.zero();
533 for (
index_t i = 0; i < num_samples; i++)
538 new_constraint.add(result->
psi_pred);
544 *margin += result->
delta;
549 new_constraint.scale(scale);
554 for (
index_t i=0; i < psi_size; i++)
565 for (
index_t i = 0; i < psi_size; i++)
569 cut_plane.features[k].feat_index = i;
570 cut_plane.features[k].entry = new_constraint[i];
583 int32_t i, j, new_size_active;
589 while ((i<size_active)&&(idle[i]<m_idle_iter))
593 while((j<size_active)&&(idle[j]>=m_idle_iter))
596 while (j<size_active)
601 gammaG0[i] = gammaG0[j];
602 cut_error[i] = cut_error[j];
613 while((j<size_active)&&(idle[j]>=m_idle_iter))
616 for (k=i;k<size_active;k++)
618 if (G[k]!=NULL) SG_FREE(G[k]);
626 G = SG_REALLOC(
float64_t*, G, size_active, new_size_active);
627 dXc.resize_array(new_size_active);
632 while ((i<size_active)&&(idle[i]<m_idle_iter))
636 while((j<size_active)&&(idle[j]>=m_idle_iter))
639 while (j<size_active)
642 for (k=0;k<new_size_active;k++)
647 while((j<size_active)&&(idle[j]>=m_idle_iter))
651 for (k=0;k<new_size_active;k++)
652 G[k] = SG_REALLOC(
float64_t, G[k], size_active, new_size_active);
656 return new_size_active;
659 void CCCSOSVM::init()
663 m_alpha_thrld = 1E-14;
664 m_cleanup_check = 100;
672 MSKrescodee r = MSK_RES_OK;
675 #if (MSK_VERSION_MAJOR == 6)
676 r = MSK_makeenv(&m_msk_env, NULL, NULL, NULL, NULL);
677 #elif (MSK_VERSION_MAJOR == 7)
678 r = MSK_makeenv(&m_msk_env, NULL);
680 #error "Unsupported Mosek version"
685 SG_ERROR(
"Error while creating mosek env: %d\n", r)
688 r = MSK_initenv(m_msk_env);
690 SG_ERROR("Error while initializing mosek env: %d\n", r)
696 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
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
void resize_vector(int32_t n)
static float32_t sqrt(float32_t x)
x^0.5
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