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

SHOGUN Machine Learning Toolbox - Documentation