Tron.cpp

Go to the documentation of this file.
00001 #include <math.h>
00002 #include <stdio.h>
00003 #include <string.h>
00004 #include <stdarg.h>
00005 
00006 #include "lib/config.h"
00007 
00008 #ifdef HAVE_LAPACK
00009 #include "lib/Mathematics.h"
00010 #include "classifier/svm/Tron.h"
00011 
00012 CTron::CTron(const function *f, double e, int it)
00013 : CSGObject()
00014 {
00015     this->fun_obj=const_cast<function *>(f);
00016     this->eps=e;
00017     this->max_iter=it;
00018 }
00019 
00020 CTron::~CTron()
00021 {
00022 }
00023 
00024 void CTron::tron(double *w)
00025 {
00026     // Parameters for updating the iterates.
00027     double eta0 = 1e-4, eta1 = 0.25, eta2 = 0.75;
00028 
00029     // Parameters for updating the trust region size delta.
00030     double sigma1 = 0.25, sigma2 = 0.5, sigma3 = 4;
00031 
00032     int n = fun_obj->get_nr_variable();
00033     int i, cg_iter;
00034     double delta, snorm, one=1.0;
00035     double alpha, f, fnew, prered, actred, gs;
00036     int search = 1, iter = 1, inc = 1;
00037     double *s = new double[n];
00038     double *r = new double[n];
00039     double *w_new = new double[n];
00040     double *g = new double[n];
00041 
00042     for (i=0; i<n; i++)
00043         w[i] = 0;
00044 
00045         f = fun_obj->fun(w);
00046     fun_obj->grad(w, g);
00047     delta = cblas_dnrm2(n, g, inc);
00048     double gnorm1 = delta;
00049     double gnorm = gnorm1;
00050 
00051     if (gnorm <= eps*gnorm1)
00052         search = 0;
00053 
00054     iter = 1;
00055 
00056     while (iter <= max_iter && search)
00057     {
00058         cg_iter = trcg(delta, g, s, r);
00059 
00060         memcpy(w_new, w, sizeof(double)*n);
00061         cblas_daxpy(n, one, s, inc, w_new, inc);
00062 
00063         gs = cblas_ddot(n, g, inc, s, inc);
00064         prered = -0.5*(gs-cblas_ddot(n, s, inc, r, inc));
00065                 fnew = fun_obj->fun(w_new);
00066 
00067         // Compute the actual reduction.
00068             actred = f - fnew;
00069 
00070         // On the first iteration, adjust the initial step bound.
00071         snorm = cblas_dnrm2(n, s, inc);
00072         if (iter == 1)
00073             delta = CMath::min(delta, snorm);
00074 
00075         // Compute prediction alpha*snorm of the step.
00076         if (fnew - f - gs <= 0)
00077             alpha = sigma3;
00078         else
00079             alpha = CMath::max(sigma1, -0.5*(gs/(fnew - f - gs)));
00080 
00081         // Update the trust region bound according to the ratio of actual to predicted reduction.
00082         if (actred < eta0*prered)
00083             delta = CMath::min(CMath::max(alpha, sigma1)*snorm, sigma2*delta);
00084         else if (actred < eta1*prered)
00085             delta = CMath::max(sigma1*delta, CMath::min(alpha*snorm, sigma2*delta));
00086         else if (actred < eta2*prered)
00087             delta = CMath::max(sigma1*delta, CMath::min(alpha*snorm, sigma3*delta));
00088         else
00089             delta = CMath::max(delta, CMath::min(alpha*snorm, sigma3*delta));
00090 
00091         SG_INFO("iter %2d act %5.3e pre %5.3e delta %5.3e f %5.3e |g| %5.3e CG %3d\n", iter, actred, prered, delta, f, gnorm, cg_iter);
00092 
00093         if (actred > eta0*prered)
00094         {
00095             iter++;
00096             memcpy(w, w_new, sizeof(double)*n);
00097             f = fnew;
00098                 fun_obj->grad(w, g);
00099 
00100             gnorm = cblas_dnrm2(n, g, inc);
00101             if (gnorm < eps*gnorm1)
00102                 break;
00103         }
00104         if (f < -1.0e+32)
00105         {
00106             SG_WARNING("f < -1.0e+32\n");
00107             break;
00108         }
00109         if (CMath::abs(actred) <= 0 && CMath::abs(prered) <= 0)
00110         {
00111             SG_WARNING("actred and prered <= 0\n");
00112             break;
00113         }
00114         if (CMath::abs(actred) <= 1.0e-12*CMath::abs(f) &&
00115             CMath::abs(prered) <= 1.0e-12*CMath::abs(f))
00116         {
00117             SG_WARNING("actred and prered too small\n");
00118             break;
00119         }
00120     }
00121 
00122     delete[] g;
00123     delete[] r;
00124     delete[] w_new;
00125     delete[] s;
00126 }
00127 
00128 int CTron::trcg(double delta, double *g, double *s, double *r)
00129 {
00130     int i, inc = 1;
00131     int n = fun_obj->get_nr_variable();
00132     double one = 1;
00133     double *d = new double[n];
00134     double *Hd = new double[n];
00135     double rTr, rnewTrnew, alpha, beta, cgtol;
00136 
00137     for (i=0; i<n; i++)
00138     {
00139         s[i] = 0;
00140         r[i] = -g[i];
00141         d[i] = r[i];
00142     }
00143     cgtol = 0.1*cblas_dnrm2(n, g, inc);
00144 
00145     int cg_iter = 0;
00146     rTr = cblas_ddot(n, r, inc, r, inc);
00147     while (1)
00148     {
00149         if (cblas_dnrm2(n, r, inc) <= cgtol)
00150             break;
00151         cg_iter++;
00152         fun_obj->Hv(d, Hd);
00153 
00154         alpha = rTr/cblas_ddot(n, d, inc, Hd, inc);
00155         cblas_daxpy(n, alpha, d, inc, s, inc);
00156         if (cblas_dnrm2(n, s, inc) > delta)
00157         {
00158             SG_INFO("cg reaches trust region boundary\n");
00159             alpha = -alpha;
00160             cblas_daxpy(n, alpha, d, inc, s, inc);
00161 
00162             double std = cblas_ddot(n, s, inc, d, inc);
00163             double sts = cblas_ddot(n, s, inc, s, inc);
00164             double dtd = cblas_ddot(n, d, inc, d, inc);
00165             double dsq = delta*delta;
00166             double rad = sqrt(std*std + dtd*(dsq-sts));
00167             if (std >= 0)
00168                 alpha = (dsq - sts)/(std + rad);
00169             else
00170                 alpha = (rad - std)/dtd;
00171             cblas_daxpy(n, alpha, d, inc, s, inc);
00172             alpha = -alpha;
00173             cblas_daxpy(n, alpha, Hd, inc, r, inc);
00174             break;
00175         }
00176         alpha = -alpha;
00177         cblas_daxpy(n, alpha, Hd, inc, r, inc);
00178         rnewTrnew = cblas_ddot(n, r, inc, r, inc);
00179         beta = rnewTrnew/rTr;
00180         cblas_dscal(n, beta, d, inc);
00181         cblas_daxpy(n, one, r, inc, d, inc);
00182         rTr = rnewTrnew;
00183     }
00184 
00185     delete[] d;
00186     delete[] Hd;
00187 
00188     return(cg_iter);
00189 }
00190 
00191 double CTron::norm_inf(int n, double *x)
00192 {
00193     double dmax = CMath::abs(x[0]);
00194     for (int i=1; i<n; i++)
00195         if (CMath::abs(x[i]) >= dmax)
00196             dmax = CMath::abs(x[i]);
00197     return(dmax);
00198 }
00199 #endif //HAVE_LAPACK

SHOGUN Machine Learning Toolbox - Documentation