SHOGUN  4.2.0
GMNPLib.cpp
Go to the documentation of this file.
1 /*-----------------------------------------------------------------------
2  *
3  * This program is free software; you can redistribute it and/or modify
4  * it under the terms of the GNU General Public License as published by
5  * the Free Software Foundation; either version 3 of the License, or
6  * (at your option) any later version.
7  *
8  * Library of solvers for Generalized Nearest Point Problem (GNPP).
9  *
10  * Written (W) 1999-2008 Vojtech Franc, xfrancv@cmp.felk.cvut.cz
11  * Copyright (C) 1999-2008 Center for Machine Perception, CTU FEL Prague
12  *
13  *
14 gmnplib.c: Library of solvers for Generalized Minimal Norm Problem (GMNP).
15
16  Generalized Minimal Norm Problem to solve is
17
18  min 0.5*alpha'*H*alpha + c'*alpha
19
20  subject to sum(alpha) = 1, alpha(i) >= 0
21
22  H [dim x dim] is symmetric positive definite matrix.
23  c [dim x 1] is an arbitrary vector.
24
25  The precision of the found solution is given by
26  the parameters tmax, tolabs and tolrel which
27  define the stopping conditions:
28
29  UB-LB <= tolabs -> exit_flag = 1 Abs. tolerance.
30  UB-LB <= UB*tolrel -> exit_flag = 2 Relative tolerance.
31  LB > th -> exit_flag = 3 Threshold on lower bound.
32  t >= tmax -> exit_flag = 0 Number of iterations.
33
34  UB ... Upper bound on the optimal solution.
35  LB ... Lower bound on the optimal solution.
36  t ... Number of iterations.
37  History ... Value of LB and UB wrt. number of iterations.
38
39
40  The following algorithms are implemented:
41  ..............................................
42
43  - GMNP solver based on improved MDM algorithm 1 (u fixed v optimized)
44  exitflag = gmnp_imdm( &get_col, diag_H, vector_c, dim,
45  tmax, tolabs, tolrel, th, &alpha, &t, &History, verb );
46
47  For more info refer to V.Franc: Optimization Algorithms for Kernel
48  Methods. Research report. CTU-CMP-2005-22. CTU FEL Prague. 2005.
49  ftp://cmp.felk.cvut.cz/pub/cmp/articles/franc/Franc-PhD.pdf .
50
51  Modifications:
52  09-sep-2005, VF
53  24-jan-2005, VF
54  26-nov-2004, VF
55  25-nov-2004, VF
56  21-nov-2004, VF
57  20-nov-2004, VF
58  31-may-2004, VF
59  23-Jan-2004, VF
60
61 -------------------------------------------------------------------- */
62
65
66 #include <string.h>
67 #include <limits.h>
68
69 using namespace shogun;
70
71 #define HISTORY_BUF 1000000
72
73 #define MINUS_INF INT_MIN
74 #define PLUS_INF INT_MAX
75
76 #define INDEX(ROW,COL,DIM) ((COL*DIM)+ROW)
77 #define KDELTA(A,B) (A==B)
78 #define KDELTA4(A1,A2,A3,A4) ((A1==A2)||(A1==A3)||(A1==A4)||(A2==A3)||(A2==A4)||(A3==A4))
79
81 {
82  SG_UNSTABLE("CGMNPLib::CGMNPLib()", "\n")
83
84  diag_H = NULL;
85  kernel_columns = NULL;
86  cache_index = NULL;
87  first_kernel_inx = 0;
88  Cache_Size = 0;
89  m_num_data = 0;
90  m_reg_const = 0;
91  m_vector_y = 0;
92  m_kernel = NULL;
93
94  first_virt_inx = 0;
95  memset(&virt_columns, 0, sizeof (virt_columns));
96  m_num_virt_data = 0;
97  m_num_classes = 0;
98 }
99
101  float64_t* vector_y, CKernel* kernel, int32_t num_data,
102  int32_t num_virt_data, int32_t num_classes, float64_t reg_const)
103 : CSGObject()
104 {
105  m_num_classes=num_classes;
106  m_num_virt_data=num_virt_data;
107  m_reg_const = reg_const;
108  m_num_data = num_data;
109  m_vector_y = vector_y;
110  m_kernel = kernel;
111
112  Cache_Size = ((int64_t) kernel->get_cache_size())*1024*1024/(sizeof(float64_t)*num_data);
113  Cache_Size = CMath::min(Cache_Size, (int64_t) num_data);
114
115  SG_INFO("using %d kernel cache lines\n", Cache_Size)
116  ASSERT(Cache_Size>=2)
117
118  /* allocates memory for kernel cache */
119  kernel_columns = SG_MALLOC(float64_t*, Cache_Size);
120  cache_index = SG_MALLOC(float64_t, Cache_Size);
121
122  for(int32_t i = 0; i < Cache_Size; i++ )
123  {
124  kernel_columns[i] = SG_MALLOC(float64_t, num_data);
125  cache_index[i] = -2;
126  }
127  first_kernel_inx = 0;
128
129
130
131  for(int32_t i = 0; i < 3; i++ )
132  {
133  virt_columns[i] = SG_MALLOC(float64_t, num_virt_data);
134  }
135  first_virt_inx = 0;
136
137  diag_H = SG_MALLOC(float64_t, num_virt_data);
138
139  for(int32_t i = 0; i < num_virt_data; i++ )
140  diag_H[i] = kernel_fce(i,i);
141 }
142
144 {
145  for(int32_t i = 0; i < Cache_Size; i++ )
146  SG_FREE(kernel_columns[i]);
147
148  for(int32_t i = 0; i < 3; i++ )
149  SG_FREE(virt_columns[i]);
150
151  SG_FREE(cache_index);
152  SG_FREE(kernel_columns);
153
154  SG_FREE(diag_H);
155 }
156
157 /* ------------------------------------------------------------
158  Returns pointer at a-th column of the kernel matrix.
159  This function maintains FIFO cache of kernel columns.
160 ------------------------------------------------------------ */
162 {
163  float64_t *col_ptr;
164  int32_t i;
165  int32_t inx;
166
167  inx = -1;
168  for( i=0; i < Cache_Size; i++ ) {
169  if( cache_index[i] == a ) { inx = i; break; }
170  }
171
172  if( inx != -1 ) {
173  col_ptr = kernel_columns[inx];
174  return( col_ptr );
175  }
176
177  col_ptr = kernel_columns[first_kernel_inx];
179
181  if( first_kernel_inx >= Cache_Size ) first_kernel_inx = 0;
182
183  for( i=0; i < m_num_data; i++ ) {
184  col_ptr[i] = m_kernel->kernel(i,a);
185  }
186
187  return( col_ptr );
188 }
189
190 /* ------------------------------------------------------------
191  Computes index of input example and its class label from
192  index of virtual "single-class" example.
193 ------------------------------------------------------------ */
194 void CGMNPLib::get_indices2( int32_t *index, int32_t *c, int32_t i )
195 {
196  *index = i / (m_num_classes-1);
197
198  *c= (i % (m_num_classes-1))+1;
199  if( *c>= m_vector_y[ *index ]) (*c)++;
200
201  return;
202 }
203
204
205 /* ------------------------------------------------------------
206  Returns pointer at the a-th column of the virtual K matrix.
207
208  (note: the b-th column must be preserved in the cache during
209  updating but b is from (a(t-2), a(t-1)) where a=a(t) and
210  thus FIFO with three columns does not have to take care od b.)
211 ------------------------------------------------------------ */
212 float64_t* CGMNPLib::get_col( int32_t a, int32_t b )
213 {
214  int32_t i;
215  float64_t *col_ptr;
216  float64_t *ker_ptr;
217  float64_t value;
218  int32_t i1,c1,i2,c2;
219
220  col_ptr = virt_columns[first_virt_inx++];
221  if( first_virt_inx >= 3 ) first_virt_inx = 0;
222
223  get_indices2( &i1, &c1, a );
224  ker_ptr = (float64_t*) get_kernel_col( i1 );
225
226  for( i=0; i < m_num_virt_data; i++ ) {
227  get_indices2( &i2, &c2, i );
228
229  if( KDELTA4(m_vector_y[i1],m_vector_y[i2],c1,c2) ) {
230  value = (+KDELTA(m_vector_y[i1],m_vector_y[i2])
231  -KDELTA(m_vector_y[i1],c2)
232  -KDELTA(m_vector_y[i2],c1)
233  +KDELTA(c1,c2)
234  )*(ker_ptr[i2]+1);
235  }
236  else
237  {
238  value = 0;
239  }
240
241  if(a==i) value += m_reg_const;
242
243  col_ptr[i] = value;
244  }
245
246  return( col_ptr );
247 }
248
249
250 /* --------------------------------------------------------------
251  GMNP solver based on improved MDM algorithm 1.
252
253  Search strategy: u determined by common rule and v is
254  optimized.
255
256  Usage: exitflag = gmnp_imdm( &get_col, diag_H, vector_c, dim,
257  tmax, tolabs, tolrel, th, &alpha, &t, &History );
258 -------------------------------------------------------------- */
259
261  int32_t dim,
262  int32_t tmax,
263  float64_t tolabs,
264  float64_t tolrel,
265  float64_t th,
266  float64_t *alpha,
267  int32_t *ptr_t,
268  float64_t **ptr_History,
269  int32_t verb)
270 {
271  float64_t LB;
272  float64_t UB;
273  float64_t aHa, ac;
274  float64_t tmp, tmp1;
275  float64_t Huu, Huv, Hvv;
276  float64_t min_beta, beta;
277  float64_t max_improv, improv;
278  float64_t lambda;
279  float64_t *History;
280  float64_t *Ha;
281  float64_t *tmp_ptr;
282  float64_t *col_u, *col_v;
283  int32_t u=0;
284  int32_t v=0;
285  int32_t new_u=0;
286  int32_t i;
287  int32_t t;
288  int32_t History_size;
289  int8_t exitflag;
290
291  /* ------------------------------------------------------------ */
292  /* Initialization */
293  /* ------------------------------------------------------------ */
294
295  Ha = SG_MALLOC(float64_t, dim);
296  if( Ha == NULL ) SG_ERROR("Not enough memory.")
297
298  History_size = (tmax < HISTORY_BUF ) ? tmax+1 : HISTORY_BUF;
299  History = SG_MALLOC(float64_t, History_size*2);
300  if( History == NULL ) SG_ERROR("Not enough memory.")
301
302  /* inx = argmin(0.5*diag_H + vector_c ); */
303  for( tmp1 = PLUS_INF, i = 0; i < dim; i++ ) {
304  tmp = 0.5*diag_H[i] + vector_c[i];
305  if( tmp1 > tmp) {
306  tmp1 = tmp;
307  v = i;
308  }
309  }
310
311  col_v = (float64_t*)get_col(v,-1);
312
313  for( min_beta = PLUS_INF, i = 0; i < dim; i++ )
314  {
315  alpha[i] = 0;
316  Ha[i] = col_v[i];
317
318  beta = Ha[i] + vector_c[i];
319  if( beta < min_beta ) {
320  min_beta = beta;
321  u = i;
322  }
323  }
324
325  alpha[v] = 1;
326  aHa = diag_H[v];
327  ac = vector_c[v];
328
329  UB = 0.5*aHa + ac;
330  LB = min_beta - 0.5*aHa;
331
332  t = 0;
333  History[INDEX(0,0,2)] = LB;
334  History[INDEX(1,0,2)] = UB;
335
336  if( verb ) {
337  SG_PRINT("Init: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
338  UB, LB, UB-LB,(UB-LB)/UB);
339  }
340
341  /* Stopping conditions */
342  if( UB-LB <= tolabs ) exitflag = 1;
343  else if(UB-LB <= CMath::abs(UB)*tolrel ) exitflag = 2;
344  else if(LB > th ) exitflag = 3;
345  else exitflag = -1;
346
347  /* ------------------------------------------------------------ */
348  /* Main optimization loop */
349  /* ------------------------------------------------------------ */
350
351  col_u = (float64_t*)get_col(u,-1);
352  while( exitflag == -1 )
353  {
354  t++;
355
356  col_v = (float64_t*)get_col(v,u);
357
358  /* Adaptation rule and update */
359  Huu = diag_H[u];
360  Hvv = diag_H[v];
361  Huv = col_u[v];
362
363  lambda = (Ha[v]-Ha[u]+vector_c[v]-vector_c[u])/(alpha[v]*(Huu-2*Huv+Hvv));
364  if( lambda < 0 ) lambda = 0; else if (lambda > 1) lambda = 1;
365
366  aHa = aHa + 2*alpha[v]*lambda*(Ha[u]-Ha[v])+
367  lambda*lambda*alpha[v]*alpha[v]*(Huu-2*Huv+Hvv);
368
369  ac = ac + lambda*alpha[v]*(vector_c[u]-vector_c[v]);
370
371  tmp = alpha[v];
372  alpha[u]=alpha[u]+lambda*alpha[v];
373  alpha[v]=alpha[v]-lambda*alpha[v];
374
375  UB = 0.5*aHa + ac;
376
377 /* max_beta = MINUS_INF;*/
378  for( min_beta = PLUS_INF, i = 0; i < dim; i++ )
379  {
380  Ha[i] = Ha[i] + lambda*tmp*(col_u[i] - col_v[i]);
381
382  beta = Ha[i]+ vector_c[i];
383
384  if( beta < min_beta )
385  {
386  new_u = i;
387  min_beta = beta;
388  }
389  }
390
391  LB = min_beta - 0.5*aHa;
392  u = new_u;
393  col_u = (float64_t*)get_col(u,-1);
394
395  /* search for optimal v while u is fixed */
396  for( max_improv = MINUS_INF, i = 0; i < dim; i++ ) {
397
398  if( alpha[i] != 0 ) {
399  beta = Ha[i] + vector_c[i];
400
401  if( beta >= min_beta ) {
402
403  tmp = diag_H[u] - 2*col_u[i] + diag_H[i];
404  if( tmp != 0 ) {
405  improv = (0.5*(beta-min_beta)*(beta-min_beta))/tmp;
406
407  if( improv > max_improv ) {
408  max_improv = improv;
409  v = i;
410  }
411  }
412  }
413  }
414  }
415
416  /* Stopping conditions */
417  if( UB-LB <= tolabs ) exitflag = 1;
418  else if( UB-LB <= CMath::abs(UB)*tolrel) exitflag = 2;
419  else if(LB > th ) exitflag = 3;
420  else if(t >= tmax) exitflag = 0;
421
422  /* print info */
423  SG_ABS_PROGRESS(CMath::abs((UB-LB)/UB),
424  -CMath::log10(CMath::abs(UB-LB)),
425  -CMath::log10(1.0),
426  -CMath::log10(tolrel), 6);
427  if(verb && (t % verb) == 0 ) {
428  SG_PRINT("%d: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
429  t, UB, LB, UB-LB,(UB-LB)/UB);
430  }
431
432  /* Store selected values */
433  if( t < History_size ) {
434  History[INDEX(0,t,2)] = LB;
435  History[INDEX(1,t,2)] = UB;
436  }
437  else {
438  tmp_ptr = SG_MALLOC(float64_t, (History_size+HISTORY_BUF)*2);
439  if( tmp_ptr == NULL ) SG_ERROR("Not enough memory.")
440  for( i = 0; i < History_size; i++ ) {
441  tmp_ptr[INDEX(0,i,2)] = History[INDEX(0,i,2)];
442  tmp_ptr[INDEX(1,i,2)] = History[INDEX(1,i,2)];
443  }
444  tmp_ptr[INDEX(0,t,2)] = LB;
445  tmp_ptr[INDEX(1,t,2)] = UB;
446
447  History_size += HISTORY_BUF;
448  SG_FREE(History);
449  History = tmp_ptr;
450  }
451  }
452
453  /* print info about last iteration*/
454  SG_DONE()
455  if(verb && (t % verb) ) {
456  SG_PRINT("exit: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
457  UB, LB, UB-LB,(UB-LB)/UB);
458  }
459
460
461  /*------------------------------------------------------- */
462  /* Set outputs */
463  /*------------------------------------------------------- */
464  (*ptr_t) = t;
465  (*ptr_History) = History;
466
467  /* Free memory */
468  SG_FREE(Ha);
469
470  return( exitflag );
471 }
472
473 /* ------------------------------------------------------------
474  Retures (a,b)-th element of the virtual kernel matrix
475  of size [num_virt_data x num_virt_data].
476 ------------------------------------------------------------ */
477 float64_t CGMNPLib::kernel_fce( int32_t a, int32_t b )
478 {
479  float64_t value;
480  int32_t i1,c1,i2,c2;
481
482  get_indices2( &i1, &c1, a );
483  get_indices2( &i2, &c2, b );
484
485  if( KDELTA4(m_vector_y[i1],m_vector_y[i2],c1,c2) ) {
486  value = (+KDELTA(m_vector_y[i1],m_vector_y[i2])
487  -KDELTA(m_vector_y[i1],c2)
488  -KDELTA(m_vector_y[i2],c1)
489  +KDELTA(c1,c2)
490  )*(m_kernel->kernel( i1, i2 )+1);
491  }
492  else
493  {
494  value = 0;
495  }
496
497  if(a==b) value += m_reg_const;
498
499  return( value );
500 }
#define SG_INFO(...)
Definition: SGIO.h:118
#define SG_DONE()
Definition: SGIO.h:157
int32_t first_kernel_inx
Definition: GMNPLib.h:148
float64_t * virt_columns[3]
Definition: GMNPLib.h:163
int32_t first_virt_inx
Definition: GMNPLib.h:161
float64_t * m_vector_y
Definition: GMNPLib.h:156
#define PLUS_INF
Definition: GMNPLib.cpp:74
float64_t m_reg_const
Definition: GMNPLib.h:154
static float64_t log10(float64_t v)
Definition: Math.h:897
CKernel * m_kernel
Definition: GMNPLib.h:158
#define SG_ERROR(...)
Definition: SGIO.h:129
float64_t ** kernel_columns
Definition: GMNPLib.h:144
float64_t kernel(int32_t idx_a, int32_t idx_b)
Definition: Kernel.h:207
float64_t * get_col(int32_t a, int32_t b)
Definition: GMNPLib.cpp:212
int32_t m_num_classes
Definition: GMNPLib.h:167
float64_t * cache_index
Definition: GMNPLib.h:146
float64_t * diag_H
Definition: GMNPLib.h:142
virtual ~CGMNPLib()
Definition: GMNPLib.cpp:143
void get_indices2(int32_t *index, int32_t *c, int32_t i)
Definition: GMNPLib.cpp:194
#define HISTORY_BUF
Definition: GMNPLib.cpp:71
#define SG_PRINT(...)
Definition: SGIO.h:137
int32_t m_num_data
Definition: GMNPLib.h:152
#define ASSERT(x)
Definition: SGIO.h:201
float64_t kernel_fce(int32_t a, int32_t b)
Definition: GMNPLib.cpp:477
Class SGObject is the base class of all shogun objects.
Definition: SGObject.h:115
double float64_t
Definition: common.h:50
int8_t gmnp_imdm(float64_t *vector_c, int32_t dim, int32_t tmax, float64_t tolabs, float64_t tolrel, float64_t th, float64_t *alpha, int32_t *ptr_t, float64_t **ptr_History, int32_t verb)
Definition: GMNPLib.cpp:260
#define INDEX(ROW, COL, DIM)
Definition: GMNPLib.cpp:76
#define MINUS_INF
Definition: GMNPLib.cpp:73
#define SG_ABS_PROGRESS(...)
Definition: SGIO.h:152
int32_t m_num_virt_data
Definition: GMNPLib.h:165
#define KDELTA(A, B)
Definition: GMNPLib.cpp:77
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
static T min(T a, T b)
Definition: Math.h:157
The Kernel base class.
Definition: Kernel.h:159
int32_t get_cache_size()
Definition: Kernel.h:599
float64_t * get_kernel_col(int32_t a)
Definition: GMNPLib.cpp:161
int64_t Cache_Size
Definition: GMNPLib.h:150
#define SG_UNSTABLE(func,...)
Definition: SGIO.h:132
#define KDELTA4(A1, A2, A3, A4)
Definition: GMNPLib.cpp:78
static T abs(T a)
Definition: Math.h:179

SHOGUN Machine Learning Toolbox - Documentation