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
00027 double eta0 = 1e-4, eta1 = 0.25, eta2 = 0.75;
00028
00029
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
00068 actred = f - fnew;
00069
00070
00071 snorm = cblas_dnrm2(n, s, inc);
00072 if (iter == 1)
00073 delta = CMath::min(delta, snorm);
00074
00075
00076 if (fnew - f - gs <= 0)
00077 alpha = sigma3;
00078 else
00079 alpha = CMath::max(sigma1, -0.5*(gs/(fnew - f - gs)));
00080
00081
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