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

SHOGUN Machine Learning Toolbox - Documentation