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

SHOGUN Machine Learning Toolbox - Documentation