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) 2012 Fernando José Iglesias García
8  * Copyright (C) 2012 Fernando José Iglesias García
9  */
10 
11 #ifdef USE_MOSEK
12 
14 #include <shogun/lib/List.h>
17 #include <shogun/loss/HingeLoss.h>
18 
19 using namespace shogun;
20 
21 CPrimalMosekSOSVM::CPrimalMosekSOSVM()
23  po_value(0.0)
24 {
25  init();
26 }
27 
28 CPrimalMosekSOSVM::CPrimalMosekSOSVM(
29  CStructuredModel* model,
30  CStructuredLabels* labs)
31 : CLinearStructuredOutputMachine(model, labs),
32  po_value(0.0)
33 {
34  init();
35 }
36 
37 void CPrimalMosekSOSVM::init()
38 {
39  SG_ADD(&m_slacks, "slacks", "Slacks vector", MS_NOT_AVAILABLE);
40  //FIXME model selection available for SO machines
41  SG_ADD(&m_regularization, "regularization", "Regularization constant", MS_NOT_AVAILABLE);
42  SG_ADD(&m_epsilon, "epsilon", "Violation tolerance", MS_NOT_AVAILABLE);
43 
44  m_regularization = 1.0;
45  m_epsilon = 0.0;
46 }
47 
48 CPrimalMosekSOSVM::~CPrimalMosekSOSVM()
49 {
50 }
51 
52 bool CPrimalMosekSOSVM::train_machine(CFeatures* data)
53 {
54  SG_DEBUG("Entering CPrimalMosekSOSVM::train_machine.\n");
55  if (data)
56  set_features(data);
57 
58  CFeatures* model_features = get_features();
59  // Initialize the model for training
60  m_model->init_training();
61  // Check that the scenary is correct to start with training
62  m_model->check_training_setup();
63  SG_DEBUG("The training setup is correct.\n");
64 
65  // Dimensionality of the joint feature space
66  int32_t M = m_model->get_dim();
67  // Number of auxiliary variables in the optimization vector
68  int32_t num_aux = m_model->get_num_aux();
69  // Number of auxiliary constraints
70  int32_t num_aux_con = m_model->get_num_aux_con();
71  // Number of training examples
72  int32_t N = model_features->get_num_vectors();
73 
74  SG_DEBUG("M=%d, N =%d, num_aux=%d, num_aux_con=%d.\n", M, N, num_aux, num_aux_con);
75 
76  // Interface with MOSEK
77  CMosek* mosek = new CMosek(0, M+num_aux+N);
78  SG_REF(mosek);
79  REQUIRE(mosek->get_rescode() == MSK_RES_OK, "Mosek object could not be properly created in PrimalMosekSOSVM training.\n");
80 
81  // Initialize the terms of the optimization problem
82  SGMatrix< float64_t > A, B, C;
83  SGVector< float64_t > a, b, lb, ub;
84  m_model->init_primal_opt(m_regularization, A, a, B, b, lb, ub, C);
85 
86  SG_DEBUG("Regularization used in PrimalMosekSOSVM equal to %.2f.\n", m_regularization);
87 
88  // Input terms of the problem that do not change between iterations
89  REQUIRE(mosek->init_sosvm(M, N, num_aux, num_aux_con, C, lb, ub, A, b) == MSK_RES_OK,
90  "Mosek error in PrimalMosekSOSVM initializing SO-SVM.\n")
91 
92  // Initialize the weight vector
93  m_w = SGVector< float64_t >(M);
94  m_w.zero();
95 
96  m_slacks = SGVector< float64_t >(N);
97  m_slacks.zero();
98 
99  // Initialize the list of constraints
100  // Each element in results is a list of CResultSet with the constraints
101  // associated to each training example
102  CDynamicObjectArray* results = new CDynamicObjectArray(N);
103  SG_REF(results);
104  for ( int32_t i = 0 ; i < N ; ++i )
105  {
106  CList* list = new CList(true);
107  results->push_back(list);
108  }
109 
110  // Initialize variables used in the loop
111  int32_t num_con = num_aux_con; // number of constraints
112  int32_t old_num_con = num_con;
113  bool exception = false;
114  index_t iteration = 0;
115 
116  SGVector< float64_t > sol(M+num_aux+N);
117  sol.zero();
118 
119  SGVector< float64_t > aux(num_aux);
120 
121  do
122  {
123  SG_DEBUG("Iteration #%d: Cutting plane training with num_con=%d and old_num_con=%d.\n",
124  iteration, num_con, old_num_con);
125 
126  old_num_con = num_con;
127 
128  for ( int32_t i = 0 ; i < N ; ++i )
129  {
130  // Predict the result of the ith training example (loss-aug)
131  CResultSet* result = m_model->argmax(m_w, i);
132 
133  // Compute the loss associated with the prediction (surrogate loss, max(0, \tilde{H}))
134  float64_t slack = CHingeLoss().loss( compute_loss_arg(result) );
135  CList* cur_list = (CList*) results->get_element(i);
136 
137  // Update the list of constraints
138  if ( cur_list->get_num_elements() > 0 )
139  {
140  // Find the maximum loss within the elements of
141  // the list of constraints
142  CResultSet* cur_res = (CResultSet*) cur_list->get_first_element();
143  float64_t max_slack = -CMath::INFTY;
144 
145  while ( cur_res != NULL )
146  {
147  max_slack = CMath::max(max_slack, CHingeLoss().loss( compute_loss_arg(cur_res) ));
148 
149  SG_UNREF(cur_res);
150  cur_res = (CResultSet*) cur_list->get_next_element();
151  }
152 
153  if ( slack > max_slack + m_epsilon )
154  {
155  // The current training example is a
156  // violated constraint
157  if ( ! insert_result(cur_list, result) )
158  {
159  exception = true;
160  break;
161  }
162 
163  add_constraint(mosek, result, num_con, i);
164  ++num_con;
165  }
166  }
167  else
168  {
169  // First iteration of do ... while, add constraint
170  if ( ! insert_result(cur_list, result) )
171  {
172  exception = true;
173  break;
174  }
175 
176  add_constraint(mosek, result, num_con, i);
177  ++num_con;
178  }
179 
180  SG_UNREF(cur_list);
181  SG_UNREF(result);
182  }
183 
184  // Solve the QP
185  SG_DEBUG("Entering Mosek QP solver.\n");
186 
187  mosek->optimize(sol);
188  for ( int32_t i = 0 ; i < M+num_aux+N ; ++i )
189  {
190  if ( i < M )
191  m_w[i] = sol[i];
192  else if ( i < M+num_aux )
193  aux[i-M] = sol[i];
194  else
195  m_slacks[i-M-num_aux] = sol[i];
196  }
197 
198  SG_DEBUG("QP solved. The primal objective value is %.4f.\n", mosek->get_primal_objective_value());
199 
200  ++iteration;
201 
202  } while ( old_num_con != num_con && ! exception );
203 
204  po_value = mosek->get_primal_objective_value();
205 
206  // Free resources
207  SG_UNREF(results);
208  SG_UNREF(mosek);
209  SG_UNREF(model_features);
210  return true;
211 }
212 
213 float64_t CPrimalMosekSOSVM::compute_loss_arg(CResultSet* result) const
214 {
215  // Dimensionality of the joint feature space
216  int32_t M = m_w.vlen;
217 
218  return SGVector< float64_t >::dot(m_w.vector, result->psi_pred.vector, M) +
219  result->delta -
220  SGVector< float64_t >::dot(m_w.vector, result->psi_truth.vector, M);
221 }
222 
223 bool CPrimalMosekSOSVM::insert_result(CList* result_list, CResultSet* result) const
224 {
225  bool succeed = result_list->insert_element(result);
226 
227  if ( ! succeed )
228  {
229  SG_PRINT("ResultSet could not be inserted in the list..."
230  "aborting training of PrimalMosekSOSVM\n");
231  }
232 
233  return succeed;
234 }
235 
236 bool CPrimalMosekSOSVM::add_constraint(
237  CMosek* mosek,
238  CResultSet* result,
239  index_t con_idx,
240  index_t train_idx) const
241 {
242  int32_t M = m_model->get_dim();
243  SGVector< float64_t > dPsi(M);
244 
245  for ( int i = 0 ; i < M ; ++i )
246  dPsi[i] = result->psi_pred[i] - result->psi_truth[i]; // -dPsi(y)
247 
248  return ( mosek->add_constraint_sosvm(dPsi, con_idx, train_idx,
249  m_model->get_num_aux(), -result->delta) == MSK_RES_OK );
250 }
251 
252 
253 float64_t CPrimalMosekSOSVM::compute_primal_objective() const
254 {
255  return po_value;
256 }
257 
258 EMachineType CPrimalMosekSOSVM::get_classifier_type()
259 {
260  return CT_PRIMALMOSEKSOSVM;
261 }
262 
263 void CPrimalMosekSOSVM::set_regularization(float64_t C)
264 {
265  m_regularization = C;
266 }
267 
269 {
270  m_epsilon = epsilon;
271 }
272 
273 #endif /* USE_MOSEK */

SHOGUN Machine Learning Toolbox - Documentation