SHOGUN  v3.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
libbmrm.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  * libbmrm.h: Implementation of the BMRM solver for SO training
8  *
9  * Copyright (C) 2012 Michal Uricar, uricamic@cmp.felk.cvut.cz
10  *
11  * Implementation of the BMRM solver
12  *--------------------------------------------------------------------- */
13 
15 #include <shogun/lib/external/libqp.h>
16 #include <shogun/lib/Time.h>
17 
18 namespace shogun
19 {
20 static const uint32_t QPSolverMaxIter=0xFFFFFFFF;
21 static const float64_t epsilon=0.0;
22 
23 static float64_t *H;
24 static uint32_t BufSize;
25 
27  bmrm_ll** tail,
28  bool* map,
29  float64_t* A,
30  uint32_t free_idx,
31  float64_t* cp_data,
32  uint32_t dim)
33 {
34  ASSERT(map[free_idx])
35 
36  LIBBMRM_MEMCPY(A+free_idx*dim, cp_data, dim*sizeof(float64_t));
37  map[free_idx]=false;
38 
40 
41  if (cp==NULL)
42  {
43  SG_SERROR("Out of memory.\n")
44  return;
45  }
46 
47  cp->address=A+(free_idx*dim);
48  cp->prev=*tail;
49  cp->next=NULL;
50  cp->idx=free_idx;
51  (*tail)->next=cp;
52  *tail=cp;
53 }
54 
56  bmrm_ll** head,
57  bmrm_ll** tail,
58  bool* map,
59  float64_t* icp)
60 {
61  bmrm_ll *cp_list_ptr=*head;
62 
63  while(cp_list_ptr->address != icp)
64  {
65  cp_list_ptr=cp_list_ptr->next;
66  }
67 
68  if (cp_list_ptr==*head)
69  {
70  *head=(*head)->next;
71  cp_list_ptr->next->prev=NULL;
72  }
73  else if (cp_list_ptr==*tail)
74  {
75  *tail=(*tail)->prev;
76  cp_list_ptr->prev->next=NULL;
77  }
78  else
79  {
80  cp_list_ptr->prev->next=cp_list_ptr->next;
81  cp_list_ptr->next->prev=cp_list_ptr->prev;
82  }
83 
84  map[cp_list_ptr->idx]=true;
85  LIBBMRM_FREE(cp_list_ptr);
86 }
87 
88 void clean_icp(ICP_stats* icp_stats,
89  BmrmStatistics& bmrm,
90  bmrm_ll** head,
91  bmrm_ll** tail,
92  float64_t*& Hmat,
93  float64_t*& diag_H,
94  float64_t*& beta,
95  bool*& map,
96  uint32_t cleanAfter,
97  float64_t*& b,
98  uint32_t*& I,
99  uint32_t cp_models
100  )
101 {
102  /* find ICP */
103  uint32_t cntICP=0;
104  uint32_t cntACP=0;
105  bmrm_ll* cp_ptr=*head;
106  uint32_t tmp_idx=0;
107 
108  while (cp_ptr != *tail)
109  {
110  if (icp_stats->ICPcounter[tmp_idx++]>=cleanAfter)
111  {
112  icp_stats->ICPs[cntICP++]=cp_ptr->address;
113  }
114  else
115  {
116  icp_stats->ACPs[cntACP++]=tmp_idx-1;
117  }
118 
119  cp_ptr=cp_ptr->next;
120  }
121 
122  /* do ICP removal */
123  if (cntICP > 0)
124  {
125  uint32_t nCP_new=bmrm.nCP-cntICP;
126 
127  for (uint32_t i=0; i<cntICP; ++i)
128  {
129  tmp_idx=0;
130  cp_ptr=*head;
131 
132  while(cp_ptr->address != icp_stats->ICPs[i])
133  {
134  cp_ptr=cp_ptr->next;
135  tmp_idx++;
136  }
137 
138  remove_cutting_plane(head, tail, map, icp_stats->ICPs[i]);
139 
140  LIBBMRM_MEMMOVE(b+tmp_idx, b+tmp_idx+1,
141  (bmrm.nCP+cp_models-tmp_idx)*sizeof(float64_t));
142  LIBBMRM_MEMMOVE(beta+tmp_idx, beta+tmp_idx+1,
143  (bmrm.nCP-tmp_idx)*sizeof(float64_t));
144  LIBBMRM_MEMMOVE(diag_H+tmp_idx, diag_H+tmp_idx+1,
145  (bmrm.nCP-tmp_idx)*sizeof(float64_t));
146  LIBBMRM_MEMMOVE(I+tmp_idx, I+tmp_idx+1,
147  (bmrm.nCP-tmp_idx)*sizeof(uint32_t));
148  LIBBMRM_MEMMOVE(icp_stats->ICPcounter+tmp_idx, icp_stats->ICPcounter+tmp_idx+1,
149  (bmrm.nCP-tmp_idx)*sizeof(uint32_t));
150  }
151 
152  /* H */
153  for (uint32_t i=0; i < nCP_new; ++i)
154  {
155  for (uint32_t j=0; j < nCP_new; ++j)
156  {
157  icp_stats->H_buff[LIBBMRM_INDEX(i, j, icp_stats->maxCPs)]=
158  Hmat[LIBBMRM_INDEX(icp_stats->ACPs[i], icp_stats->ACPs[j], icp_stats->maxCPs)];
159  }
160  }
161 
162  for (uint32_t i=0; i<nCP_new; ++i)
163  for (uint32_t j=0; j<nCP_new; ++j)
164  Hmat[LIBBMRM_INDEX(i, j, icp_stats->maxCPs)]=
165  icp_stats->H_buff[LIBBMRM_INDEX(i, j, icp_stats->maxCPs)];
166 
167  bmrm.nCP=nCP_new;
168  }
169 }
170 
171 /*----------------------------------------------------------------------
172  Returns pointer at i-th column of Hessian matrix.
173  ----------------------------------------------------------------------*/
174 static const float64_t *get_col( uint32_t i)
175 {
176  return( &H[ BufSize*i ] );
177 }
178 
180  CDualLibQPBMSOSVM *machine,
181  float64_t* W,
182  float64_t TolRel,
183  float64_t TolAbs,
184  float64_t _lambda,
185  uint32_t _BufSize,
186  bool cleanICP,
187  uint32_t cleanAfter,
188  float64_t K,
189  uint32_t Tmax,
190  bool verbose)
191 {
192  BmrmStatistics bmrm;
193  libqp_state_T qp_exitflag={0, 0, 0, 0};
194  float64_t *b, *beta, *diag_H, *prevW;
195  float64_t R, *subgrad, *A, QPSolverTolRel, C=1.0, wdist=0.0;
196  floatmax_t rsum, sq_norm_W, sq_norm_Wdiff=0.0;
197  uint32_t *I;
198  uint8_t S=1;
199  CStructuredModel* model=machine->get_model();
200  uint32_t nDim=model->get_dim();
201  SG_UNREF(model);
202 
203  CTime ttime;
204  float64_t tstart, tstop;
205 
206  bmrm_ll *CPList_head, *CPList_tail, *cp_ptr, *cp_ptr2, *cp_list=NULL;
207  float64_t *A_1=NULL, *A_2=NULL;
208  bool *map=NULL;
209 
210 
211  tstart=ttime.cur_time_diff(false);
212 
213  BufSize=_BufSize;
214  QPSolverTolRel=1e-9;
215 
216  H=NULL;
217  b=NULL;
218  beta=NULL;
219  A=NULL;
220  subgrad=NULL;
221  diag_H=NULL;
222  I=NULL;
223  prevW=NULL;
224 
226 
227  if (H==NULL)
228  {
229  bmrm.exitflag=-2;
230  goto cleanup;
231  }
232 
233  A= (float64_t*) LIBBMRM_CALLOC(nDim*BufSize, float64_t);
234 
235  if (A==NULL)
236  {
237  bmrm.exitflag=-2;
238  goto cleanup;
239  }
240 
241  b= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
242 
243  if (b==NULL)
244  {
245  bmrm.exitflag=-2;
246  goto cleanup;
247  }
248 
249  beta= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
250 
251  if (beta==NULL)
252  {
253  bmrm.exitflag=-2;
254  goto cleanup;
255  }
256 
257  subgrad= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t);
258 
259  if (subgrad==NULL)
260  {
261  bmrm.exitflag=-2;
262  goto cleanup;
263  }
264 
265  diag_H= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
266 
267  if (diag_H==NULL)
268  {
269  bmrm.exitflag=-2;
270  goto cleanup;
271  }
272 
273  I= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
274 
275  if (I==NULL)
276  {
277  bmrm.exitflag=-2;
278  goto cleanup;
279  }
280 
281  ICP_stats icp_stats;
282  icp_stats.maxCPs = BufSize;
283 
284  icp_stats.ICPcounter= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
285  if (icp_stats.ICPcounter==NULL)
286  {
287  bmrm.exitflag=-2;
288  goto cleanup;
289  }
290 
291  icp_stats.ICPs= (float64_t**) LIBBMRM_CALLOC(BufSize, float64_t*);
292  if (icp_stats.ICPs==NULL)
293  {
294  bmrm.exitflag=-2;
295  goto cleanup;
296  }
297 
298  icp_stats.ACPs= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
299  if (icp_stats.ACPs==NULL)
300  {
301  bmrm.exitflag=-2;
302  goto cleanup;
303  }
304 
305  /* Temporary buffers for ICP removal */
306  icp_stats.H_buff= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, float64_t);
307  if (icp_stats.H_buff==NULL)
308  {
309  bmrm.exitflag=-2;
310  goto cleanup;
311  }
312 
313  map= (bool*) LIBBMRM_CALLOC(BufSize, bool);
314 
315  if (map==NULL)
316  {
317  bmrm.exitflag=-2;
318  goto cleanup;
319  }
320 
321  memset( (bool*) map, true, BufSize);
322 
323  cp_list= (bmrm_ll*) LIBBMRM_CALLOC(1, bmrm_ll);
324 
325  if (cp_list==NULL)
326  {
327  bmrm.exitflag=-2;
328  goto cleanup;
329  }
330 
331  prevW= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t);
332 
333  if (prevW==NULL)
334  {
335  bmrm.exitflag=-2;
336  goto cleanup;
337  }
338 
342 
343  /* Iinitial solution */
344  R=machine->risk(subgrad, W);
345 
346  bmrm.nCP=0;
347  bmrm.nIter=0;
348  bmrm.exitflag=0;
349 
350  b[0]=-R;
351 
352  /* Cutting plane auxiliary double linked list */
353 
354  LIBBMRM_MEMCPY(A, subgrad, nDim*sizeof(float64_t));
355  map[0]=false;
356  cp_list->address=&A[0];
357  cp_list->idx=0;
358  cp_list->prev=NULL;
359  cp_list->next=NULL;
360  CPList_head=cp_list;
361  CPList_tail=cp_list;
362 
363  /* Compute initial value of Fp, Fd, assuming that W is zero vector */
364 
365  sq_norm_W=0;
366  bmrm.Fp=R+0.5*_lambda*sq_norm_W;
367  bmrm.Fd=-LIBBMRM_PLUS_INF;
368 
369  tstop=ttime.cur_time_diff(false);
370 
371  /* Verbose output */
372 
373  if (verbose)
374  SG_SPRINT("%4d: tim=%.3lf, Fp=%lf, Fd=%lf, R=%lf\n",
375  bmrm.nIter, tstop-tstart, bmrm.Fp, bmrm.Fd, R);
376 
377  /* store Fp, Fd and wdist history */
378  bmrm.hist_Fp[0]=bmrm.Fp;
379  bmrm.hist_Fd[0]=bmrm.Fd;
380  bmrm.hist_wdist[0]=0.0;
381 
382  /* main loop */
383 
384  while (bmrm.exitflag==0)
385  {
386  tstart=ttime.cur_time_diff(false);
387  bmrm.nIter++;
388 
389  /* Update H */
390 
391  if (bmrm.nCP>0)
392  {
393  A_2=get_cutting_plane(CPList_tail);
394  cp_ptr=CPList_head;
395 
396  for (uint32_t i=0; i<bmrm.nCP; ++i)
397  {
398  A_1=get_cutting_plane(cp_ptr);
399  cp_ptr=cp_ptr->next;
400  rsum= SGVector<float64_t>::dot(A_1, A_2, nDim);
401 
402  H[LIBBMRM_INDEX(bmrm.nCP, i, BufSize)]
403  = H[LIBBMRM_INDEX(i, bmrm.nCP, BufSize)]
404  = rsum/_lambda;
405  }
406  }
407 
408  A_2=get_cutting_plane(CPList_tail);
409  rsum = SGVector<float64_t>::dot(A_2, A_2, nDim);
410 
411  H[LIBBMRM_INDEX(bmrm.nCP, bmrm.nCP, BufSize)]=rsum/_lambda;
412 
413  diag_H[bmrm.nCP]=H[LIBBMRM_INDEX(bmrm.nCP, bmrm.nCP, BufSize)];
414  I[bmrm.nCP]=1;
415 
416  bmrm.nCP++;
417  beta[bmrm.nCP]=0.0; // [beta; 0]
418 
419 #if 0
420  /* TODO: scaling...*/
421  float64_t scale = SGVector<float64_t>::max(diag_H, BufSize)/(1000.0*_lambda);
422  SGVector<float64_t> sb(bmrm.nCP);
423  sb.zero();
424  sb.vec1_plus_scalar_times_vec2(sb.vector, 1/scale, b, bmrm.nCP);
425 
426  SGVector<float64_t> sh(bmrm.nCP);
427  sh.zero();
428  sb.vec1_plus_scalar_times_vec2(sh.vector, 1/scale, diag_H, bmrm.nCP);
429 
430  qp_exitflag =
431  libqp_splx_solver(&get_col, sh.vector, sb.vector, &C, I, &S, beta,
432  bmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
433 #else
434  /* call QP solver */
435  qp_exitflag=libqp_splx_solver(&get_col, diag_H, b, &C, I, &S, beta,
436  bmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
437 #endif
438 
439  bmrm.qp_exitflag=qp_exitflag.exitflag;
440 
441  /* Update ICPcounter (add one to unused and reset used)
442  * + compute number of active CPs */
443  bmrm.nzA=0;
444 
445  for (uint32_t aaa=0; aaa<bmrm.nCP; ++aaa)
446  {
447  if (beta[aaa]>epsilon)
448  {
449  ++bmrm.nzA;
450  icp_stats.ICPcounter[aaa]=0;
451  }
452  else
453  {
454  icp_stats.ICPcounter[aaa]+=1;
455  }
456  }
457 
458  /* W update */
459  memset(W, 0, sizeof(float64_t)*nDim);
460  cp_ptr=CPList_head;
461  for (uint32_t j=0; j<bmrm.nCP; ++j)
462  {
463  A_1=get_cutting_plane(cp_ptr);
464  cp_ptr=cp_ptr->next;
465  SGVector<float64_t>::vec1_plus_scalar_times_vec2(W, -beta[j]/_lambda, A_1, nDim);
466  }
467 
468  /* risk and subgradient computation */
469  R = machine->risk(subgrad, W);
470  add_cutting_plane(&CPList_tail, map, A,
471  find_free_idx(map, BufSize), subgrad, nDim);
472 
473  sq_norm_W=SGVector<float64_t>::dot(W, W, nDim);
474  b[bmrm.nCP]=SGVector<float64_t>::dot(subgrad, W, nDim) - R;
475 
476  sq_norm_Wdiff=0.0;
477  for (uint32_t j=0; j<nDim; ++j)
478  {
479  sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
480  }
481 
482  bmrm.Fp=R+0.5*_lambda*sq_norm_W;
483  bmrm.Fd=-qp_exitflag.QP;
484  wdist=CMath::sqrt(sq_norm_Wdiff);
485 
486  /* Stopping conditions */
487 
488  if (bmrm.Fp - bmrm.Fd <= TolRel*LIBBMRM_ABS(bmrm.Fp))
489  bmrm.exitflag=1;
490 
491  if (bmrm.Fp - bmrm.Fd <= TolAbs)
492  bmrm.exitflag=2;
493 
494  if (bmrm.nCP >= BufSize)
495  bmrm.exitflag=-1;
496 
497  tstop=ttime.cur_time_diff(false);
498 
499  /* Verbose output */
500 
501  if (verbose)
502  SG_SPRINT("%4d: tim=%.3lf, Fp=%lf, Fd=%lf, (Fp-Fd)=%lf, (Fp-Fd)/Fp=%lf, R=%lf, nCP=%d, nzA=%d, QPexitflag=%d\n",
503  bmrm.nIter, tstop-tstart, bmrm.Fp, bmrm.Fd, bmrm.Fp-bmrm.Fd,
504  (bmrm.Fp-bmrm.Fd)/bmrm.Fp, R, bmrm.nCP, bmrm.nzA, qp_exitflag.exitflag);
505 
506  /* Keep Fp, Fd and w_dist history */
507  bmrm.hist_Fp[bmrm.nIter]=bmrm.Fp;
508  bmrm.hist_Fd[bmrm.nIter]=bmrm.Fd;
509  bmrm.hist_wdist[bmrm.nIter]=wdist;
510 
511  /* Check size of Buffer */
512 
513  if (bmrm.nCP>=BufSize)
514  {
515  bmrm.exitflag=-2;
516  SG_SERROR("Buffer exceeded.\n")
517  }
518 
519  /* keep W (for wdist history track) */
520  LIBBMRM_MEMCPY(prevW, W, nDim*sizeof(float64_t));
521 
522  /* Inactive Cutting Planes (ICP) removal */
523  if (cleanICP)
524  {
525  clean_icp(&icp_stats, bmrm, &CPList_head, &CPList_tail, H, diag_H, beta, map, cleanAfter, b, I);
526  }
527  } /* end of main loop */
528 
529  bmrm.hist_Fp.resize_vector(bmrm.nIter);
530  bmrm.hist_Fd.resize_vector(bmrm.nIter);
531  bmrm.hist_wdist.resize_vector(bmrm.nIter);
532 
533  cp_ptr=CPList_head;
534 
535  while(cp_ptr!=NULL)
536  {
537  cp_ptr2=cp_ptr;
538  cp_ptr=cp_ptr->next;
539  LIBBMRM_FREE(cp_ptr2);
540  cp_ptr2=NULL;
541  }
542 
543  cp_list=NULL;
544 
545 cleanup:
546 
547  LIBBMRM_FREE(H);
548  LIBBMRM_FREE(b);
549  LIBBMRM_FREE(beta);
550  LIBBMRM_FREE(A);
551  LIBBMRM_FREE(subgrad);
552  LIBBMRM_FREE(diag_H);
553  LIBBMRM_FREE(I);
554  LIBBMRM_FREE(icp_stats.ICPcounter);
555  LIBBMRM_FREE(icp_stats.ICPs);
556  LIBBMRM_FREE(icp_stats.ACPs);
557  LIBBMRM_FREE(icp_stats.H_buff);
558  LIBBMRM_FREE(map);
559  LIBBMRM_FREE(prevW);
560 
561  if (cp_list)
562  LIBBMRM_FREE(cp_list);
563 
564  return(bmrm);
565 }
566 }

SHOGUN Machine Learning Toolbox - Documentation