QpSolver.h
Go to the documentation of this file.
1 //===========================================================================
2 /*!
3  *
4  *
5  * \brief General and specialized quadratic program classes and a generic solver.
6  *
7  *
8  *
9  * \author T. Glasmachers, O.Krause
10  * \date 2007-2013
11  *
12  *
13  * \par Copyright 1995-2015 Shark Development Team
14  *
15  * <BR><HR>
16  * This file is part of Shark.
17  * <http://image.diku.dk/shark/>
18  *
19  * Shark is free software: you can redistribute it and/or modify
20  * it under the terms of the GNU Lesser General Public License as published
21  * by the Free Software Foundation, either version 3 of the License, or
22  * (at your option) any later version.
23  *
24  * Shark is distributed in the hope that it will be useful,
25  * but WITHOUT ANY WARRANTY; without even the implied warranty of
26  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27  * GNU Lesser General Public License for more details.
28  *
29  * You should have received a copy of the GNU Lesser General Public License
30  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
31  *
32  */
33 //===========================================================================
34 
35 
36 #ifndef SHARK_ALGORITHMS_QP_QPSOLVER_H
37 #define SHARK_ALGORITHMS_QP_QPSOLVER_H
38 
39 #include <shark/Core/Timer.h>
41 #include <shark/Data/Dataset.h>
42 
43 namespace shark{
44 
45 /// \brief Most gneral problem formualtion, needs to be configured by hand.
46 template<class MatrixT>
48 public:
49  typedef MatrixT MatrixType;
50  typedef typename MatrixType::QpFloatType QpFloatType;
51 
52  //Setup only using kernel matrix, labels and regularization parameter
54  : quadratic(quadratic)
55  , linear(quadratic.size(),0)
56  , alpha(quadratic.size(),0)
57  , diagonal(quadratic.size())
58  , permutation(quadratic.size())
59  , boxMin(quadratic.size(),0)
60  , boxMax(quadratic.size(),0)
61  {
62  for(std::size_t i = 0; i!= dimensions(); ++i){
63  permutation[i] = i;
64  diagonal(i) = quadratic.entry(i, i);
65  }
66  }
67 
68  /// \brief constructor which initializes a C-SVM problem with weighted datapoints and different regularizers for every class
70  MatrixType& quadratic,
71  Data<unsigned int> const& labels,
72  Data<double> const& weights,
73  RealVector const& regularizers
74  ): quadratic(quadratic)
75  , linear(quadratic.size())
76  , alpha(quadratic.size(),0)
77  , diagonal(quadratic.size())
78  , permutation(quadratic.size())
79  , boxMin(quadratic.size())
80  , boxMax(quadratic.size())
81  {
82  SIZE_CHECK(dimensions() == linear.size());
83  SIZE_CHECK(dimensions() == quadratic.size());
84  SIZE_CHECK(dimensions() == labels.numberOfElements());
85  SIZE_CHECK(dimensions() == weights.numberOfElements());
86  SIZE_CHECK(regularizers.size() > 0);
87  SIZE_CHECK(regularizers.size() <= 2);
88 
89  double Cn = regularizers[0];
90  double Cp = regularizers[0];
91  if(regularizers.size() == 2)
92  Cp = regularizers[1];
93 
94  for(std::size_t i = 0; i!= dimensions(); ++i){
95  unsigned int label = labels.element(i);
96  double weight = weights.element(i);
97  permutation[i] = i;
98  diagonal(i) = quadratic.entry(i, i);
99  linear(i) = label? 1.0:-1.0;
100  boxMin(i) = label? 0.0:-Cn*weight;
101  boxMax(i) = label? Cp*weight : 0.0;
102  }
103  }
104 
105  std::size_t dimensions()const{
106  return quadratic.size();
107  }
108 
109  /// exchange two variables via the permutation
110  void flipCoordinates(std::size_t i, std::size_t j)
111  {
112  if (i == j) return;
113 
114  // notify the matrix cache
115  quadratic.flipColumnsAndRows(i, j);
116  std::swap( alpha[i], alpha[j]);
117  std::swap( linear[i], linear[j]);
118  std::swap( diagonal[i], diagonal[j]);
119  std::swap( boxMin[i], boxMin[j]);
120  std::swap( boxMax[i], boxMax[j]);
122  }
123 
124  /// \brief Scales all box constraints by a constant factor and adapts the solution by scaling it by the same factor.
125  void scaleBoxConstraints(double factor){
126  alpha *= factor;
127  boxMin *=factor;
128  boxMax *=factor;
129  }
130 
131 
132  /// representation of the quadratic part of the objective function
133  MatrixType& quadratic;
134 
135  ///\brief Linear part of the problem
136  RealVector linear;
137 
138  /// Solution candidate
139  RealVector alpha;
140 
141  /// diagonal matrix entries
142  /// The diagonal array is of fixed size and not subject to shrinking.
143  RealVector diagonal;
144 
145  /// permutation of the variables alpha, gradient, etc.
146  std::vector<std::size_t> permutation;
147 
148  RealVector boxMin;
149 
150  RealVector boxMax;
151 };
152 
153 ///\brief Boxed problem for alpha in [lower,upper]^n and equality constraints.
154 ///
155 ///It is assumed for the initial alpha value that there exists a sum to one constraint and lower <= 1/n <= upper
156 template<class MatrixT>
158 public:
159  typedef MatrixT MatrixType;
160  typedef typename MatrixType::QpFloatType QpFloatType;
161 
162  //Setup only using kernel matrix, labels and regularization parameter
163  BoxedSVMProblem(MatrixType& quadratic, RealVector const& linear, double lower, double upper)
164  : quadratic(quadratic)
165  , linear(linear)
166  , alpha(quadratic.size(),1.0/quadratic.size())
167  , diagonal(quadratic.size())
168  , permutation(quadratic.size())
169  , m_lower(lower)
170  , m_upper(upper)
171  {
172  SIZE_CHECK(dimensions() == linear.size());
173  SIZE_CHECK(dimensions() == quadratic.size());
174 
175  for(std::size_t i = 0; i!= dimensions(); ++i){
176  permutation[i] = i;
177  diagonal(i) = quadratic.entry(i, i);
178  }
179  }
180 
181  std::size_t dimensions()const{
182  return quadratic.size();
183  }
184 
185  double boxMin(std::size_t i)const{
186  return m_lower;
187  }
188  double boxMax(std::size_t i)const{
189  return m_upper;
190  }
191 
192  /// representation of the quadratic part of the objective function
193  MatrixType& quadratic;
194 
195  ///\brief Linear part of the problem
196  RealVector linear;
197 
198  /// Solution candidate
199  RealVector alpha;
200 
201  /// diagonal matrix entries
202  /// The diagonal array is of fixed size and not subject to shrinking.
203  RealVector diagonal;
204 
205  /// exchange two variables via the permutation
206  void flipCoordinates(std::size_t i, std::size_t j)
207  {
208  if (i == j) return;
209 
210  // notify the matrix cache
211  quadratic.flipColumnsAndRows(i, j);
212  std::swap( alpha[i], alpha[j]);
213  std::swap( linear[i], linear[j]);
214  std::swap( diagonal[i], diagonal[j]);
216  }
217 
218  /// \brief Scales all box constraints by a constant factor and adapts the solution by scaling it by the same factor.
219  void scaleBoxConstraints(double factor){
220  m_lower*=factor;
221  m_upper*=factor;
222  alpha *= factor;
223  }
224 
225  /// permutation of the variables alpha, gradient, etc.
226  std::vector<std::size_t> permutation;
227 private:
228  double m_lower;
229  double m_upper;
230 };
231 
232 
233 /// \brief Problem formulation for binary C-SVM problems
234 template<class MatrixT>
236 public:
237  typedef MatrixT MatrixType;
238  typedef typename MatrixType::QpFloatType QpFloatType;
239 
240  /// \brief Setup only using kernel matrix, labels and regularization parameter
241  CSVMProblem(MatrixType& quadratic, Data<unsigned int> const& labels, double C)
242  : quadratic(quadratic)
243  , linear(quadratic.size())
244  , alpha(quadratic.size(),0)
245  , diagonal(quadratic.size())
246  , permutation(quadratic.size())
247  , positive(quadratic.size())
248  , m_Cp(C)
249  , m_Cn(C)
250  {
251  SIZE_CHECK(dimensions() == linear.size());
252  SIZE_CHECK(dimensions() == quadratic.size());
253  SIZE_CHECK(dimensions() == labels.numberOfElements());
254 
255  for(std::size_t i = 0; i!= dimensions(); ++i){
256  permutation[i] = i;
257  diagonal(i) = quadratic.entry(i, i);
258  linear(i) = labels.element(i)? 1.0:-1.0;
259  positive[i] = linear(i) > 0;
260  }
261  }
262  ///\brief Setup using kernel matrix, labels and different regularization parameters for positive and negative classes
263  CSVMProblem(MatrixType& quadratic, Data<unsigned int> const& labels, RealVector const& regularizers)
264  : quadratic(quadratic)
265  , linear(quadratic.size())
266  , alpha(quadratic.size(),0)
267  , diagonal(quadratic.size())
268  , permutation(quadratic.size())
269  , positive(quadratic.size())
270  {
271  SIZE_CHECK(dimensions() == linear.size());
272  SIZE_CHECK(dimensions() == quadratic.size());
273  SIZE_CHECK(dimensions() == labels.numberOfElements());
274 
275  SIZE_CHECK(regularizers.size() > 0);
276  SIZE_CHECK(regularizers.size() <= 2);
277  m_Cp = m_Cn = regularizers[0];
278  if(regularizers.size() == 2)
279  m_Cp = regularizers[1];
280 
281 
282  for(std::size_t i = 0; i!= dimensions(); ++i){
283  permutation[i] = i;
284  diagonal(i) = quadratic.entry(i, i);
285  linear(i) = labels.element(i)? 1.0:-1.0;
286  positive[i] = linear(i) > 0;
287  }
288  }
289 
290  //Setup with changed linear part
291  CSVMProblem(MatrixType& quadratic, RealVector linear, Data<unsigned int> const& labels, double C)
292  : quadratic(quadratic)
293  , linear(linear)
294  , alpha(quadratic.size(),0)
295  , diagonal(quadratic.size())
296  , permutation(quadratic.size())
297  , positive(quadratic.size())
298  , m_Cp(C)
299  , m_Cn(C)
300  {
301  SIZE_CHECK(dimensions() == quadratic.size());
302  SIZE_CHECK(dimensions() == linear.size());
303  SIZE_CHECK(dimensions() == labels.numberOfElements());
304 
305  for(std::size_t i = 0; i!= dimensions(); ++i){
306  permutation[i] = i;
307  diagonal(i) = quadratic.entry(i, i);
308  positive[i] = labels.element(i) ? 1: 0;
309  }
310  }
311 
312  std::size_t dimensions()const{
313  return quadratic.size();
314  }
315 
316  double boxMin(std::size_t i)const{
317  return positive[i] ? 0.0 : -m_Cn;
318  }
319  double boxMax(std::size_t i)const{
320  return positive[i] ? m_Cp : 0.0;
321  }
322 
323  /// representation of the quadratic part of the objective function
324  MatrixType& quadratic;
325 
326  ///\brief Linear part of the problem
327  RealVector linear;
328 
329  /// Solution candidate
330  RealVector alpha;
331 
332  /// diagonal matrix entries
333  /// The diagonal array is of fixed size and not subject to shrinking.
334  RealVector diagonal;
335 
336  /// exchange two variables via the permutation
337  void flipCoordinates(std::size_t i, std::size_t j)
338  {
339  if (i == j) return;
340 
341  // notify the matrix cache
342  quadratic.flipColumnsAndRows(i, j);
343  std::swap( linear[i], linear[j]);
344  std::swap( alpha[i], alpha[j]);
345  std::swap( diagonal[i], diagonal[j]);
347  std::swap( positive[i], positive[j]);
348  }
349 
350  /// \brief Scales all box constraints by a constant factor and adapts the solution by scaling it by the same factor.
351  void scaleBoxConstraints(double factor, double variableScalingFactor){
352  bool sameFactor = factor == variableScalingFactor;
353  double newCp = m_Cp*factor;
354  double newCn = m_Cn*factor;
355  for(std::size_t i = 0; i != dimensions(); ++i){
356  if(sameFactor && alpha(i)== m_Cp)
357  alpha(i) = newCp;
358  else if(sameFactor && alpha(i) == -m_Cn)
359  alpha(i) = -newCn;
360  else
361  alpha(i) *= variableScalingFactor;
362  }
363  m_Cp = newCp;
364  m_Cn = newCn;
365  }
366 
367  /// permutation of the variables alpha, gradient, etc.
368  std::vector<std::size_t> permutation;
369 private:
370  ///\brief whether the label of the point is positive
371  std::vector<char> positive;
372 
373  ///\brief Regularization constant of the positive class
374  double m_Cp;
375  ///\brief Regularization constant of the negative class
376  double m_Cn;
377 };
378 
383  AlphaDeactivated = 3//also: AlphaUpperBound and AlphaLowerBound
384 };
385 
386 ///
387 /// \brief Quadratic program solver
388 ///
389 /// todo: new documentation
390 template<class Problem, class SelectionStrategy = typename Problem::PreferedSelectionStrategy >
391 class QpSolver
392 {
393 public:
395  Problem& problem
396  ):m_problem(problem){}
397 
398  void solve(
399  QpStoppingCondition& stop,
400  QpSolutionProperties* prop = NULL
401  ){
402  double start_time = Timer::now();
403  unsigned long long iter = 0;
404  unsigned long long shrinkCounter = 0;
405 
406  SelectionStrategy workingSet;
407 
408  // decomposition loop
409  for(;;){
410  //stop if iterations exceed
411  if( iter == stop.maxIterations){
412  if (prop != NULL) prop->type = QpMaxIterationsReached;
413  break;
414  }
415  //stop if the maximum running time is exceeded
416  if (stop.maxSeconds < 1e100 && (iter+1) % 1000 == 0 ){
417  double current_time = Timer::now();
418  if (current_time - start_time > stop.maxSeconds){
419  if (prop != NULL) prop->type = QpTimeout;
420  break;
421  }
422  }
423  // select a working set and check for optimality
424  std::size_t i = 0, j = 0;
425  if (workingSet(m_problem,i, j) < stop.minAccuracy){
426  m_problem.unshrink();
427  if(m_problem.checkKKT() < stop.minAccuracy){
428  if (prop != NULL) prop->type = QpAccuracyReached;
429  break;
430  }
431  m_problem.shrink(stop.minAccuracy);
432  workingSet(m_problem,i,j);
433  workingSet.reset();
434  }
435 
436  //update smo with the selected working set
437  m_problem.updateSMO(i,j);
438 
439  //do a shrinking every 1000 iterations. if variables got shrink
440  //notify working set selection
441  if(shrinkCounter == 0 && m_problem.shrink(stop.minAccuracy)){
442  shrinkCounter = std::max<std::size_t>(1000,m_problem.dimensions());
443  workingSet.reset();
444  }
445  iter++;
446  shrinkCounter--;
447  }
448 
449  if (prop != NULL)
450  {
451  double finish_time = Timer::now();
452 
453  std::size_t i = 0, j = 0;
454  prop->accuracy = workingSet(m_problem,i, j);
455  prop->value = m_problem.functionValue();
456  prop->iterations = iter;
457  prop->seconds = finish_time - start_time;
458  }
459  }
460 
461 protected:
462  Problem& m_problem;
463 };
464 
465 }
466 #endif