SHOGUN  v2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
libqp_gsmo.cpp
Go to the documentation of this file.
1 /*-----------------------------------------------------------------------
2  * libqp_gsmo.c: implementation of the Generalized SMO algorithm.
3  *
4  * DESCRIPTION
5  * The library provides function which solves the following instance of
6  * a convex Quadratic Programming task:
7  *
8  * min QP(x) := 0.5*x'*H*x + f'*x
9  * x
10  *
11  * s.t. a'*x = b
12  * LB[i] <= x[i] <= UB[i] for all i=1..n
13  *
14  * A precision of the found solution is controlled by the input argument
15  * TolKKT which defines tightness of the relaxed Karush-Kuhn-Tucker
16  * stopping conditions.
17  *
18  * INPUT ARGUMENTS
19  * get_col function which returns pointer to the i-th column of H.
20  * diag_H [float64_t n x 1] vector containing values on the diagonal of H.
21  * f [float64_t n x 1] vector.
22  * a [float64_t n x 1] Vector which must not contain zero entries.
23  * b [float64_t 1 x 1] Scalar.
24  * LB [float64_t n x 1] Lower bound; -inf is allowed.
25  * UB [float64_t n x 1] Upper bound; inf is allowed.
26  * x [float64_t n x 1] solution vector; must be feasible.
27  * n [uint32_t 1 x 1] dimension of H.
28  * MaxIter [uint32_t 1 x 1] max number of iterations.
29  * TolKKT [float64_t 1 x 1] Tightness of KKT stopping conditions.
30  * print_state print function; if == NULL it is not called.
31  *
32  * RETURN VALUE
33  * structure [libqp_state_T]
34  * .QP [1x1] Primal objective value.
35  * .exitflag [1 x 1] Indicates which stopping condition was used:
36  * -3 ... initial solution vector does not satisfy equality constraint
37  * -2 ... initial solution vector does not satisfy bounds
38  * -1 ... not enough memory
39  * 0 ... Maximal number of iterations reached: nIter >= MaxIter.
40  * 4 ... Relaxed KKT conditions satisfied.
41  * .nIter [1x1] Number of iterations.
42  *
43  * REFERENCE
44  * S.S. Keerthi, E.G. Gilbert. Convergence of a generalized SMO algorithm
45  * for SVM classier design. Technical Report CD-00-01, Control Division,
46  * Dept. of Mechanical and Production Engineering, National University
47  * of Singapore, 2000.
48  * http://citeseer.ist.psu.edu/keerthi00convergence.html
49  *
50  *
51  * Copyright (C) 2006-2008 Vojtech Franc, xfrancv@cmp.felk.cvut.cz
52  * Center for Machine Perception, CTU FEL Prague
53  *
54  * This program is free software; you can redistribute it and/or
55  * modify it under the terms of the GNU General Public
56  * License as published by the Free Software Foundation;
57  * Version 3, 29 June 2007
58  *-------------------------------------------------------------------- */
59 
60 #include <math.h>
61 #include <stdlib.h>
62 #include <stdio.h>
63 #include <string.h>
64 #include <stdint.h>
65 #include <limits.h>
66 
67 #include <shogun/lib/common.h>
69 
70 namespace shogun
71 {
72 
73 libqp_state_T libqp_gsmo_solver(const float64_t* (*get_col)(uint32_t),
74  float64_t *diag_H,
75  float64_t *f,
76  float64_t *a,
77  float64_t b,
78  float64_t *LB,
79  float64_t *UB,
80  float64_t *x,
81  uint32_t n,
82  uint32_t MaxIter,
83  float64_t TolKKT,
84  void (*print_state)(libqp_state_T state))
85 {
86  float64_t *col_u;
87  float64_t *col_v;
88  float64_t *Nabla;
89  float64_t minF_up;
90  float64_t maxF_low;
91  float64_t tau;
92  float64_t F_i;
93  float64_t tau_ub, tau_lb;
94  uint32_t i, j;
95  uint32_t u=0, v=0;
96  libqp_state_T state;
97  float64_t atx = 0.0;
98 
99  Nabla = NULL;
100 
101  /* ------------------------------------------------------------ */
102  /* Initialization */
103  /* ------------------------------------------------------------ */
104 
105  // check bounds of initial guess
106  for (i=0; i<n; i++)
107  {
108  if (x[i]>UB[i])
109  {
110  state.exitflag = -2;
111  goto cleanup;
112  }
113  if (x[i]<LB[i])
114  {
115  state.exitflag = -2;
116  goto cleanup;
117  }
118  }
119 
120  // check equality constraint
121  for (i=0; i<n; i++)
122  atx += a[i]*x[i];
123  if (fabs(b-atx)>1e-9)
124  {
125  printf("%f \ne %f\n",b,atx);
126  state.exitflag = -3;
127  goto cleanup;
128  }
129 
130  /* Nabla = H*x + f is gradient*/
131  Nabla = (float64_t*)LIBQP_CALLOC(n, sizeof(float64_t));
132  if( Nabla == NULL )
133  {
134  state.exitflag=-1;
135  goto cleanup;
136  }
137 
138  /* compute gradient */
139  for( i=0; i < n; i++ )
140  {
141  Nabla[i] += f[i];
142  if( x[i] != 0 ) {
143  col_u = (float64_t*)get_col(i);
144  for( j=0; j < n; j++ ) {
145  Nabla[j] += col_u[j]*x[i];
146  }
147  }
148  }
149 
150  if( print_state != NULL)
151  {
152  state.QP = 0;
153  for(i = 0; i < n; i++ )
154  state.QP += 0.5*(x[i]*Nabla[i]+x[i]*f[i]);
155 
156  print_state( state );
157  }
158 
159 
160  /* ------------------------------------------------------------ */
161  /* Main optimization loop */
162  /* ------------------------------------------------------------ */
163 
164  state.nIter = 0;
165  state.exitflag = 100;
166  while( state.exitflag == 100 )
167  {
168  state.nIter ++;
169 
170  /* find the most violating pair of variables */
171  minF_up = LIBQP_PLUS_INF;
172  maxF_low = -LIBQP_PLUS_INF;
173  for(i = 0; i < n; i++ )
174  {
175 
176  F_i = Nabla[i]/a[i];
177 
178  if(LB[i] < x[i] && x[i] < UB[i])
179  { /* i is from I_0 */
180  if( minF_up > F_i) { minF_up = F_i; u = i; }
181  if( maxF_low < F_i) { maxF_low = F_i; v = i; }
182  }
183  else if((a[i] > 0 && x[i] == LB[i]) || (a[i] < 0 && x[i] == UB[i]))
184  { /* i is from I_1 or I_2 */
185  if( minF_up > F_i) { minF_up = F_i; u = i; }
186  }
187  else if((a[i] > 0 && x[i] == UB[i]) || (a[i] < 0 && x[i] == LB[i]))
188  { /* i is from I_3 or I_4 */
189  if( maxF_low < F_i) { maxF_low = F_i; v = i; }
190  }
191  }
192 
193  /* check KKT conditions */
194  if( maxF_low - minF_up <= TolKKT )
195  state.exitflag = 4;
196  else
197  {
198  /* SMO update of the most violating pair */
199  col_u = (float64_t*)get_col(u);
200  col_v = (float64_t*)get_col(v);
201 
202  if( a[u] > 0 )
203  { tau_lb = (LB[u]-x[u])*a[u]; tau_ub = (UB[u]-x[u])*a[u]; }
204  else
205  { tau_ub = (LB[u]-x[u])*a[u]; tau_lb = (UB[u]-x[u])*a[u]; }
206 
207  if( a[v] > 0 )
208  { tau_lb = LIBQP_MAX(tau_lb,(x[v]-UB[v])*a[v]); tau_ub = LIBQP_MIN(tau_ub,(x[v]-LB[v])*a[v]); }
209  else
210  { tau_lb = LIBQP_MAX(tau_lb,(x[v]-LB[v])*a[v]); tau_ub = LIBQP_MIN(tau_ub,(x[v]-UB[v])*a[v]); }
211 
212  tau = (Nabla[v]/a[v]-Nabla[u]/a[u])/
213  (diag_H[u]/(a[u]*a[u]) + diag_H[v]/(a[v]*a[v]) - 2*col_u[v]/(a[u]*a[v]));
214 
215  tau = LIBQP_MIN(LIBQP_MAX(tau,tau_lb),tau_ub);
216 
217  x[u] += tau/a[u];
218  x[v] -= tau/a[v];
219 
220  /* update Nabla */
221  for(i = 0; i < n; i++ )
222  Nabla[i] += col_u[i]*tau/a[u] - col_v[i]*tau/a[v];
223 
224  }
225 
226  if( state.nIter >= MaxIter )
227  state.exitflag = 0;
228 
229  if( print_state != NULL)
230  {
231  state.QP = 0;
232  for(i = 0; i < n; i++ )
233  state.QP += 0.5*(x[i]*Nabla[i]+x[i]*f[i]);
234 
235  print_state( state );
236  }
237 
238  }
239 
240  /* compute primal objective value */
241  state.QP = 0;
242  for(i = 0; i < n; i++ )
243  state.QP += 0.5*(x[i]*Nabla[i]+x[i]*f[i]);
244 
245 cleanup:
246 
247  LIBQP_FREE(Nabla);
248 
249  return( state );
250 }
251 
252 } /* shogun namespace */
253 

SHOGUN Machine Learning Toolbox - Documentation