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