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

SHOGUN Machine Learning Toolbox - Documentation