43 MSK_deleteenv(&m_msk_env);
63 MSKidxt *qsubi,*qsubj;
69 aptrb = (MSKlidxt*) SG_MALLOC(MSKlidxt, k);
70 aptre = (MSKlidxt*) SG_MALLOC(MSKlidxt, k);
71 asub = (MSKidxt*) SG_MALLOC(MSKidxt, k);
72 bkx = (MSKboundkeye*) SG_MALLOC(MSKboundkeye, k);
73 qsubi = (MSKidxt*) SG_MALLOC(MSKidxt, Q_size);
74 qsubj = (MSKidxt*) SG_MALLOC(MSKidxt, Q_size);
93 for (int32_t i=0; i < k;i++)
102 bux[i] = MSK_INFINITY;
114 r = MSK_maketask(m_msk_env, 1, k, &task);
117 SG_ERROR(
"Could not create MOSEK task: %d\n", r)
119 r = MSK_inputdata(task,
129 SG_ERROR("Error setting input data: %d\n", r)
133 for (int32_t i=0;i<k;i++)
135 for (int32_t j=0;j<=i;j++)
139 qval[t] = G[i][j]/(1+rho);
144 r = MSK_putqobj(task, k*(k+1)/2, qsubi, qsubj, qval);
146 SG_ERROR(
"Error MSK_putqobj: %d\n", r)
159 MSK_putdouparam(task, MSK_DPAR_INTPNT_TOL_REL_GAP, 1E-8);
161 r = MSK_optimize(task);
164 SG_ERROR("Error MSK_optimize: %d\n", r)
166 MSK_getsolutionslice(task,
174 MSK_getprimalobj(task, MSK_SOL_ITR, dual_obj);
176 MSK_deletetask(&task);
204 int32_t iter, size_active;
215 float64_t expected_descent, primal_obj_b=-1, reg_master_obj;
219 float64_t proximal_term, primal_lower_bound;
231 int32_t suff_decrease_cond=0;
244 new_constraint = find_cutting_plane(&margin);
248 primal_lower_bound = 0;
249 expected_descent = -primal_obj_b;
250 initial_primal_obj = primal_obj_b;
252 SG_INFO(
"Running CCCP inner loop solver: ")
254 while ((!suff_decrease_cond) && (expected_descent<-m_eps) && (iter<m_max_iter))
264 dXc[size_active - 1] = new_constraint;
274 delta[size_active-1] = margin;
276 alpha[size_active-1] = 0.0;
278 idle[size_active-1] = 0;
290 G = SG_REALLOC(
float64_t*, G, size_active-1, size_active);
291 G[size_active - 1] = NULL;
292 for (
index_t j=0; j < size_active;j++)
294 G[j] = SG_REALLOC(
float64_t, G[j], size_active-1, size_active);
296 for (
index_t j=0; j < size_active-1; j++)
298 G[size_active-1][j] = dXc[size_active-1].sparse_dot(dXc[j]);
299 G[j][size_active-1] = G[size_active-1][j];
301 G[size_active-1][size_active-1] = dXc[size_active-1].sparse_dot(dXc[size_active-1]);
306 gammaG0[size_active-1] = dXc[size_active-1].dense_dot(1.0, w_b.
vector, w_b.
vlen, 0);
310 for (
index_t i = 0; i < size_active; i++)
311 gammaG0[i] = dXc[i].dense_dot(1.0, w_b.
vector, w_b.
vlen, 0);
315 for (
index_t i = 0; i < size_active; i++)
320 proximal_rhs[i] = delta[i] - rho/(1+rho)*gammaG0[i];
323 proximal_rhs[i] = (1+rho)*delta[i] - rho*gammaG0[i];
326 SG_ERROR(
"Invalid QPType: %d\n", m_qp_type)
335 mosek_qp_optimize(G, proximal_rhs.
vector, alpha.
vector, size_active, &dual_obj, rho);
356 SG_ERROR(
"Invalid QPType: %d\n", m_qp_type)
365 for (
index_t j = 0; j < size_active; j++)
367 if (alpha[j]>m_C*m_alpha_thrld)
371 for (
index_t k = 0; k < dXc[j].num_feat_entries; k++)
373 index_t idx = dXc[j].features[k].feat_index;
374 m_w[idx] += alpha[j]/(1+rho)*dXc[j].features[k].entry;
383 for (int32_t j=0;j<size_active;j++)
384 dual_obj -= proximal_rhs[j]/(1+rho)*alpha[j];
393 for (
index_t j = 0; j < size_active; j++)
395 sigma_k += alpha[j]*cut_error[j];
402 for (
index_t j = 0; j < size_active; j++)
403 SG_DEBUG(
"alpha[%d]: %.8g, cut_error[%d]: %.8g\n", j, alpha[j], j, cut_error[j])
404 SG_DEBUG(
"sigma_k: %.8g\n", sigma_k)
405 SG_DEBUG(
"alphasum: %.8g\n", alphasum)
409 for (
index_t j = 0; j < size_active; j++)
411 if (alpha[j]<m_alpha_thrld*m_C)
417 new_constraint = find_cutting_plane(&margin);
423 SG_DEBUG(
"ITER PRIMAL_OBJ %.4f\n", primal_obj)
428 proximal_term += (
m_w[i]-w_b[i])*(
m_w[i]-w_b[i]);
430 reg_master_obj = -dual_obj+0.5*rho*temp_var/(1+rho);
431 expected_descent = reg_master_obj - primal_obj_b;
433 v_k = (reg_master_obj - proximal_term*rho/2) - primal_obj_b;
435 primal_lower_bound =
CMath::max(primal_lower_bound, reg_master_obj - 0.5*rho*(1+rho)*proximal_term);
437 SG_DEBUG(
"ITER REG_MASTER_OBJ: %.4f\n", reg_master_obj)
438 SG_DEBUG(
"ITER EXPECTED_DESCENT: %.4f\n", expected_descent)
439 SG_DEBUG(
"ITER PRIMLA_OBJ_B: %.4f\n", primal_obj_b)
441 SG_DEBUG(
"ITER ||w-w_b||^2: %.4f\n", proximal_term)
442 SG_DEBUG(
"ITER PRIMAL_LOWER_BOUND: %.4f\n", primal_lower_bound)
444 SG_DEBUG(
"ITER margin: %.4f\n", margin)
445 SG_DEBUG(
"ITER psi*-psi: %.4f\n", value-margin)
447 obj_difference = primal_obj - primal_obj_b;
449 if (primal_obj<primal_obj_b+kappa*expected_descent)
452 if ((gTd>m2*v_k)||(rho<min_rho+1E-8))
457 for (
index_t i = 0; i < size_active; i++)
460 cut_error[i] -= m_C*dXc[i].dense_dot(1.0, w_b.
vector, w_b.
vlen, 0);
464 primal_obj_b = primal_obj;
476 SG_DEBUG(
"NULL STEP: SS(ii) FAILS.\n")
486 if ((cut_error[size_active-1]>m3*last_sigma_k)&&(
CMath::abs(obj_difference)>last_z_k_norm+last_sigma_k))
488 SG_DEBUG(
"NULL STEP: NS(ii) FAILS.\n")
495 last_sigma_k = sigma_k;
496 last_z_k_norm = z_k_norm;
500 if (primal_obj_b/initial_primal_obj<1-decrease_proportion)
501 suff_decrease_cond = 1;
504 if (iter % m_cleanup_check == 0)
506 size_active = resize_cleanup(size_active, idle, alpha, delta, gammaG0, proximal_rhs, &G, dXc, cut_error);
511 SG_INFO(
" Inner loop optimization finished.\n")
513 for (
index_t j = 0; j < size_active; j++)
521 m_primal_obj = primal_obj_b;
534 new_constraint.zero();
535 for (
index_t i = 0; i < num_samples; i++)
542 new_constraint.add(result->
psi_pred);
547 new_constraint.vlen);
549 new_constraint.vlen);
553 SG_ERROR(
"model(%s) should have either of psi_computed or psi_computed_sparse"
561 *margin += result->
delta;
566 new_constraint.scale(scale);
571 for (
index_t i=0; i < psi_size; i++)
582 for (
index_t i = 0; i < psi_size; i++)
586 cut_plane.features[k].feat_index = i;
587 cut_plane.features[k].entry = new_constraint[i];
600 int32_t i, j, new_size_active;
606 while ((i<size_active)&&(idle[i]<m_idle_iter))
610 while((j<size_active)&&(idle[j]>=m_idle_iter))
613 while (j<size_active)
618 gammaG0[i] = gammaG0[j];
619 cut_error[i] = cut_error[j];
630 while((j<size_active)&&(idle[j]>=m_idle_iter))
633 for (k=i;k<size_active;k++)
635 if (G[k]!=NULL) SG_FREE(G[k]);
643 G = SG_REALLOC(
float64_t*, G, size_active, new_size_active);
644 dXc.resize_array(new_size_active);
649 while ((i<size_active)&&(idle[i]<m_idle_iter))
653 while((j<size_active)&&(idle[j]>=m_idle_iter))
656 while (j<size_active)
659 for (k=0;k<new_size_active;k++)
664 while((j<size_active)&&(idle[j]>=m_idle_iter))
668 for (k=0;k<new_size_active;k++)
669 G[k] = SG_REALLOC(
float64_t, G[k], size_active, new_size_active);
673 return new_size_active;
676 void CCCSOSVM::init()
680 m_alpha_thrld = 1E-14;
681 m_cleanup_check = 100;
689 MSKrescodee r = MSK_RES_OK;
692 #if (MSK_VERSION_MAJOR == 6)
693 r = MSK_makeenv(&m_msk_env, NULL, NULL, NULL, NULL);
694 #elif (MSK_VERSION_MAJOR == 7)
695 r = MSK_makeenv(&m_msk_env, NULL);
697 #error "Unsupported Mosek version"
702 SG_ERROR(
"Error while creating mosek env: %d\n", r)
705 r = MSK_initenv(m_msk_env);
707 SG_ERROR("Error while initializing mosek env: %d\n", r)
713 SG_ADD(&m_cleanup_check,
"m_cleanup_check",
"Cleanup after given number of iterations",
MS_NOT_AVAILABLE);
SGVector< float64_t > psi_truth
static const float64_t INFTY
infinity
virtual int32_t get_num_vectors() const =0
CStructuredModel * m_model
virtual int32_t get_dim() const =0
void scale(T alpha)
Scale vector inplace.
bool train_machine(CFeatures *data=NULL)
void set_features(CFeatures *f)
void add_to_dense(T alpha, T *vec, int32_t dim, bool abs_val=false)
Template Dynamic array class that creates an array that can be used like a list or an array...
static float64_t dot(const bool *v1, const bool *v2, int32_t n)
Compute dot product between v1 and v2 (blas optimized)
bool resize_array(int32_t n, bool exact_resize=false)
static void vec1_plus_scalar_times_vec2(T *vec1, const T scalar, const T *vec2, int32_t n)
x=x+alpha*y
T dense_dot(T alpha, T *vec, int32_t dim, T b)
Class CStructuredModel that represents the application specific model and contains most of the applic...
void set_w(SGVector< float64_t > W)
all of classes and functions are contained in the shogun namespace
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.
void scale(Matrix A, Matrix B, typename Matrix::Scalar alpha)
virtual EMachineType get_classifier_type()
SGVector< T > clone() const
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)
virtual const char * get_name() const
SGSparseVector< float64_t > psi_pred_sparse
CFeatures * get_features()