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

SHOGUN 机器学习工具包 - 项目文档