SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
libppbm.cpp
Go to the documentation of this file.
00001 /*
00002  * This program is free software; you can redistribute it and/or modify
00003  * it under the terms of the GNU General Public License as published by
00004  * the Free Software Foundation; either version 3 of the License, or
00005  * (at your option) any later version.
00006  *
00007  * libp3bm.h: Implementation of the Proximal Point P-BMRM solver for SO training
00008  *
00009  * Copyright (C) 2012 Michal Uricar, uricamic@cmp.felk.cvut.cz
00010  *
00011  * Implementation of the Proximal Point P-BMRM
00012  *--------------------------------------------------------------------- */
00013 
00014 #include <shogun/structure/libp3bm.h>
00015 #include <shogun/lib/external/libqp.h>
00016 #include <shogun/lib/Time.h>
00017 
00018 namespace shogun
00019 {
00020 static const uint32_t QPSolverMaxIter=0xFFFFFFFF;
00021 static const float64_t epsilon=0.0;
00022 
00023 static float64_t *H, *H2;
00024 static uint32_t BufSize;
00025 
00026 /*----------------------------------------------------------------------
00027   Returns pointer at i-th column of Hessian matrix.
00028   ----------------------------------------------------------------------*/
00029 static const float64_t *get_col( uint32_t i)
00030 {
00031     return( &H2[ BufSize*i ] );
00032 }
00033 
00034 BmrmStatistics svm_ppbm_solver(
00035         CDualLibQPBMSOSVM *machine,
00036         float64_t*      W,
00037         float64_t       TolRel,
00038         float64_t       TolAbs,
00039         float64_t       _lambda,
00040         uint32_t        _BufSize,
00041         bool            cleanICP,
00042         uint32_t        cleanAfter,
00043         float64_t       K,
00044         uint32_t        Tmax,
00045         bool            verbose)
00046 {
00047     BmrmStatistics ppbmrm;
00048     libqp_state_T qp_exitflag={0, 0, 0, 0}, qp_exitflag_good={0, 0, 0, 0};
00049     float64_t *b, *b2, *beta, *beta_good, *beta_start, *diag_H, *diag_H2;
00050     float64_t R, *subgrad, *A, QPSolverTolRel, C=1.0;
00051     float64_t *prevW, *wt, alpha, alpha_start, alpha_good=0.0, Fd_alpha0=0.0;
00052     float64_t lastFp, wdist, gamma=0.0;
00053     floatmax_t rsum, sq_norm_W, sq_norm_Wdiff, sq_norm_prevW, eps;
00054     uint32_t *I, *I2, *I_start, *I_good;
00055     uint8_t S=1;
00056     CStructuredModel* model=machine->get_model();
00057     uint32_t nDim=model->get_dim();
00058     CSOSVMHelper* helper = NULL;
00059     uint32_t qp_cnt=0;
00060     bmrm_ll *CPList_head, *CPList_tail, *cp_ptr, *cp_ptr2, *cp_list=NULL;
00061     float64_t *A_1=NULL, *A_2=NULL;
00062     bool *map=NULL, tuneAlpha=true, flag=true, alphaChanged=false, isThereGoodSolution=false;
00063 
00064     CTime ttime;
00065     float64_t tstart, tstop;
00066 
00067 
00068     tstart=ttime.cur_time_diff(false);
00069 
00070     BufSize=_BufSize;
00071     QPSolverTolRel=1e-9;
00072 
00073     H=NULL;
00074     b=NULL;
00075     beta=NULL;
00076     A=NULL;
00077     subgrad=NULL;
00078     diag_H=NULL;
00079     I=NULL;
00080     prevW=NULL;
00081     wt=NULL;
00082     diag_H2=NULL;
00083     b2=NULL;
00084     I2=NULL;
00085     H2=NULL;
00086     I_good=NULL;
00087     I_start=NULL;
00088     beta_start=NULL;
00089     beta_good=NULL;
00090 
00091     alpha=0.0;
00092 
00093     H= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, float64_t);
00094 
00095     A= (float64_t*) LIBBMRM_CALLOC(nDim*BufSize, float64_t);
00096 
00097     b= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
00098 
00099     beta= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
00100 
00101     subgrad= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t);
00102 
00103     diag_H= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
00104 
00105     I= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
00106 
00107     /* structure for maintaining inactive CPs info */
00108     ICP_stats icp_stats;
00109     icp_stats.maxCPs = BufSize;
00110     icp_stats.ICPcounter= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
00111     icp_stats.ICPs= (float64_t**) LIBBMRM_CALLOC(BufSize, float64_t*);
00112     icp_stats.ACPs= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
00113 
00114     cp_list= (bmrm_ll*) LIBBMRM_CALLOC(1, bmrm_ll);
00115 
00116     prevW= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t);
00117 
00118     wt= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t);
00119 
00120     if (H==NULL || A==NULL || b==NULL || beta==NULL || subgrad==NULL ||
00121             diag_H==NULL || I==NULL || icp_stats.ICPcounter==NULL ||
00122             icp_stats.ICPs==NULL || icp_stats.ACPs==NULL ||
00123             cp_list==NULL || prevW==NULL || wt==NULL)
00124     {
00125         ppbmrm.exitflag=-2;
00126         goto cleanup;
00127     }
00128 
00129     map= (bool*) LIBBMRM_CALLOC(BufSize, bool);
00130 
00131     if (map==NULL)
00132     {
00133         ppbmrm.exitflag=-2;
00134         goto cleanup;
00135     }
00136 
00137     memset( (bool*) map, true, BufSize);
00138 
00139     /* Temporary buffers for ICP removal */
00140     icp_stats.H_buff= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, float64_t);
00141 
00142     if (icp_stats.H_buff==NULL)
00143     {
00144         ppbmrm.exitflag=-2;
00145         goto cleanup;
00146     }
00147 
00148     /* Temporary buffers */
00149     beta_start= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
00150 
00151     beta_good= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
00152 
00153     b2= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
00154 
00155     diag_H2= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
00156 
00157     H2= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, float64_t);
00158 
00159     I_start= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
00160 
00161     I_good= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
00162 
00163     I2= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
00164 
00165     if (beta_start==NULL || beta_good==NULL || b2==NULL || diag_H2==NULL ||
00166             I_start==NULL || I_good==NULL || I2==NULL || H2==NULL)
00167     {
00168         ppbmrm.exitflag=-2;
00169         goto cleanup;
00170     }
00171 
00172     ppbmrm.hist_Fp.resize_vector(BufSize);
00173     ppbmrm.hist_Fd.resize_vector(BufSize);
00174     ppbmrm.hist_wdist.resize_vector(BufSize);
00175 
00176     /* Iinitial solution */
00177     R = machine->risk(subgrad, W);
00178 
00179     ppbmrm.nCP=0;
00180     ppbmrm.nIter=0;
00181     ppbmrm.exitflag=0;
00182 
00183     b[0]=-R;
00184 
00185     /* Cutting plane auxiliary double linked list */
00186     LIBBMRM_MEMCPY(A, subgrad, nDim*sizeof(float64_t));
00187     map[0]=false;
00188     cp_list->address=&A[0];
00189     cp_list->idx=0;
00190     cp_list->prev=NULL;
00191     cp_list->next=NULL;
00192     CPList_head=cp_list;
00193     CPList_tail=cp_list;
00194 
00195     /* Compute initial value of Fp, Fd, assuming that W is zero vector */
00196     sq_norm_Wdiff=0.0;
00197 
00198     b[0] = SGVector<float64_t>::dot(subgrad, W, nDim);
00199     sq_norm_W = SGVector<float64_t>::dot(W, W, nDim);
00200     for (uint32_t j=0; j<nDim; ++j)
00201     {
00202         sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
00203     }
00204 
00205     ppbmrm.Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff;
00206     ppbmrm.Fd=-LIBBMRM_PLUS_INF;
00207     lastFp=ppbmrm.Fp;
00208     wdist=CMath::sqrt(sq_norm_Wdiff);
00209 
00210     K = (sq_norm_W == 0.0) ? 0.4 : 0.01*CMath::sqrt(sq_norm_W);
00211 
00212     LIBBMRM_MEMCPY(prevW, W, nDim*sizeof(float64_t));
00213 
00214     tstop=ttime.cur_time_diff(false);
00215 
00216     /* Keep history of Fp, Fd, wdist */
00217     ppbmrm.hist_Fp[0]=ppbmrm.Fp;
00218     ppbmrm.hist_Fd[0]=ppbmrm.Fd;
00219     ppbmrm.hist_wdist[0]=wdist;
00220 
00221     /* Verbose output */
00222 
00223     if (verbose)
00224         SG_SDEBUG("%4d: tim=%.3lf, Fp=%lf, Fd=%lf, R=%lf, K=%lf\n",
00225                 ppbmrm.nIter, tstop-tstart, ppbmrm.Fp, ppbmrm.Fd, R, K);
00226 
00227     if (verbose)
00228         helper = machine->get_helper();
00229 
00230     /* main loop */
00231 
00232     while (ppbmrm.exitflag==0)
00233     {
00234         tstart=ttime.cur_time_diff(false);
00235         ppbmrm.nIter++;
00236 
00237         /* Update H */
00238 
00239         if (ppbmrm.nCP>0)
00240         {
00241             A_2=get_cutting_plane(CPList_tail);
00242             cp_ptr=CPList_head;
00243 
00244             for (uint32_t i=0; i<ppbmrm.nCP; ++i)
00245             {
00246                 A_1=get_cutting_plane(cp_ptr);
00247                 cp_ptr=cp_ptr->next;
00248                 rsum=SGVector<float64_t>::dot(A_1, A_2, nDim);
00249 
00250                 H[LIBBMRM_INDEX(ppbmrm.nCP, i, BufSize)]
00251                     = H[LIBBMRM_INDEX(i, ppbmrm.nCP, BufSize)]
00252                     = rsum;
00253             }
00254         }
00255 
00256         A_2=get_cutting_plane(CPList_tail);
00257         rsum = SGVector<float64_t>::dot(A_2, A_2, nDim);
00258 
00259         H[LIBBMRM_INDEX(ppbmrm.nCP, ppbmrm.nCP, BufSize)]=rsum;
00260 
00261         diag_H[ppbmrm.nCP]=H[LIBBMRM_INDEX(ppbmrm.nCP, ppbmrm.nCP, BufSize)];
00262         I[ppbmrm.nCP]=1;
00263 
00264         ppbmrm.nCP++;
00265         beta[ppbmrm.nCP]=0.0; // [beta; 0]
00266 
00267         /* tune alpha cycle */
00268         /* ---------------------------------------------------------------------- */
00269 
00270         flag=true;
00271         isThereGoodSolution=false;
00272         LIBBMRM_MEMCPY(beta_start, beta, ppbmrm.nCP*sizeof(float64_t));
00273         LIBBMRM_MEMCPY(I_start, I, ppbmrm.nCP*sizeof(uint32_t));
00274         qp_cnt=0;
00275         alpha_good=alpha;
00276 
00277         if (tuneAlpha)
00278         {
00279             alpha_start=alpha; alpha=0.0;
00280             beta[ppbmrm.nCP]=0.0;
00281             LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*sizeof(uint32_t));
00282             I2[ppbmrm.nCP]=1;
00283 
00284             /* add alpha-dependent terms to H, diag_h and b */
00285             cp_ptr=CPList_head;
00286 
00287             for (uint32_t i=0; i<ppbmrm.nCP; ++i)
00288             {
00289                 A_1=get_cutting_plane(cp_ptr);
00290                 cp_ptr=cp_ptr->next;
00291 
00292                 rsum = SGVector<float64_t>::dot(A_1, prevW, nDim);
00293 
00294                 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
00295                 diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
00296 
00297                 for (uint32_t j=0; j<ppbmrm.nCP; ++j)
00298                     H2[LIBBMRM_INDEX(i, j, BufSize)]=
00299                         H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha);
00300             }
00301 
00302             /* solve QP with current alpha */
00303             qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta,
00304                     ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
00305             ppbmrm.qp_exitflag=qp_exitflag.exitflag;
00306             qp_cnt++;
00307             Fd_alpha0=-qp_exitflag.QP;
00308 
00309             /* obtain w_t and check if norm(w_{t+1} -w_t) <= K */
00310             for (uint32_t i=0; i<nDim; ++i)
00311             {
00312                 rsum=0.0;
00313                 cp_ptr=CPList_head;
00314 
00315                 for (uint32_t j=0; j<ppbmrm.nCP; ++j)
00316                 {
00317                     A_1=get_cutting_plane(cp_ptr);
00318                     cp_ptr=cp_ptr->next;
00319                     rsum+=A_1[i]*beta[j];
00320                 }
00321 
00322                 wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha);
00323             }
00324 
00325             sq_norm_Wdiff=0.0;
00326 
00327             for (uint32_t i=0; i<nDim; ++i)
00328                 sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]);
00329 
00330             if (CMath::sqrt(sq_norm_Wdiff) <= K)
00331             {
00332                 flag=false;
00333 
00334                 if (alpha!=alpha_start)
00335                     alphaChanged=true;
00336             }
00337             else
00338             {
00339                 alpha=alpha_start;
00340             }
00341 
00342             while(flag)
00343             {
00344                 LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*sizeof(uint32_t));
00345                 LIBBMRM_MEMCPY(beta, beta_start, ppbmrm.nCP*sizeof(float64_t));
00346                 I2[ppbmrm.nCP]=1;
00347                 beta[ppbmrm.nCP]=0.0;
00348 
00349                 /* add alpha-dependent terms to H, diag_h and b */
00350                 cp_ptr=CPList_head;
00351 
00352                 for (uint32_t i=0; i<ppbmrm.nCP; ++i)
00353                 {
00354                     A_1=get_cutting_plane(cp_ptr);
00355                     cp_ptr=cp_ptr->next;
00356 
00357                     rsum = SGVector<float64_t>::dot(A_1, prevW, nDim);
00358 
00359                     b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
00360                     diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
00361 
00362                     for (uint32_t j=0; j<ppbmrm.nCP; ++j)
00363                         H2[LIBBMRM_INDEX(i, j, BufSize)]=H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha);
00364                 }
00365 
00366                 /* solve QP with current alpha */
00367                 qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta,
00368                         ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
00369                 ppbmrm.qp_exitflag=qp_exitflag.exitflag;
00370                 qp_cnt++;
00371 
00372                 /* obtain w_t and check if norm(w_{t+1}-w_t) <= K */
00373                 for (uint32_t i=0; i<nDim; ++i)
00374                 {
00375                     rsum=0.0;
00376                     cp_ptr=CPList_head;
00377 
00378                     for (uint32_t j=0; j<ppbmrm.nCP; ++j)
00379                     {
00380                         A_1=get_cutting_plane(cp_ptr);
00381                         cp_ptr=cp_ptr->next;
00382                         rsum+=A_1[i]*beta[j];
00383                     }
00384 
00385                     wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha);
00386                 }
00387 
00388                 sq_norm_Wdiff=0.0;
00389                 for (uint32_t i=0; i<nDim; ++i)
00390                     sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]);
00391 
00392                 if (CMath::sqrt(sq_norm_Wdiff) > K)
00393                 {
00394                     /* if there is a record of some good solution
00395                      * (i.e. adjust alpha by division by 2) */
00396 
00397                     if (isThereGoodSolution)
00398                     {
00399                         LIBBMRM_MEMCPY(beta, beta_good, ppbmrm.nCP*sizeof(float64_t));
00400                         LIBBMRM_MEMCPY(I2, I_good, ppbmrm.nCP*sizeof(uint32_t));
00401                         alpha=alpha_good;
00402                         qp_exitflag=qp_exitflag_good;
00403                         flag=false;
00404                     }
00405                     else
00406                     {
00407                         if (alpha == 0)
00408                         {
00409                             alpha=1.0;
00410                             alphaChanged=true;
00411                         }
00412                         else
00413                         {
00414                             alpha*=2;
00415                             alphaChanged=true;
00416                         }
00417                     }
00418                 }
00419                 else
00420                 {
00421                     if (alpha > 0)
00422                     {
00423                         /* keep good solution and try for alpha /= 2 if previous alpha was 1 */
00424                         LIBBMRM_MEMCPY(beta_good, beta, ppbmrm.nCP*sizeof(float64_t));
00425                         LIBBMRM_MEMCPY(I_good, I2, ppbmrm.nCP*sizeof(uint32_t));
00426                         alpha_good=alpha;
00427                         qp_exitflag_good=qp_exitflag;
00428                         isThereGoodSolution=true;
00429 
00430                         if (alpha!=1.0)
00431                         {
00432                             alpha/=2.0;
00433                             alphaChanged=true;
00434                         }
00435                         else
00436                         {
00437                             alpha=0.0;
00438                             alphaChanged=true;
00439                         }
00440                     }
00441                     else
00442                     {
00443                         flag=false;
00444                     }
00445                 }
00446             }
00447         }
00448         else
00449         {
00450             alphaChanged=false;
00451             LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*sizeof(uint32_t));
00452             LIBBMRM_MEMCPY(beta, beta_start, ppbmrm.nCP*sizeof(float64_t));
00453 
00454             /* add alpha-dependent terms to H, diag_h and b */
00455             cp_ptr=CPList_head;
00456 
00457             for (uint32_t i=0; i<ppbmrm.nCP; ++i)
00458             {
00459                 A_1=get_cutting_plane(cp_ptr);
00460                 cp_ptr=cp_ptr->next;
00461 
00462                 rsum = SGVector<float64_t>::dot(A_1, prevW, nDim);
00463 
00464                 b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
00465                 diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
00466 
00467                 for (uint32_t j=0; j<ppbmrm.nCP; ++j)
00468                     H2[LIBBMRM_INDEX(i, j, BufSize)]=
00469                         H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha);
00470             }
00471             /* solve QP with current alpha */
00472             qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta,
00473                     ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
00474             ppbmrm.qp_exitflag=qp_exitflag.exitflag;
00475             qp_cnt++;
00476         }
00477 
00478         /* ----------------------------------------------------------------------------------------------- */
00479 
00480         /* Update ICPcounter (add one to unused and reset used) + compute number of active CPs */
00481         ppbmrm.nzA=0;
00482 
00483         for (uint32_t aaa=0; aaa<ppbmrm.nCP; ++aaa)
00484         {
00485             if (beta[aaa]>epsilon)
00486             {
00487                 ++ppbmrm.nzA;
00488                 icp_stats.ICPcounter[aaa]=0;
00489             }
00490             else
00491             {
00492                 icp_stats.ICPcounter[aaa]+=1;
00493             }
00494         }
00495 
00496         /* W update */
00497         for (uint32_t i=0; i<nDim; ++i)
00498         {
00499             rsum=0.0;
00500             cp_ptr=CPList_head;
00501 
00502             for (uint32_t j=0; j<ppbmrm.nCP; ++j)
00503             {
00504                 A_1=get_cutting_plane(cp_ptr);
00505                 cp_ptr=cp_ptr->next;
00506                 rsum+=A_1[i]*beta[j];
00507             }
00508 
00509             W[i]=(2*alpha*prevW[i]-rsum)/(_lambda+2*alpha);
00510         }
00511 
00512         /* risk and subgradient computation */
00513         R = machine->risk(subgrad, W);
00514         add_cutting_plane(&CPList_tail, map, A,
00515                 find_free_idx(map, BufSize), subgrad, nDim);
00516 
00517         sq_norm_W=SGVector<float64_t>::dot(W, W, nDim);
00518         sq_norm_prevW=SGVector<float64_t>::dot(prevW, prevW, nDim);
00519         b[ppbmrm.nCP]=SGVector<float64_t>::dot(subgrad, W, nDim) - R;
00520 
00521         sq_norm_Wdiff=0.0;
00522         for (uint32_t j=0; j<nDim; ++j)
00523         {
00524             sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
00525         }
00526 
00527         /* compute Fp and Fd */
00528         ppbmrm.Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff;
00529         ppbmrm.Fd=-qp_exitflag.QP+((alpha*_lambda)/(_lambda + 2*alpha))*sq_norm_prevW;
00530 
00531         /* gamma + tuneAlpha flag */
00532         if (alphaChanged)
00533         {
00534             eps=1.0-(ppbmrm.Fd/ppbmrm.Fp);
00535             gamma=(lastFp*(1-eps)-Fd_alpha0)/(Tmax*(1-eps));
00536         }
00537 
00538         if ((lastFp-ppbmrm.Fp) <= gamma)
00539         {
00540             tuneAlpha=true;
00541         }
00542         else
00543         {
00544             tuneAlpha=false;
00545         }
00546 
00547         /* Stopping conditions - set only with nonzero alpha */
00548         if (alpha==0.0)
00549         {
00550             if (ppbmrm.Fp-ppbmrm.Fd<=TolRel*LIBBMRM_ABS(ppbmrm.Fp))
00551                 ppbmrm.exitflag=1;
00552 
00553             if (ppbmrm.Fp-ppbmrm.Fd<=TolAbs)
00554                 ppbmrm.exitflag=2;
00555         }
00556 
00557         if (ppbmrm.nCP>=BufSize)
00558             ppbmrm.exitflag=-1;
00559 
00560         tstop=ttime.cur_time_diff(false);
00561 
00562         /* compute wdist (= || W_{t+1} - W_{t} || ) */
00563         sq_norm_Wdiff=0.0;
00564 
00565         for (uint32_t i=0; i<nDim; ++i)
00566         {
00567             sq_norm_Wdiff+=(W[i]-prevW[i])*(W[i]-prevW[i]);
00568         }
00569 
00570         wdist=CMath::sqrt(sq_norm_Wdiff);
00571 
00572         /* Keep history of Fp, Fd, wdist */
00573         ppbmrm.hist_Fp[ppbmrm.nIter]=ppbmrm.Fp;
00574         ppbmrm.hist_Fd[ppbmrm.nIter]=ppbmrm.Fd;
00575         ppbmrm.hist_wdist[ppbmrm.nIter]=wdist;
00576 
00577         /* Verbose output */
00578         if (verbose)
00579             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",
00580                     ppbmrm.nIter, tstop-tstart, ppbmrm.Fp, ppbmrm.Fd, ppbmrm.Fp-ppbmrm.Fd,
00581                     (ppbmrm.Fp-ppbmrm.Fd)/ppbmrm.Fp, R, ppbmrm.nCP, ppbmrm.nzA, wdist, alpha,
00582                     qp_cnt, gamma, tuneAlpha);
00583 
00584         /* Check size of Buffer */
00585         if (ppbmrm.nCP>=BufSize)
00586         {
00587             ppbmrm.exitflag=-2;
00588             SG_SERROR("Buffer exceeded.\n")
00589         }
00590 
00591         /* keep w_t + Fp */
00592         LIBBMRM_MEMCPY(prevW, W, nDim*sizeof(float64_t));
00593         lastFp=ppbmrm.Fp;
00594 
00595         /* Inactive Cutting Planes (ICP) removal */
00596         if (cleanICP)
00597         {
00598             clean_icp(&icp_stats, ppbmrm, &CPList_head, &CPList_tail, H, diag_H, beta, map, cleanAfter, b, I);
00599         }
00600 
00601         /* Debug: compute objective and training error */
00602         if (verbose)
00603         {
00604             SGVector<float64_t> w_debug(W, nDim, false);
00605             float64_t primal = CSOSVMHelper::primal_objective(w_debug, model, _lambda);
00606             float64_t train_error = CSOSVMHelper::average_loss(w_debug, model);
00607             helper->add_debug_info(primal, ppbmrm.nIter, train_error);
00608         }
00609     } /* end of main loop */
00610 
00611     if (verbose)
00612     {
00613         helper->terminate();
00614         SG_UNREF(helper);
00615     }
00616 
00617     ppbmrm.hist_Fp.resize_vector(ppbmrm.nIter);
00618     ppbmrm.hist_Fd.resize_vector(ppbmrm.nIter);
00619     ppbmrm.hist_wdist.resize_vector(ppbmrm.nIter);
00620 
00621     cp_ptr=CPList_head;
00622 
00623     while(cp_ptr!=NULL)
00624     {
00625         cp_ptr2=cp_ptr;
00626         cp_ptr=cp_ptr->next;
00627         LIBBMRM_FREE(cp_ptr2);
00628         cp_ptr2=NULL;
00629     }
00630 
00631     cp_list=NULL;
00632 
00633 cleanup:
00634 
00635     LIBBMRM_FREE(H);
00636     LIBBMRM_FREE(b);
00637     LIBBMRM_FREE(beta);
00638     LIBBMRM_FREE(A);
00639     LIBBMRM_FREE(subgrad);
00640     LIBBMRM_FREE(diag_H);
00641     LIBBMRM_FREE(I);
00642     LIBBMRM_FREE(icp_stats.ICPcounter);
00643     LIBBMRM_FREE(icp_stats.ICPs);
00644     LIBBMRM_FREE(icp_stats.ACPs);
00645     LIBBMRM_FREE(icp_stats.H_buff);
00646     LIBBMRM_FREE(map);
00647     LIBBMRM_FREE(prevW);
00648     LIBBMRM_FREE(wt);
00649     LIBBMRM_FREE(beta_start);
00650     LIBBMRM_FREE(beta_good);
00651     LIBBMRM_FREE(I_start);
00652     LIBBMRM_FREE(I_good);
00653     LIBBMRM_FREE(I2);
00654     LIBBMRM_FREE(b2);
00655     LIBBMRM_FREE(diag_H2);
00656     LIBBMRM_FREE(H2);
00657 
00658     if (cp_list)
00659         LIBBMRM_FREE(cp_list);
00660 
00661     SG_UNREF(model);
00662 
00663     return(ppbmrm);
00664 }
00665 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation