SHOGUN
v2.0.0
|
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