SHOGUN  v2.0.0
libqp_splx.cpp
Go to the documentation of this file.
1 /*-----------------------------------------------------------------------
3  * simplex constraints.
4  *
5  * DESCRIPTION
6  * The library provides function which solves the following instance of
8  *
9  * min QP(x):= 0.5*x'*H*x + f'*x
10  * x
11  *
12  * subject to:
13  * sum_{i in I_k} x[i] == b[k] for all k such that S[k] == 0
14  * sum_{i in I_k} x[i] <= b[k] for all k such that S[k] == 1
15  * x(i) >= 0 for all i=1:n
16  *
17  * where I_k = { i | I[i] == k}, k={1,...,m}.
18  *
19  * A precision of the found solution is controled by the input argumens
20  * MaxIter, TolAbs, QP_TH and MaxIter which define the stopping conditions:
21  *
22  * nIter >= MaxIter -> exitflag = 0 Number of iterations
23  * QP-QD <= TolAbs -> exitflag = 1 Abs. tolerance (duality gap)
24  * QP-QD <= QP*TolRel -> exitflag = 2 Relative tolerance
25  * QP <= QP_TH -> exitflag = 3 Threshold on objective value
26  *
27  * where QP and QD are primal respectively dual objective values.
28  *
29  * INPUT ARGUMENTS
30  * get_col function which returns pointer to the i-th column of H.
31  * diag_H [float64_t n x 1] vector containing values on the diagonal of H.
32  * f [float64_t n x 1] vector.
33  * b [float64_t n x 1] vector of positive numbers.
34  * I [uint16_T n x 1] vector containing numbers 1...m.
35  * S [uint8_T n x 1] vector containing numbers 0 and 1.
36  * x [float64_t n x 1] solution vector; must be feasible.
37  * n [uint32_t 1 x 1] dimension of H.
38  * MaxIter [uint32_t 1 x 1] max number of iterations.
39  * TolAbs [float64_t 1 x 1] Absolute tolerance.
40  * TolRel [float64_t 1 x 1] Relative tolerance.
41  * QP_TH [float64_t 1 x 1] Threshold on the primal value.
42  * print_state print function; if == NULL it is not called.
43  *
44  * RETURN VALUE
45  * structure [libqp_state_T]
46  * .QP [1 x 1] Primal objective value.
47  * .QD [1 x 1] Dual objective value.
48  * .nIter [1 x 1] Number of iterations.
49  * .exitflag [1 x 1] Indicates which stopping condition was used:
50  * -1 ... Not enough memory.
51  * 0 ... Maximal number of iteations reached: nIter >= MaxIter.
52  * 1 ... Relarive tolerance reached: QP-QD <= abs(QP)*TolRel
53  * 2 ... Absolute tolerance reached: QP-QD <= TolAbs
54  * 3 ... Objective value reached threshold: QP <= QP_TH.
55  *
56  * REFERENCE
57  * The algorithm is described in:
58  * V. Franc, V. Hlavac. A Novel Algorithm for Learning Support Vector Machines
59  * with Structured Output Spaces. Research Report K333 22/06, CTU-CMP-2006-04.
60  * May, 2006. ftp://cmp.felk.cvut.cz/pub/cmp/articles/franc/Franc-TR-2006-04.ps
61  *
62  * Copyright (C) 2006-2008 Vojtech Franc, xfrancv@cmp.felk.cvut.cz
63  * Center for Machine Perception, CTU FEL Prague
64  *
65  * This program is free software; you can redistribute it and/or
66  * modify it under the terms of the GNU General Public
68  * Version 3, 29 June 2007
69  *-------------------------------------------------------------------- */
70
71 #include <math.h>
72 #include <stdlib.h>
73 #include <stdio.h>
74 #include <string.h>
75 #include <stdint.h>
76 #include <limits.h>
77
78 #include <shogun/lib/common.h>
80 namespace shogun
81 {
82
83 libqp_state_T libqp_splx_solver(const float64_t* (*get_col)(uint32_t),
84  float64_t *diag_H,
85  float64_t *f,
86  float64_t *b,
87  uint32_t *I,
88  uint8_t *S,
89  float64_t *x,
90  uint32_t n,
91  uint32_t MaxIter,
92  float64_t TolAbs,
93  float64_t TolRel,
94  float64_t QP_TH,
95  void (*print_state)(libqp_state_T state))
96 {
97  float64_t *d;
98  float64_t *col_u, *col_v;
99  float64_t *x_neq;
100  float64_t tmp;
101  float64_t improv;
102  float64_t tmp_num;
103  float64_t tmp_den=0;
104  float64_t tau=0;
106  uint32_t *inx;
107  uint32_t *nk;
108  uint32_t m;
109  int32_t u=0;
110  int32_t v=0;
111  uint32_t k;
112  uint32_t i, j;
113  libqp_state_T state;
114
115
116  /* ------------------------------------------------------------
117  Initialization
118  ------------------------------------------------------------ */
119  state.nIter = 0;
120  state.QP = LIBQP_PLUS_INF;
121  state.QD = -LIBQP_PLUS_INF;
122  state.exitflag = 100;
123
124  inx=NULL;
125  nk=NULL;
126  d=NULL;
127  x_neq = NULL;
128
129  /* count number of constraints */
130  for( i=0, m=0; i < n; i++ )
131  m = LIBQP_MAX(m,I[i]);
132
133  /* auxciliary variables for tranforming equalities to inequalities */
134  x_neq = (float64_t*) LIBQP_CALLOC(m, sizeof(float64_t));
135  if( x_neq == NULL )
136  {
137  state.exitflag=-1;
138  goto cleanup;
139  }
140
141  /* inx is translation table between variable index i and its contraint */
142  inx = (uint32_t*) LIBQP_CALLOC(m*n, sizeof(uint32_t));
143  if( inx == NULL )
144  {
145  state.exitflag=-1;
146  goto cleanup;
147  }
148
149  /* nk is the number of variables coupled by i-th linear constraint */
150  nk = (uint32_t*) LIBQP_CALLOC(m, sizeof(uint32_t));
151  if( nk == NULL )
152  {
153  state.exitflag=-1;
154  goto cleanup;
155  }
156
157  /* setup auxciliary variables */
158  for( i=0; i < m; i++ )
159  x_neq[i] = b[i];
160
161
162  /* create inx and nk */
163  for( i=0; i < n; i++ ) {
164  k = I[i]-1;
165  inx[LIBQP_INDEX(nk[k],k,n)] = i;
166  nk[k]++;
167
168  if(S[k] != 0)
169  x_neq[k] -= x[i];
170  }
171
172  /* d = H*x + f is gradient*/
173  d = (float64_t*) LIBQP_CALLOC(n, sizeof(float64_t));
174  if( d == NULL )
175  {
176  state.exitflag=-1;
177  goto cleanup;
178  }
179
181  for( i=0; i < n; i++ )
182  {
183  d[i] += f[i];
184  if( x[i] > 0 ) {
185  col_u = (float64_t*)get_col(i);
186  for( j=0; j < n; j++ ) {
187  d[j] += col_u[j]*x[i];
188  }
189  }
190  }
191
192  /* compute state.QP = 0.5*x'*(f+d);
193  state.QD = 0.5*x'*(f-d); */
194  for( i=0, state.QP = 0, state.QD=0; i < n; i++)
195  {
196  state.QP += x[i]*(f[i]+d[i]);
197  state.QD += x[i]*(f[i]-d[i]);
198  }
199  state.QP = 0.5*state.QP;
200  state.QD = 0.5*state.QD;
201
202  for( i=0; i < m; i++ )
203  {
204  for( j=0, tmp = LIBQP_PLUS_INF; j < nk[i]; j++ )
205  tmp = LIBQP_MIN(tmp, d[inx[LIBQP_INDEX(j,i,n)]]);
206
207  if(S[i] == 0)
208  state.QD += b[i]*tmp;
209  else
210  state.QD += b[i]*LIBQP_MIN(tmp,0);
211  }
212
213  /* print initial state */
214  if( print_state != NULL)
215  print_state( state );
216
217  /* ------------------------------------------------------------
218  Main optimization loop
219  ------------------------------------------------------------ */
220  while( state.exitflag == 100 )
221  {
222  state.nIter ++;
223
224  /* go over blocks of variables coupled by lin. constraint */
225  for( k=0; k < m; k++ )
226  {
227
228  /* compute u = argmin_{i in I_k} d[i]
229  delta = sum_{i in I_k} x[i]*d[i] - b*min_{i in I_k} */
230  for( j=0, tmp = LIBQP_PLUS_INF, delta = 0; j < nk[k]; j++ )
231  {
232  i = inx[LIBQP_INDEX(j,k,n)];
233  delta += x[i]*d[i];
234  if( tmp > d[i] ) {
235  tmp = d[i];
236  u = i;
237  }
238  }
239
240  if(S[k] != 0 && d[u] > 0)
241  u = -1;
242  else
243  delta -= b[k]*d[u];
244
245  /* if satisfied then k-th block of variables needs update */
246  if( delta > TolAbs/m && delta > TolRel*LIBQP_ABS(state.QP)/m)
247  {
248  /* for fixed u select v = argmax_{i in I_k} Improvement(i) */
249  if( u != -1 )
250  {
251  col_u = (float64_t*)get_col(u);
252  improv = -LIBQP_PLUS_INF;
253  for( j=0; j < nk[k]; j++ )
254  {
255  i = inx[LIBQP_INDEX(j,k,n)];
256
257  if(x[i] > 0 && i != uint32_t(u))
258  {
259  tmp_num = x[i]*(d[i] - d[u]);
260  tmp_den = x[i]*x[i]*(diag_H[u] - 2*col_u[i] + diag_H[i]);
261  if( tmp_den > 0 )
262  {
263  if( tmp_num < tmp_den )
264  tmp = tmp_num*tmp_num / tmp_den;
265  else
266  tmp = tmp_num - 0.5 * tmp_den;
267
268  if( tmp > improv )
269  {
270  improv = tmp;
271  tau = LIBQP_MIN(1,tmp_num/tmp_den);
272  v = i;
273  }
274  }
275  }
276  }
277
278  /* check if virtual variable can be for updated */
279  if(x_neq[k] > 0 && S[k] != 0)
280  {
281  tmp_num = -x_neq[k]*d[u];
282  tmp_den = x_neq[k]*x_neq[k]*diag_H[u];
283  if( tmp_den > 0 )
284  {
285  if( tmp_num < tmp_den )
286  tmp = tmp_num*tmp_num / tmp_den;
287  else
288  tmp = tmp_num - 0.5 * tmp_den;
289
290  if( tmp > improv )
291  {
292  improv = tmp;
293  tau = LIBQP_MIN(1,tmp_num/tmp_den);
294  v = -1;
295  }
296  }
297  }
298
299  /* minimize objective w.r.t variable u and v */
300  if(v != -1)
301  {
302  tmp = x[v]*tau;
303  x[u] += tmp;
304  x[v] -= tmp;
305
306  /* update d = H*x + f */
307  col_v = (float64_t*)get_col(v);
308  for(i = 0; i < n; i++ )
309  d[i] += tmp*(col_u[i]-col_v[i]);
310  }
311  else
312  {
313  tmp = x_neq[k]*tau;
314  x[u] += tmp;
315  x_neq[k] -= tmp;
316
317  /* update d = H*x + f */
318  for(i = 0; i < n; i++ )
319  d[i] += tmp*col_u[i];
320  }
321  }
322  else
323  {
324  improv = -LIBQP_PLUS_INF;
325  for( j=0; j < nk[k]; j++ )
326  {
327  i = inx[LIBQP_INDEX(j,k,n)];
328
329  if(x[i] > 0)
330  {
331  tmp_num = x[i]*d[i];
332  tmp_den = x[i]*x[i]*diag_H[i];
333  if( tmp_den > 0 )
334  {
335  if( tmp_num < tmp_den )
336  tmp = tmp_num*tmp_num / tmp_den;
337  else
338  tmp = tmp_num - 0.5 * tmp_den;
339
340  if( tmp > improv )
341  {
342  improv = tmp;
343  tau = LIBQP_MIN(1,tmp_num/tmp_den);
344  v = i;
345  }
346  }
347  }
348  }
349
350  tmp = x[v]*tau;
351  x_neq[k] += tmp;
352  x[v] -= tmp;
353
354  /* update d = H*x + f */
355  col_v = (float64_t*)get_col(v);
356  for(i = 0; i < n; i++ )
357  d[i] -= tmp*col_v[i];
358  }
359
360  /* update objective value */
361  state.QP = state.QP - improv;
362  }
363  }
364
365  /* Compute primal and dual objectives */
366  for( i=0, state.QP = 0, state.QD=0; i < n; i++)
367  {
368  state.QP += x[i]*(f[i]+d[i]);
369  state.QD += x[i]*(f[i]-d[i]);
370  }
371  state.QP = 0.5*state.QP;
372  state.QD = 0.5*state.QD;
373
374  for( k=0; k < m; k++ )
375  {
376  for( j=0,tmp = LIBQP_PLUS_INF; j < nk[k]; j++ ) {
377  i = inx[LIBQP_INDEX(j,k,n)];
378  tmp = LIBQP_MIN(tmp, d[i]);
379  }
380
381  if(S[k] == 0)
382  state.QD += b[k]*tmp;
383  else
384  state.QD += b[k]*LIBQP_MIN(tmp,0);
385  }
386
387  /* print state */
388  if( print_state != NULL)
389  print_state( state );
390
391  /* check stopping conditions */
392  if(state.QP-state.QD <= LIBQP_ABS(state.QP)*TolRel ) state.exitflag = 1;
393  else if( state.QP-state.QD <= TolAbs ) state.exitflag = 2;
394  else if( state.QP <= QP_TH ) state.exitflag = 3;
395  else if( state.nIter >= MaxIter) state.exitflag = 0;
396  }
397
398  /*----------------------------------------------------------
399  Clean up
400  ---------------------------------------------------------- */
401 cleanup:
402  LIBQP_FREE( d );
403  LIBQP_FREE( inx );
404  LIBQP_FREE( nk );
405  LIBQP_FREE( x_neq );
406
407  return( state );
408 }
409 }
410

SHOGUN Machine Learning Toolbox - Documentation