SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PrimalMosekSOSVM.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2014 Shell Hu
8  * Written (W) 2012 Fernando José Iglesias García
9  * Copyright (C) 2012 Fernando José Iglesias García
10  */
11 
12 #ifdef USE_MOSEK
13 
15 #include <shogun/lib/List.h>
18 #include <shogun/loss/HingeLoss.h>
19 
20 using namespace shogun;
21 
22 CPrimalMosekSOSVM::CPrimalMosekSOSVM()
24  po_value(0.0)
25 {
26  init();
27 }
28 
29 CPrimalMosekSOSVM::CPrimalMosekSOSVM(
30  CStructuredModel* model,
31  CStructuredLabels* labs)
32 : CLinearStructuredOutputMachine(model, labs),
33  po_value(0.0)
34 {
35  init();
36 }
37 
38 void CPrimalMosekSOSVM::init()
39 {
40  SG_ADD(&m_slacks, "slacks", "Slacks vector", MS_NOT_AVAILABLE);
41  //FIXME model selection available for SO machines
42  SG_ADD(&m_regularization, "regularization", "Regularization constant", MS_NOT_AVAILABLE);
43  SG_ADD(&m_epsilon, "epsilon", "Violation tolerance", MS_NOT_AVAILABLE);
44  SG_ADD(&m_lb, "lb", "Lower bounds", MS_NOT_AVAILABLE);
45  SG_ADD(&m_ub, "ub", "Upper bounds", MS_NOT_AVAILABLE);
46 
47  m_regularization = 1.0;
48  m_epsilon = 0.0;
49 }
50 
51 CPrimalMosekSOSVM::~CPrimalMosekSOSVM()
52 {
53 }
54 
55 bool CPrimalMosekSOSVM::train_machine(CFeatures* data)
56 {
57  SG_DEBUG("Entering CPrimalMosekSOSVM::train_machine.\n");
58  if (data)
59  set_features(data);
60 
61  CFeatures* model_features = get_features();
62  // Initialize the model for training
63  m_model->init_training();
64  // Check that the scenary is correct to start with training
65  m_model->check_training_setup();
66  SG_DEBUG("The training setup is correct.\n");
67 
68  // Dimensionality of the joint feature space
69  int32_t M = m_model->get_dim();
70  // Number of auxiliary variables in the optimization vector
71  int32_t num_aux = m_model->get_num_aux();
72  // Number of auxiliary constraints
73  int32_t num_aux_con = m_model->get_num_aux_con();
74  // Number of training examples
75  int32_t N = model_features->get_num_vectors();
76 
77  SG_DEBUG("M=%d, N =%d, num_aux=%d, num_aux_con=%d.\n", M, N, num_aux, num_aux_con);
78 
79  // Interface with MOSEK
80  CMosek* mosek = new CMosek(0, M+num_aux+N);
81  SG_REF(mosek);
82  REQUIRE(mosek->get_rescode() == MSK_RES_OK, "Mosek object could not be properly created in PrimalMosekSOSVM training.\n");
83 
84  // Initialize the terms of the optimization problem
85  SGMatrix< float64_t > A, B, C;
86  SGVector< float64_t > a, b, lb, ub;
87  m_model->init_primal_opt(m_regularization, A, a, B, b, lb, ub, C);
88 
89  REQUIRE(lb.vlen == 0 || lb.vlen == M,
90  "%s::train_machine(): lb.vlen can only be 0 or w.vlen!\n", get_name());
91 
92  REQUIRE(ub.vlen == 0 || ub.vlen == M,
93  "%s::train_machine(): ub.vlen can only be 0 or w.vlen!\n", get_name());
94 
95  if (lb.vlen == M)
96  set_lower_bounds(lb);
97 
98  if (ub.vlen == M)
99  set_upper_bounds(ub);
100 
101  SG_DEBUG("Regularization used in PrimalMosekSOSVM equal to %.2f.\n", m_regularization);
102 
103  // Input terms of the problem that do not change between iterations
104  REQUIRE(mosek->init_sosvm(M, N, num_aux, num_aux_con, C, m_lb, m_ub, A, b) == MSK_RES_OK,
105  "Mosek error in PrimalMosekSOSVM initializing SO-SVM.\n")
106 
107  // Initialize the weight vector
108  m_w = SGVector< float64_t >(M);
109  m_w.zero();
110 
111  m_slacks = SGVector< float64_t >(N);
112  m_slacks.zero();
113 
114  // Initialize the list of constraints
115  // Each element in results is a list of CResultSet with the constraints
116  // associated to each training example
117  CDynamicObjectArray* results = new CDynamicObjectArray(N);
118  SG_REF(results);
119  for ( int32_t i = 0 ; i < N ; ++i )
120  {
121  CList* list = new CList(true);
122  results->push_back(list);
123  }
124 
125  // Initialize variables used in the loop
126  int32_t num_con = num_aux_con; // number of constraints
127  int32_t old_num_con = num_con;
128  bool exception = false;
129  index_t iteration = 0;
130 
131  SGVector< float64_t > sol(M+num_aux+N);
132  sol.zero();
133 
134  SGVector< float64_t > aux(num_aux);
135 
136  do
137  {
138  SG_DEBUG("Iteration #%d: Cutting plane training with num_con=%d and old_num_con=%d.\n",
139  iteration, num_con, old_num_con);
140 
141  old_num_con = num_con;
142 
143  for ( int32_t i = 0 ; i < N ; ++i )
144  {
145  // Predict the result of the ith training example (loss-aug)
146  CResultSet* result = m_model->argmax(m_w, i);
147 
148  // Compute the loss associated with the prediction (surrogate loss, max(0, \tilde{H}))
149  float64_t slack = CHingeLoss().loss( compute_loss_arg(result) );
150  CList* cur_list = (CList*) results->get_element(i);
151 
152  // Update the list of constraints
153  if ( cur_list->get_num_elements() > 0 )
154  {
155  // Find the maximum loss within the elements of
156  // the list of constraints
157  CResultSet* cur_res = (CResultSet*) cur_list->get_first_element();
158  float64_t max_slack = -CMath::INFTY;
159 
160  while ( cur_res != NULL )
161  {
162  max_slack = CMath::max(max_slack, CHingeLoss().loss( compute_loss_arg(cur_res) ));
163 
164  SG_UNREF(cur_res);
165  cur_res = (CResultSet*) cur_list->get_next_element();
166  }
167 
168  if ( slack > max_slack + m_epsilon )
169  {
170  // The current training example is a
171  // violated constraint
172  if ( ! insert_result(cur_list, result) )
173  {
174  exception = true;
175  break;
176  }
177 
178  add_constraint(mosek, result, num_con, i);
179  ++num_con;
180  }
181  }
182  else
183  {
184  // First iteration of do ... while, add constraint
185  if ( ! insert_result(cur_list, result) )
186  {
187  exception = true;
188  break;
189  }
190 
191  add_constraint(mosek, result, num_con, i);
192  ++num_con;
193  }
194 
195  SG_UNREF(cur_list);
196  SG_UNREF(result);
197  }
198 
199  // Solve the QP
200  SG_DEBUG("Entering Mosek QP solver.\n");
201 
202  mosek->optimize(sol);
203  for ( int32_t i = 0 ; i < M+num_aux+N ; ++i )
204  {
205  if ( i < M )
206  m_w[i] = sol[i];
207  else if ( i < M+num_aux )
208  aux[i-M] = sol[i];
209  else
210  m_slacks[i-M-num_aux] = sol[i];
211  }
212 
213  SG_DEBUG("QP solved. The primal objective value is %.4f.\n", mosek->get_primal_objective_value());
214 
215  ++iteration;
216 
217  } while ( old_num_con != num_con && ! exception );
218 
219  po_value = mosek->get_primal_objective_value();
220 
221  // Free resources
222  SG_UNREF(results);
223  SG_UNREF(mosek);
224  SG_UNREF(model_features);
225  return true;
226 }
227 
228 float64_t CPrimalMosekSOSVM::compute_loss_arg(CResultSet* result) const
229 {
230  // Dimensionality of the joint feature space
231  int32_t M = m_w.vlen;
232 
233  if(result->psi_computed)
234  {
235  return SGVector< float64_t >::dot(m_w.vector, result->psi_pred.vector, M) +
236  result->delta -
237  SGVector< float64_t >::dot(m_w.vector, result->psi_truth.vector, M);
238  }
239  else if(result->psi_computed_sparse)
240  {
241  return result->psi_pred_sparse.dense_dot(1.0, m_w.vector, m_w.vlen, 0) +
242  result->delta -
243  result->psi_truth_sparse.dense_dot(1.0, m_w.vector, m_w.vlen, 0);
244  }
245  else
246  {
247  SG_ERROR("model(%s) should have either of psi_computed or psi_computed_sparse"
248  "to be set true\n", m_model->get_name());
249  return 0;
250  }
251 }
252 
253 bool CPrimalMosekSOSVM::insert_result(CList* result_list, CResultSet* result) const
254 {
255  bool succeed = result_list->insert_element(result);
256 
257  if ( ! succeed )
258  {
259  SG_PRINT("ResultSet could not be inserted in the list..."
260  "aborting training of PrimalMosekSOSVM\n");
261  }
262 
263  return succeed;
264 }
265 
266 bool CPrimalMosekSOSVM::add_constraint(
267  CMosek* mosek,
268  CResultSet* result,
269  index_t con_idx,
270  index_t train_idx) const
271 {
272  int32_t M = m_model->get_dim();
273  SGVector< float64_t > dPsi(M);
274 
275  if (result->psi_computed)
276  {
277  for ( int i = 0 ; i < M ; ++i )
278  dPsi[i] = result->psi_pred[i] - result->psi_truth[i]; // -dPsi(y)
279  }
280  else if(result->psi_computed_sparse)
281  {
282  dPsi.zero();
283  result->psi_pred_sparse.add_to_dense(1.0, dPsi.vector, dPsi.vlen);
284  result->psi_truth_sparse.add_to_dense(-1.0, dPsi.vector, dPsi.vlen);
285  }
286  else
287  {
288  SG_ERROR("model(%s) should have either of psi_computed or psi_computed_sparse"
289  "to be set true\n", m_model->get_name());
290  }
291 
292  return ( mosek->add_constraint_sosvm(dPsi, con_idx, train_idx,
293  m_model->get_num_aux(), -result->delta) == MSK_RES_OK );
294 }
295 
296 
297 float64_t CPrimalMosekSOSVM::compute_primal_objective() const
298 {
299  return po_value;
300 }
301 
302 EMachineType CPrimalMosekSOSVM::get_classifier_type()
303 {
304  return CT_PRIMALMOSEKSOSVM;
305 }
306 
307 void CPrimalMosekSOSVM::set_regularization(float64_t C)
308 {
309  m_regularization = C;
310 }
311 
313 {
314  m_epsilon = epsilon;
315 }
316 
317 void CPrimalMosekSOSVM::set_lower_bounds(SGVector< float64_t > lb)
318 {
319  m_lb = lb.clone();
320 }
321 
322 void CPrimalMosekSOSVM::set_upper_bounds(SGVector< float64_t > ub)
323 {
324  m_ub = ub.clone();
325 }
326 
327 #endif /* USE_MOSEK */

SHOGUN Machine Learning Toolbox - Documentation