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

SHOGUN Machine Learning Toolbox - Documentation