SHOGUN  4.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
libppbm.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  * libp3bm.h: Implementation of the Proximal Point P-BMRM solver for SO training
8  *
9  * Copyright (C) 2012 Michal Uricar, uricamic@cmp.felk.cvut.cz
10  *
11  * Implementation of the Proximal Point P-BMRM
12  *--------------------------------------------------------------------- */
13 
14 #include <shogun/structure/libppbm.h>
15 #ifdef USE_GPL_SHOGUN
16 #include <shogun/lib/external/libqp.h>
17  #include <shogun/mathematics/Math.h>
18 #include <shogun/lib/Time.h>
19 
20 namespace shogun
21 {
22 static const uint32_t QPSolverMaxIter=0xFFFFFFFF;
23 static const float64_t epsilon=0.0;
24 
25 static float64_t *H, *H2;
26 
27 /*----------------------------------------------------------------------
28  Returns pointer at i-th column of Hessian matrix.
29  ----------------------------------------------------------------------*/
30 static const float64_t *get_col( uint32_t i)
31 {
32  return( &H2[ BufSize*i ] );
33 }
34 
35 BmrmStatistics svm_ppbm_solver(
36  CDualLibQPBMSOSVM *machine,
37  float64_t* W,
38  float64_t TolRel,
39  float64_t TolAbs,
40  float64_t _lambda,
41  uint32_t _BufSize,
42  bool cleanICP,
43  uint32_t cleanAfter,
44  float64_t K,
45  uint32_t Tmax,
46  bool verbose)
47 {
48  BmrmStatistics ppbmrm;
49  libqp_state_T qp_exitflag={0, 0, 0, 0}, qp_exitflag_good={0, 0, 0, 0};
50  float64_t *b, *b2, *beta, *beta_good, *beta_start, *diag_H, *diag_H2;
51  float64_t R, *subgrad, *A, QPSolverTolRel, C=1.0;
52  float64_t *prevW, *wt, alpha, alpha_start, alpha_good=0.0, Fd_alpha0=0.0;
53  float64_t lastFp, wdist, gamma=0.0;
54  floatmax_t rsum, sq_norm_W, sq_norm_Wdiff, sq_norm_prevW, eps;
55  uint32_t *I, *I2, *I_start, *I_good;
56  uint8_t S=1;
57  CStructuredModel* model=machine->get_model();
58  uint32_t nDim=model->get_dim();
59  CSOSVMHelper* helper = NULL;
60  uint32_t qp_cnt=0;
61  bmrm_ll *CPList_head, *CPList_tail, *cp_ptr, *cp_ptr2, *cp_list=NULL;
62  float64_t *A_1=NULL, *A_2=NULL;
63  bool *map=NULL, tuneAlpha=true, flag=true, alphaChanged=false, isThereGoodSolution=false;
64 
65  CTime ttime;
66  float64_t tstart, tstop;
67 
68 
69  tstart=ttime.cur_time_diff(false);
70 
71  BufSize=_BufSize;
72  QPSolverTolRel=1e-9;
73 
74  H=NULL;
75  b=NULL;
76  beta=NULL;
77  A=NULL;
78  subgrad=NULL;
79  diag_H=NULL;
80  I=NULL;
81  prevW=NULL;
82  wt=NULL;
83  diag_H2=NULL;
84  b2=NULL;
85  I2=NULL;
86  H2=NULL;
87  I_good=NULL;
88  I_start=NULL;
89  beta_start=NULL;
90  beta_good=NULL;
91 
92  alpha=0.0;
93 
94  H= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, float64_t);
95 
96  ASSERT(nDim > 0);
97  ASSERT(BufSize > 0);
98  REQUIRE(BufSize < (std::numeric_limits<size_t>::max() / nDim),
99  "overflow: %u * %u > %u -- biggest possible BufSize=%u or nDim=%u\n",
100  BufSize, nDim, std::numeric_limits<size_t>::max(),
101  (std::numeric_limits<size_t>::max() / nDim),
102  (std::numeric_limits<size_t>::max() / BufSize));
103 
104  A= (float64_t*) LIBBMRM_CALLOC(nDim*BufSize, float64_t);
105 
106  b= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
107 
108  beta= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
109 
110  subgrad= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t);
111 
112  diag_H= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
113 
114  I= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
115 
116  /* structure for maintaining inactive CPs info */
117  ICP_stats icp_stats;
118  icp_stats.maxCPs = BufSize;
119  icp_stats.ICPcounter= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
120  icp_stats.ICPs= (float64_t**) LIBBMRM_CALLOC(BufSize, float64_t*);
121  icp_stats.ACPs= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
122 
123  cp_list= (bmrm_ll*) LIBBMRM_CALLOC(1, bmrm_ll);
124 
125  prevW= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t);
126 
127  wt= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t);
128 
129  if (H==NULL || A==NULL || b==NULL || beta==NULL || subgrad==NULL ||
130  diag_H==NULL || I==NULL || icp_stats.ICPcounter==NULL ||
131  icp_stats.ICPs==NULL || icp_stats.ACPs==NULL ||
132  cp_list==NULL || prevW==NULL || wt==NULL)
133  {
134  ppbmrm.exitflag=-2;
135  goto cleanup;
136  }
137 
138  map= (bool*) LIBBMRM_CALLOC(BufSize, bool);
139 
140  if (map==NULL)
141  {
142  ppbmrm.exitflag=-2;
143  goto cleanup;
144  }
145 
146  memset( (bool*) map, true, BufSize);
147 
148  /* Temporary buffers for ICP removal */
149  icp_stats.H_buff= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, float64_t);
150 
151  if (icp_stats.H_buff==NULL)
152  {
153  ppbmrm.exitflag=-2;
154  goto cleanup;
155  }
156 
157  /* Temporary buffers */
158  beta_start= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
159 
160  beta_good= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
161 
162  b2= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
163 
164  diag_H2= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
165 
166  H2= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, float64_t);
167 
168  I_start= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
169 
170  I_good= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
171 
172  I2= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
173 
174  if (beta_start==NULL || beta_good==NULL || b2==NULL || diag_H2==NULL ||
175  I_start==NULL || I_good==NULL || I2==NULL || H2==NULL)
176  {
177  ppbmrm.exitflag=-2;
178  goto cleanup;
179  }
180 
181  ppbmrm.hist_Fp.resize_vector(BufSize);
182  ppbmrm.hist_Fd.resize_vector(BufSize);
183  ppbmrm.hist_wdist.resize_vector(BufSize);
184 
185  /* Iinitial solution */
186  R = machine->risk(subgrad, W);
187 
188  ppbmrm.nCP=0;
189  ppbmrm.nIter=0;
190  ppbmrm.exitflag=0;
191 
192  b[0]=-R;
193 
194  /* Cutting plane auxiliary double linked list */
195  LIBBMRM_MEMCPY(A, subgrad, nDim*sizeof(float64_t));
196  map[0]=false;
197  cp_list->address=&A[0];
198  cp_list->idx=0;
199  cp_list->prev=NULL;
200  cp_list->next=NULL;
201  CPList_head=cp_list;
202  CPList_tail=cp_list;
203 
204  /* Compute initial value of Fp, Fd, assuming that W is zero vector */
205  sq_norm_Wdiff=0.0;
206 
207  b[0] = CMath::dot(subgrad, W, nDim);
208  sq_norm_W = CMath::dot(W, W, nDim);
209  for (uint32_t j=0; j<nDim; ++j)
210  {
211  sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
212  }
213 
214  ppbmrm.Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff;
215  ppbmrm.Fd=-LIBBMRM_PLUS_INF;
216  lastFp=ppbmrm.Fp;
217  wdist=CMath::sqrt(sq_norm_Wdiff);
218 
219  K = (sq_norm_W == 0.0) ? 0.4 : 0.01*CMath::sqrt(sq_norm_W);
220 
221  LIBBMRM_MEMCPY(prevW, W, nDim*sizeof(float64_t));
222 
223  tstop=ttime.cur_time_diff(false);
224 
225  /* Keep history of Fp, Fd, wdist */
226  ppbmrm.hist_Fp[0]=ppbmrm.Fp;
227  ppbmrm.hist_Fd[0]=ppbmrm.Fd;
228  ppbmrm.hist_wdist[0]=wdist;
229 
230  /* Verbose output */
231 
232  if (verbose)
233  SG_SDEBUG("%4d: tim=%.3lf, Fp=%lf, Fd=%lf, R=%lf, K=%lf\n",
234  ppbmrm.nIter, tstop-tstart, ppbmrm.Fp, ppbmrm.Fd, R, K);
235 
236  if (verbose)
237  helper = machine->get_helper();
238 
239  /* main loop */
240 
241  while (ppbmrm.exitflag==0)
242  {
243  tstart=ttime.cur_time_diff(false);
244  ppbmrm.nIter++;
245 
246  /* Update H */
247 
248  if (ppbmrm.nCP>0)
249  {
250  A_2=get_cutting_plane(CPList_tail);
251  cp_ptr=CPList_head;
252 
253  for (uint32_t i=0; i<ppbmrm.nCP; ++i)
254  {
255  A_1=get_cutting_plane(cp_ptr);
256  cp_ptr=cp_ptr->next;
257  rsum=CMath::dot(A_1, A_2, nDim);
258 
259  H[LIBBMRM_INDEX(ppbmrm.nCP, i, BufSize)]
260  = H[LIBBMRM_INDEX(i, ppbmrm.nCP, BufSize)]
261  = rsum;
262  }
263  }
264 
265  A_2=get_cutting_plane(CPList_tail);
266  rsum = CMath::dot(A_2, A_2, nDim);
267 
268  H[LIBBMRM_INDEX(ppbmrm.nCP, ppbmrm.nCP, BufSize)]=rsum;
269 
270  diag_H[ppbmrm.nCP]=H[LIBBMRM_INDEX(ppbmrm.nCP, ppbmrm.nCP, BufSize)];
271  I[ppbmrm.nCP]=1;
272 
273  beta[ppbmrm.nCP]=0.0; // [beta; 0]
274  ppbmrm.nCP++;
275 
276  /* tune alpha cycle */
277  /* ---------------------------------------------------------------------- */
278 
279  flag=true;
280  isThereGoodSolution=false;
281  LIBBMRM_MEMCPY(beta_start, beta, ppbmrm.nCP*sizeof(float64_t));
282  LIBBMRM_MEMCPY(I_start, I, ppbmrm.nCP*sizeof(uint32_t));
283  qp_cnt=0;
284  alpha_good=alpha;
285 
286  if (tuneAlpha)
287  {
288  alpha_start=alpha; alpha=0.0;
289  beta[ppbmrm.nCP]=0.0;
290  LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*sizeof(uint32_t));
291  I2[ppbmrm.nCP]=1;
292 
293  /* add alpha-dependent terms to H, diag_h and b */
294  cp_ptr=CPList_head;
295 
296  for (uint32_t i=0; i<ppbmrm.nCP; ++i)
297  {
298  A_1=get_cutting_plane(cp_ptr);
299  cp_ptr=cp_ptr->next;
300 
301  rsum = CMath::dot(A_1, prevW, nDim);
302 
303  b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
304  diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
305 
306  for (uint32_t j=0; j<ppbmrm.nCP; ++j)
307  H2[LIBBMRM_INDEX(i, j, BufSize)]=
308  H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha);
309  }
310 
311  /* solve QP with current alpha */
312  qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta,
313  ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
314  ppbmrm.qp_exitflag=qp_exitflag.exitflag;
315  qp_cnt++;
316  Fd_alpha0=-qp_exitflag.QP;
317 
318  /* obtain w_t and check if norm(w_{t+1} -w_t) <= K */
319  for (uint32_t i=0; i<nDim; ++i)
320  {
321  rsum=0.0;
322  cp_ptr=CPList_head;
323 
324  for (uint32_t j=0; j<ppbmrm.nCP; ++j)
325  {
326  A_1=get_cutting_plane(cp_ptr);
327  cp_ptr=cp_ptr->next;
328  rsum+=A_1[i]*beta[j];
329  }
330 
331  wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha);
332  }
333 
334  sq_norm_Wdiff=0.0;
335 
336  for (uint32_t i=0; i<nDim; ++i)
337  sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]);
338 
339  if (CMath::sqrt(sq_norm_Wdiff) <= K)
340  {
341  flag=false;
342 
343  if (alpha!=alpha_start)
344  alphaChanged=true;
345  }
346  else
347  {
348  alpha=alpha_start;
349  }
350 
351  while(flag)
352  {
353  LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*sizeof(uint32_t));
354  LIBBMRM_MEMCPY(beta, beta_start, ppbmrm.nCP*sizeof(float64_t));
355  I2[ppbmrm.nCP]=1;
356  beta[ppbmrm.nCP]=0.0;
357 
358  /* add alpha-dependent terms to H, diag_h and b */
359  cp_ptr=CPList_head;
360 
361  for (uint32_t i=0; i<ppbmrm.nCP; ++i)
362  {
363  A_1=get_cutting_plane(cp_ptr);
364  cp_ptr=cp_ptr->next;
365 
366  rsum = CMath::dot(A_1, prevW, nDim);
367 
368  b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
369  diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
370 
371  for (uint32_t j=0; j<ppbmrm.nCP; ++j)
372  H2[LIBBMRM_INDEX(i, j, BufSize)]=H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha);
373  }
374 
375  /* solve QP with current alpha */
376  qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta,
377  ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
378  ppbmrm.qp_exitflag=qp_exitflag.exitflag;
379  qp_cnt++;
380 
381  /* obtain w_t and check if norm(w_{t+1}-w_t) <= K */
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<ppbmrm.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  wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha);
395  }
396 
397  sq_norm_Wdiff=0.0;
398  for (uint32_t i=0; i<nDim; ++i)
399  sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]);
400 
401  if (CMath::sqrt(sq_norm_Wdiff) > K)
402  {
403  /* if there is a record of some good solution
404  * (i.e. adjust alpha by division by 2) */
405 
406  if (isThereGoodSolution)
407  {
408  LIBBMRM_MEMCPY(beta, beta_good, ppbmrm.nCP*sizeof(float64_t));
409  LIBBMRM_MEMCPY(I2, I_good, ppbmrm.nCP*sizeof(uint32_t));
410  alpha=alpha_good;
411  qp_exitflag=qp_exitflag_good;
412  flag=false;
413  }
414  else
415  {
416  if (alpha == 0)
417  {
418  alpha=1.0;
419  alphaChanged=true;
420  }
421  else
422  {
423  alpha*=2;
424  alphaChanged=true;
425  }
426  }
427  }
428  else
429  {
430  if (alpha > 0)
431  {
432  /* keep good solution and try for alpha /= 2 if previous alpha was 1 */
433  LIBBMRM_MEMCPY(beta_good, beta, ppbmrm.nCP*sizeof(float64_t));
434  LIBBMRM_MEMCPY(I_good, I2, ppbmrm.nCP*sizeof(uint32_t));
435  alpha_good=alpha;
436  qp_exitflag_good=qp_exitflag;
437  isThereGoodSolution=true;
438 
439  if (alpha!=1.0)
440  {
441  alpha/=2.0;
442  alphaChanged=true;
443  }
444  else
445  {
446  alpha=0.0;
447  alphaChanged=true;
448  }
449  }
450  else
451  {
452  flag=false;
453  }
454  }
455  }
456  }
457  else
458  {
459  alphaChanged=false;
460  LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*sizeof(uint32_t));
461  LIBBMRM_MEMCPY(beta, beta_start, ppbmrm.nCP*sizeof(float64_t));
462 
463  /* add alpha-dependent terms to H, diag_h and b */
464  cp_ptr=CPList_head;
465 
466  for (uint32_t i=0; i<ppbmrm.nCP; ++i)
467  {
468  A_1=get_cutting_plane(cp_ptr);
469  cp_ptr=cp_ptr->next;
470 
471  rsum = CMath::dot(A_1, prevW, nDim);
472 
473  b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
474  diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
475 
476  for (uint32_t j=0; j<ppbmrm.nCP; ++j)
477  H2[LIBBMRM_INDEX(i, j, BufSize)]=
478  H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha);
479  }
480  /* solve QP with current alpha */
481  qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta,
482  ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
483  ppbmrm.qp_exitflag=qp_exitflag.exitflag;
484  qp_cnt++;
485  }
486 
487  /* ----------------------------------------------------------------------------------------------- */
488 
489  /* Update ICPcounter (add one to unused and reset used) + compute number of active CPs */
490  ppbmrm.nzA=0;
491 
492  for (uint32_t aaa=0; aaa<ppbmrm.nCP; ++aaa)
493  {
494  if (beta[aaa]>epsilon)
495  {
496  ++ppbmrm.nzA;
497  icp_stats.ICPcounter[aaa]=0;
498  }
499  else
500  {
501  icp_stats.ICPcounter[aaa]+=1;
502  }
503  }
504 
505  /* W update */
506  for (uint32_t i=0; i<nDim; ++i)
507  {
508  rsum=0.0;
509  cp_ptr=CPList_head;
510 
511  for (uint32_t j=0; j<ppbmrm.nCP; ++j)
512  {
513  A_1=get_cutting_plane(cp_ptr);
514  cp_ptr=cp_ptr->next;
515  rsum+=A_1[i]*beta[j];
516  }
517 
518  W[i]=(2*alpha*prevW[i]-rsum)/(_lambda+2*alpha);
519  }
520 
521  /* risk and subgradient computation */
522  R = machine->risk(subgrad, W);
523  add_cutting_plane(&CPList_tail, map, A,
524  find_free_idx(map, BufSize), subgrad, nDim);
525 
526  sq_norm_W=CMath::dot(W, W, nDim);
527  sq_norm_prevW=CMath::dot(prevW, prevW, nDim);
528  b[ppbmrm.nCP]=CMath::dot(subgrad, W, nDim) - R;
529 
530  sq_norm_Wdiff=0.0;
531  for (uint32_t j=0; j<nDim; ++j)
532  {
533  sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
534  }
535 
536  /* compute Fp and Fd */
537  ppbmrm.Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff;
538  ppbmrm.Fd=-qp_exitflag.QP+((alpha*_lambda)/(_lambda + 2*alpha))*sq_norm_prevW;
539 
540  /* gamma + tuneAlpha flag */
541  if (alphaChanged)
542  {
543  eps=1.0-(ppbmrm.Fd/ppbmrm.Fp);
544  gamma=(lastFp*(1-eps)-Fd_alpha0)/(Tmax*(1-eps));
545  }
546 
547  if ((lastFp-ppbmrm.Fp) <= gamma)
548  {
549  tuneAlpha=true;
550  }
551  else
552  {
553  tuneAlpha=false;
554  }
555 
556  /* Stopping conditions - set only with nonzero alpha */
557  if (alpha==0.0)
558  {
559  if (ppbmrm.Fp-ppbmrm.Fd<=TolRel*LIBBMRM_ABS(ppbmrm.Fp))
560  ppbmrm.exitflag=1;
561 
562  if (ppbmrm.Fp-ppbmrm.Fd<=TolAbs)
563  ppbmrm.exitflag=2;
564  }
565 
566  if (ppbmrm.nCP>=BufSize)
567  ppbmrm.exitflag=-1;
568 
569  tstop=ttime.cur_time_diff(false);
570 
571  /* compute wdist (= || W_{t+1} - W_{t} || ) */
572  sq_norm_Wdiff=0.0;
573 
574  for (uint32_t i=0; i<nDim; ++i)
575  {
576  sq_norm_Wdiff+=(W[i]-prevW[i])*(W[i]-prevW[i]);
577  }
578 
579  wdist=CMath::sqrt(sq_norm_Wdiff);
580 
581  /* Keep history of Fp, Fd, wdist */
582  ppbmrm.hist_Fp[ppbmrm.nIter]=ppbmrm.Fp;
583  ppbmrm.hist_Fd[ppbmrm.nIter]=ppbmrm.Fd;
584  ppbmrm.hist_wdist[ppbmrm.nIter]=wdist;
585 
586  /* Verbose output */
587  if (verbose)
588  SG_SDEBUG("%4d: tim=%.3lf, Fp=%lf, Fd=%lf, (Fp-Fd)=%lf, (Fp-Fd)/Fp=%lf, R=%lf, nCP=%d, nzA=%d, wdist=%lf, alpha=%lf, qp_cnt=%d, gamma=%lf, tuneAlpha=%d\n",
589  ppbmrm.nIter, tstop-tstart, ppbmrm.Fp, ppbmrm.Fd, ppbmrm.Fp-ppbmrm.Fd,
590  (ppbmrm.Fp-ppbmrm.Fd)/ppbmrm.Fp, R, ppbmrm.nCP, ppbmrm.nzA, wdist, alpha,
591  qp_cnt, gamma, tuneAlpha);
592 
593  /* Check size of Buffer */
594  if (ppbmrm.nCP>=BufSize)
595  {
596  ppbmrm.exitflag=-2;
597  SG_SERROR("Buffer exceeded.\n")
598  }
599 
600  /* keep w_t + Fp */
601  LIBBMRM_MEMCPY(prevW, W, nDim*sizeof(float64_t));
602  lastFp=ppbmrm.Fp;
603 
604  /* Inactive Cutting Planes (ICP) removal */
605  if (cleanICP)
606  {
607  clean_icp(&icp_stats, ppbmrm, &CPList_head, &CPList_tail, H, diag_H, beta, map, cleanAfter, b, I);
608  }
609 
610  // next CP would exceed BufSize
611  if (ppbmrm.nCP+1 >= BufSize)
612  ppbmrm.exitflag=-1;
613 
614  /* Debug: compute objective and training error */
615  if (verbose)
616  {
617  SGVector<float64_t> w_debug(W, nDim, false);
618  float64_t primal = CSOSVMHelper::primal_objective(w_debug, model, _lambda);
619  float64_t train_error = CSOSVMHelper::average_loss(w_debug, model);
620  helper->add_debug_info(primal, ppbmrm.nIter, train_error);
621  }
622  } /* end of main loop */
623 
624  if (verbose)
625  {
626  helper->terminate();
627  SG_UNREF(helper);
628  }
629 
630  ppbmrm.hist_Fp.resize_vector(ppbmrm.nIter);
631  ppbmrm.hist_Fd.resize_vector(ppbmrm.nIter);
632  ppbmrm.hist_wdist.resize_vector(ppbmrm.nIter);
633 
634  cp_ptr=CPList_head;
635 
636  while(cp_ptr!=NULL)
637  {
638  cp_ptr2=cp_ptr;
639  cp_ptr=cp_ptr->next;
640  LIBBMRM_FREE(cp_ptr2);
641  cp_ptr2=NULL;
642  }
643 
644  cp_list=NULL;
645 
646 cleanup:
647 
648  LIBBMRM_FREE(H);
649  LIBBMRM_FREE(b);
650  LIBBMRM_FREE(beta);
651  LIBBMRM_FREE(A);
652  LIBBMRM_FREE(subgrad);
653  LIBBMRM_FREE(diag_H);
654  LIBBMRM_FREE(I);
655  LIBBMRM_FREE(icp_stats.ICPcounter);
656  LIBBMRM_FREE(icp_stats.ICPs);
657  LIBBMRM_FREE(icp_stats.ACPs);
658  LIBBMRM_FREE(icp_stats.H_buff);
659  LIBBMRM_FREE(map);
660  LIBBMRM_FREE(prevW);
661  LIBBMRM_FREE(wt);
662  LIBBMRM_FREE(beta_start);
663  LIBBMRM_FREE(beta_good);
664  LIBBMRM_FREE(I_start);
665  LIBBMRM_FREE(I_good);
666  LIBBMRM_FREE(I2);
667  LIBBMRM_FREE(b2);
668  LIBBMRM_FREE(diag_H2);
669  LIBBMRM_FREE(H2);
670 
671  if (cp_list)
672  LIBBMRM_FREE(cp_list);
673 
674  SG_UNREF(model);
675 
676  return(ppbmrm);
677 }
678 }
679 #endif //USE_GPL_SHOGUN
REQUIRE
#define REQUIRE(x,...)
Definition: SGIO.h:206
shogun::CSOSVMHelper::primal_objective
static float64_t primal_objective(SGVector< float64_t > w, CStructuredModel *model, float64_t lbda)
Definition: SOSVMHelper.cpp:57
ASSERT
#define ASSERT(x)
Definition: SGIO.h:201
float64_t
double float64_t
Definition: common.h:50
floatmax_t
long double floatmax_t
Definition: common.h:51
shogun::CMath::dot
static float64_t dot(const bool *v1, const bool *v2, int32_t n)
Compute dot product between v1 and v2 (blas optimized)
Definition: Math.h:627
shogun::CSOSVMHelper::average_loss
static float64_t average_loss(SGVector< float64_t > w, CStructuredModel *model, bool is_ub=false)
Definition: SOSVMHelper.cpp:88
SG_UNREF
#define SG_UNREF(x)
Definition: SGObject.h:55
shogun
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
SG_SDEBUG
#define SG_SDEBUG(...)
Definition: SGIO.h:168
SG_SERROR
#define SG_SERROR(...)
Definition: SGIO.h:179
shogun::linalg::max
Matrix::Scalar max(Matrix m)
Definition: Redux.h:68
shogun::CMath::sqrt
static float32_t sqrt(float32_t x)
Definition: Math.h:459
SHOGUN Machine Learning Toolbox - Documentation