18 using namespace shogun;
63 MSKidxt *qsubi,*qsubj;
70 aptrb = (MSKlidxt*) SG_MALLOC(MSKlidxt, k);
71 aptre = (MSKlidxt*) SG_MALLOC(MSKlidxt, k);
72 asub = (MSKidxt*) SG_MALLOC(MSKidxt, k);
73 bkx = (MSKboundkeye*) SG_MALLOC(MSKboundkeye, k);
74 qsubi = (MSKidxt*) SG_MALLOC(MSKidxt, Q_size);
75 qsubj = (MSKidxt*) SG_MALLOC(MSKidxt, Q_size);
95 for (int32_t i=0; i < k;i++)
104 bux[i] = MSK_INFINITY;
111 #if (MSK_VERSION_MAJOR == 6)
112 r = MSK_makeenv(&env, NULL, NULL, NULL, NULL);
113 #elif (MSK_VERSION_MAJOR == 7)
114 r = MSK_makeenv(&env, NULL);
116 #error "Unsupported Mosek version"
127 r = MSK_initenv(env);
132 r = MSK_maketask(env,1,k,&task);
140 r = MSK_inputdata(task,
154 for (int32_t i=0;i<k;i++) {
155 for (int32_t j=0;j<=i;j++) {
158 qval[t] = G[i][j]/(1+rho);
163 r = MSK_putqobj(task, k*(k+1)/2, qsubi,qsubj,qval);
177 MSK_putdouparam(task, MSK_DPAR_INTPNT_TOL_REL_GAP, 1E-14);
180 r = MSK_optimize(task);
184 MSK_getsolutionslice(task,
191 MSK_getprimalobj(task, MSK_SOL_ITR, dual_obj);
193 MSK_deletetask(&task);
224 int32_t iter, size_active;
235 float64_t expected_descent, primal_obj_b=-1, reg_master_obj;
239 float64_t proximal_term, primal_lower_bound;
251 int32_t suff_decrease_cond=0;
264 new_constraint = find_cutting_plane(&margin);
268 primal_lower_bound = 0;
269 expected_descent = -primal_obj_b;
270 initial_primal_obj = primal_obj_b;
272 SG_INFO(
"Running CCCP inner loop solver: ")
274 while ((!suff_decrease_cond)&&(expected_descent<-m_eps)&&(iter<m_max_iter))
284 dXc[size_active - 1] = new_constraint;
294 delta[size_active-1] = margin;
296 alpha[size_active-1] = 0.0;
298 idle[size_active-1] = 0;
304 cut_error[size_active-1] += (primal_obj_b - 0.5*w_b.
dot(w_b.
vector, w_b.
vector, w_b.
vlen));
310 G = SG_REALLOC(
float64_t*, G, size_active-1, size_active);
311 G[size_active - 1] = NULL;
312 for (
index_t j=0; j < size_active;j++)
314 G[j] = SG_REALLOC(
float64_t, G[j], size_active-1, size_active);
316 for (
index_t j=0; j < size_active-1; j++)
318 G[size_active-1][j] = dXc[size_active-1].sparse_dot(dXc[j]);
319 G[j][size_active-1] = G[size_active-1][j];
321 G[size_active-1][size_active-1] = dXc[size_active-1].sparse_dot(dXc[size_active-1]);
326 gammaG0[size_active-1] = dXc[size_active-1].dense_dot(1.0, w_b.
vector, w_b.
vlen, 0);
330 for (
index_t i = 0; i < size_active; i++)
331 gammaG0[i] = dXc[i].dense_dot(1.0, w_b.
vector, w_b.
vlen, 0);
335 for (
index_t i = 0; i < size_active; i++)
340 proximal_rhs[i] = delta[i] - rho/(1+rho)*gammaG0[i];
343 proximal_rhs[i] = (1+rho)*delta[i] - rho*gammaG0[i];
346 SG_ERROR(
"Invalid QPType: %d\n", m_qp_type)
355 mosek_qp_optimize(G, proximal_rhs.
vector, alpha.
vector, size_active, &dual_obj, rho);
376 SG_ERROR(
"Invalid QPType: %d\n", m_qp_type)
385 for (
index_t j = 0; j < size_active; j++)
387 if (alpha[j]>m_C*m_alpha_thrld)
391 for (
index_t k = 0; k < dXc[j].num_feat_entries; k++)
393 index_t idx = dXc[j].features[k].feat_index;
394 m_w[idx] += alpha[j]/(1+rho)*dXc[j].features[k].entry;
403 for (int32_t j=0;j<size_active;j++)
404 dual_obj -= proximal_rhs[j]/(1+rho)*alpha[j];
413 for (
index_t j = 0; j < size_active; j++)
415 sigma_k += alpha[j]*cut_error[j];
422 for (
index_t j = 0; j < size_active; j++)
423 SG_DEBUG(
"alpha[%d]: %.8g, cut_error[%d]: %.8g\n", j, alpha[j], j, cut_error[j])
424 SG_DEBUG(
"sigma_k: %.8g\n", sigma_k)
425 SG_DEBUG(
"alphasum: %.8g\n", alphasum)
429 for (
index_t j = 0; j < size_active; j++)
431 if (alpha[j]<m_alpha_thrld*m_C)
437 new_constraint = find_cutting_plane(&margin);
443 SG_DEBUG(
"ITER PRIMAL_OBJ %.4f\n", primal_obj)
448 proximal_term += (
m_w[i]-w_b[i])*(
m_w[i]-w_b[i]);
450 reg_master_obj = -dual_obj+0.5*rho*temp_var/(1+rho);
451 expected_descent = reg_master_obj - primal_obj_b;
453 v_k = (reg_master_obj - proximal_term*rho/2) - primal_obj_b;
455 primal_lower_bound =
CMath::max(primal_lower_bound, reg_master_obj - 0.5*rho*(1+rho)*proximal_term);
457 SG_DEBUG(
"ITER REG_MASTER_OBJ: %.4f\n", reg_master_obj)
458 SG_DEBUG(
"ITER EXPECTED_DESCENT: %.4f\n", expected_descent)
459 SG_DEBUG(
"ITER PRIMLA_OBJ_B: %.4f\n", primal_obj_b)
461 SG_DEBUG(
"ITER ||w-w_b||^2: %.4f\n", proximal_term)
462 SG_DEBUG(
"ITER PRIMAL_LOWER_BOUND: %.4f\n", primal_lower_bound)
464 SG_DEBUG(
"ITER margin: %.4f\n", margin)
465 SG_DEBUG(
"ITER psi*-psi: %.4f\n", value-margin)
467 obj_difference = primal_obj - primal_obj_b;
469 if (primal_obj<primal_obj_b+kappa*expected_descent)
472 if ((gTd>m2*v_k)||(rho<min_rho+1E-8))
477 for (
index_t i = 0; i < size_active; i++)
480 cut_error[i] -= m_C*dXc[i].dense_dot(1.0, w_b.
vector, w_b.
vlen, 0);
484 primal_obj_b = primal_obj;
496 SG_DEBUG(
"NULL STEP: SS(ii) FAILS.\n")
506 if ((cut_error[size_active-1]>m3*last_sigma_k)&&(
CMath::abs(obj_difference)>last_z_k_norm+last_sigma_k))
508 SG_DEBUG(
"NULL STEP: NS(ii) FAILS.\n")
515 last_sigma_k = sigma_k;
516 last_z_k_norm = z_k_norm;
520 if (primal_obj_b/initial_primal_obj<1-decrease_proportion)
521 suff_decrease_cond = 1;
524 if (iter % m_cleanup_check == 0)
526 size_active = resize_cleanup(size_active, idle, alpha, delta, gammaG0, proximal_rhs, &G, dXc, cut_error);
531 SG_INFO(
" Inner loop optimization finished.\n")
533 for (
index_t j = 0; j < size_active; j++)
541 m_primal_obj = primal_obj_b;
554 new_constraint.zero();
555 for (
index_t i = 0; i < num_samples; i++)
560 new_constraint.add(result->
psi_pred);
566 *margin += result->
delta;
571 new_constraint.scale(scale);
576 for (
index_t i=0; i < psi_size; i++)
587 for (
index_t i = 0; i < psi_size; i++)
591 cut_plane.features[k].feat_index = i;
592 cut_plane.features[k].entry = new_constraint[i];
605 int32_t i, j, new_size_active;
611 while ((i<size_active)&&(idle[i]<m_idle_iter))
615 while((j<size_active)&&(idle[j]>=m_idle_iter))
618 while (j<size_active)
623 gammaG0[i] = gammaG0[j];
624 cut_error[i] = cut_error[j];
635 while((j<size_active)&&(idle[j]>=m_idle_iter))
638 for (k=i;k<size_active;k++)
640 if (G[k]!=NULL) free(G[k]);
648 G = SG_REALLOC(
float64_t*, G, size_active, new_size_active);
649 dXc.resize_array(new_size_active);
654 while ((i<size_active)&&(idle[i]<m_idle_iter))
658 while((j<size_active)&&(idle[j]>=m_idle_iter))
661 while (j<size_active)
664 for (k=0;k<new_size_active;k++)
669 while((j<size_active)&&(idle[j]>=m_idle_iter))
673 for (k=0;k<new_size_active;k++)
674 G[k] = SG_REALLOC(
float64_t, G[k], size_active, new_size_active);
678 return new_size_active;
681 void CCCSOSVM::init()
685 m_alpha_thrld = 1E-14;
686 m_cleanup_check = 100;
696 SG_ADD(&m_cleanup_check,
"m_cleanup_check",
"Cleanup after given number of iterations",
MS_NOT_AVAILABLE);