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

SHOGUN Machine Learning Toolbox - Documentation