From df8a7310b15f2fa132794bd74ecb52fdfb735887 Mon Sep 17 00:00:00 2001 From: Parikshit Ram Date: Fri, 25 Jan 2008 01:29:08 +0000 Subject: [PATCH] optimizer added to svn, forgot it earlier --- fastlib/u/pram/mog_l2e/optimizers.h | 474 ++++++++++++++++++++++++++++ 1 file changed, 474 insertions(+) create mode 100644 fastlib/u/pram/mog_l2e/optimizers.h diff --git a/fastlib/u/pram/mog_l2e/optimizers.h b/fastlib/u/pram/mog_l2e/optimizers.h new file mode 100644 index 0000000000..4ec569be8d --- /dev/null +++ b/fastlib/u/pram/mog_l2e/optimizers.h @@ -0,0 +1,474 @@ +/** + * @author Parikshit Ram (pram@cc.gatech.edu) + * @file optimizers.h + * + * Implements classes for two types of optimizer + * + */ + +#ifndef OPTIMIZERS_H +#define OPTIMIZERS_H + +#include + +/** + * An optimizer using the Nelder Mead method, + * also known as the polytope or the simplex + * method. + * + * It does multivariate minimization of an + * objective function. If it is optimizing in + * 'd' dimensions, it would require 'd+1' + * starting points. + * + * Example use: + * + * @code + * double init_pts[d+1][d]; + * index_t number_of_function_evaluations; + * struct datanode *opt_module = fx_submodule(NULL,"NelderMead","opt_module"); + * Matrix data; + * index_t dim_param_space; + * + * ... + * NelderMead opt; + * opt.Init(obj_function, data, dim_param_space, opt_module); + * ... + * opt.Eval(init_pts, &number_of_function_evaluations); + * // init_pts[0] contains the optimal point found + * @endcode + * + */ +class NelderMead { + + private: + index_t dimension_; + Matrix data_; + long double (*func_ptr_)(Vector&, const Matrix&); + datanode *opt_module_; + + public: + + NelderMead() { + } + + ~NelderMead() { + } + + void Init(long double (*fun)(Vector&, const Matrix&), + Matrix& data, datanode *opt_module) { + data_.Copy(data); + func_ptr_ = fun; + opt_module_ = opt_module; + dimension_ = fx_param_int_req(opt_module_, "param_space_dim"); + } + + const Matrix& data() { + return data_; + } + + index_t dimension() { + return dimension_; + } + + void Eval(double **pts) { + + index_t dim = dimension(), num_func_eval; + index_t i, j, ihi, ilo, inhi,mpts = dim + 1; + double sum, swap, *psum; + long double swap_y, rtol, ytry, ysave, TINY = 1.0e-10; + long double *y; + Vector param_passed; + long double tol = fx_param_double(opt_module_,"tolerance", 1.0e-5); + index_t NMAX = fx_param_int(opt_module_, "MAX_FUNC_EVAL", 50000); + + param_passed.Init(dim); + psum = (double*)malloc(dim * sizeof(double)); + num_func_eval = 0; + y = (long double*)malloc(mpts*sizeof(long double)); + for(i = 0; i < mpts; i++) { + param_passed.CopyValues(pts[i]); + y[i] = (*func_ptr_)(param_passed, data()); + } + + for(;;) { + ilo = 0; + ihi = y[0] > y[1] ? (inhi = 1,0) : (inhi = 0,1); + for( i = 0; i < mpts; i++ ) { + if(y[i] <= y[ilo]) ilo = i; + if(y[i] > y[ihi]) { + inhi = ihi; + ihi = i; + } + else if((y[i] > y[inhi])&&(i != ihi)) inhi = i; + } + + rtol = 2.0 * fabs(y[ihi] - y[ilo]) / ( fabs(y[ihi]) + fabs(y[ilo]) + TINY ) ; + if(rtol < tol) { + swap_y = y[0]; + y[0] = y[ilo]; + y[ilo] = swap_y; + for( i = 0; i < dim; i++ ) { + swap = pts[0][i]; + pts[0][i] = pts[ilo][i] ; + pts[ilo][i] = swap; + } + break; + } + if(num_func_eval > NMAX){ + NOTIFY("Maximum number of function evaluations exceeded"); + break; + } + num_func_eval += 2; + + // Beginning a new iteration. + // Extrapolating by a factor of -1.0 through the face of the simplex + // across from the high point, i.e, reflect the simplex from the high point + for( j = 0 ; j < dim ; j++ ){ + sum = 0.0; + for( i = 0 ; i < mpts ; i++ ) + if (i != ihi) + sum += pts[i][j]; + + psum[j] = sum / dim; + } + + ytry = ModSimplex_(pts, y, psum, ihi, -1.0); + if( ytry <= y[ilo] ) { + // result better than best point + // so additional extrapolation by a factor of 2 + ytry = ModSimplex_(pts, y, psum, ihi, 2.0); + } + else if( ytry >= y[ihi] ) { + // result worse than the worst point + // so there is a lower intermediate point, + // i.e., do a one dimensional contraction + ysave = y[ihi]; + + ytry = ModSimplex_(pts, y, psum, ihi, -0.5); + if( ytry > y[ihi] ) { + // Can't get rid of the high point, + // try to contract around the best point + for( i = 0; i < mpts; i++ ) { + if( i != ilo ) { + for( j = 0; j < dim; j++ ) { + pts[i][j] = psum[j] = 0.5 * ( pts[i][j] + pts[ilo][j] ); + } + param_passed.CopyValues(psum); + y[i] = (*func_ptr_)(param_passed, data()); + } + } + num_func_eval += dim; + for( j = 0 ; j < dim ; j++ ){ + sum = 0.0; + for( i = 0 ; i < mpts ; i++ ) + if (i != ihi) + sum += pts[i][j]; + psum[j] = sum / dim; + } + } + } + else --num_func_eval; + } + fx_format_result(opt_module_, "func_evals", "%d", num_func_eval); + return; + } + + long double ModSimplex_(double **pts, long double *y, double *psum, + index_t ihi, float fac) { + + + index_t j, dim = dimension(); + long double ytry; + double *ptry; + Vector param_passed; + + param_passed.Init(dim); + ptry = (double*) malloc (dim * sizeof(double)); + for (j = 0; j < dim; j++) { + ptry[j] = psum[j] * (1 - fac) + pts[ihi][j] * fac; + } + param_passed.CopyValues(ptry); + ytry = (*func_ptr_)(param_passed, data()); + if (ytry < y[ihi]) { + y[ihi] = ytry; + for (j = 0; j < dim; j++) { + pts[ihi][j] = ptry[j]; + } + } + return ytry; + } + +}; + +/** + * An optimizer using the Quasi Newton method, + * also known as the variable metrics + * method. + * + * It does multivariate minimization of an + * objective function using only the function + * value and the gradients. + * + * Example use: + * + * @code + * double init_pt[d]; + * index_t number_of_iters; + * struct datanode *opt_module = fx_submodule(NULL,"QuasiNewton","opt_module"); + * Matrix data; + * index_t dim_param_space; + * + * ... + * QuasiNewton opt; + * opt.Init(obj_function, data, dim_param_space, opt_module); + * ... + * opt.Eval(init_pt, &number_of_function_evaluations); + * // init_pt contains the optimal point found + * @endcode + * + */ + +class QuasiNewton { + + private: + index_t dimension_; + Matrix data_; + long double (*func_ptr_)(Vector&, const Matrix&, Vector*); + datanode *opt_module_; + + public: + + QuasiNewton(){ + } + + ~QuasiNewton(){ + } + + void Init(long double (*fun)(Vector&, const Matrix&, Vector*), + Matrix& data, datanode *opt_module){ + + data_.Copy(data); + func_ptr_ = fun; + opt_module_ = opt_module; + dimension_ = fx_param_int_req(opt_module_, "param_space_dim"); + } + + const Matrix data() { + return data_; + } + + index_t dimension() { + return dimension_; + } + + void Eval(double *pt){ + + index_t n = dimension(), iters; + index_t i, its, MAXIMUM_ITERATIONS = fx_param_int(opt_module_,"MAX_ITERS",200); + long double temp_1, temp_2, temp_3, temp_4, f_previous, f_min, + maximum_step_length, sum = 0.0, sumdg, sumxi, temp, test; + Vector dgrad, grad, hdgrad, xi; + Vector pold, pnew; + Matrix hessian; + double EPSILON = fx_param_double(opt_module_, "EPSILON", 3.0e-8); + fx_format_param(opt_module_, "TOLERANCE", "%lf", 1.0e-5); + double TOLERANCE = fx_param_double_req(opt_module_, "TOLERANCE"); + double MAX_STEP_SIZE = fx_param_double(opt_module_, "MAX_STEP_SIZE", 100.0); + double g_tol = fx_param_double(opt_module_, "gtol", 1.0e-7); + + dgrad.Init(n); + grad.Init(n); + hdgrad.Init(n); + hessian.Init(n,n); + pnew.Init(n); + xi.Init(n); + pold.Copy(pt,n); + f_previous = (*func_ptr_)(pold, data(), &grad); + Vector tmp; + tmp.Init(n); + tmp.SetAll(1.0); + hessian.SetDiagonal(tmp); + la::ScaleOverwrite(-1.0, grad, &xi); + + sum = la::Dot(pold, pold); + double fmax; + if( sqrt(sum) > (float)n ) { + fmax = sqrt(sum); + } + else { + fmax = (float)n; + } + maximum_step_length = MAX_STEP_SIZE*fmax; + + for(its = 0; its < MAXIMUM_ITERATIONS; its++) { + + dgrad.CopyValues(grad); + LineSearch_(pold, f_previous, &grad, &xi, + &pnew, &f_min, maximum_step_length); + f_previous = f_min; + la::SubOverwrite(pold, pnew, &xi); + pold.CopyValues(pnew); + + for(i = 0; i < n; i++) { + pt[i] = pold.get(i); + } + + test = 0.0; + for(i = 0; i < n; i++){ + if(fabs(pold.get(i)) > 1.0) fmax = fabs(pold.get(i)); + else{ fmax = 1.0; } + temp = fabs(xi.get(i)) / fmax; + if(temp > test) test = temp; + } + if(test < TOLERANCE) { + iters = its; + fx_format_result(opt_module_, "iters", "%d", iters); + return; + } + + test = 0.0; + if(f_min > 1.0) temp_1 = f_min; + else{ temp_1 = 1.0; } + + for(i = 0; i < n; i++) { + if(fabs(pold.get(i)) > 1.0) fmax = pold.get(i); + else{ fmax = 1.0; } + + temp = fabs(grad.get(i))*fmax / temp_1; + if(temp > test) test = temp; + } + if(test < g_tol) { + iters = its; + fx_format_result(opt_module_, "iters", "%d", iters); + return; + } + + la::SubFrom(grad, &dgrad); + la::Scale(-1.0, &dgrad); + la::MulOverwrite(hessian,dgrad, &hdgrad); + + temp_2 = la::Dot(dgrad, xi); + temp_4 = la::Dot(dgrad, hdgrad); + sumdg = la::Dot(dgrad, dgrad); + sumxi = la::Dot(xi, xi); + + if (temp_2 > sqrt(EPSILON*sumdg*sumxi)) { + temp_2 = 1.0 / temp_2; + temp_3 = 1.0 / temp_4; + + la::ScaleOverwrite(temp_2, xi, &dgrad); + la::AddExpert((-1.0*temp_3), hdgrad, &dgrad); + + Matrix co, ro, tmp; + co.AliasColVector(xi); + ro.AliasRowVector(xi); + la::MulInit(co, ro, &tmp); + la::AddExpert(temp_2, tmp, &hessian); + + co.Destruct(); + ro.Destruct(); + tmp.Destruct(); + co.AliasColVector(hdgrad); + ro.AliasRowVector(hdgrad); + la::MulInit(co, ro, &tmp); + la::AddExpert((-1.0*temp_3), tmp, &hessian); + + co.Destruct(); + ro.Destruct(); + tmp.Destruct(); + co.AliasColVector(dgrad); + ro.AliasRowVector(dgrad); + la::MulInit(co, ro, &tmp); + la::AddExpert(temp_4, tmp, &hessian); + } + la::MulOverwrite(hessian, grad, &xi); + la::Scale((-1.0), &xi); + } + NOTIFY("Too many iterations in Quasi Newton\n"); + } + + void LineSearch_(Vector pold, long double fold, Vector *grad, + Vector *xi, Vector *pnew, long double *f_min, + long double maximum_step_length){ + + index_t i, n = dimension(); + long double a, step_length, previous_step_length = 0.0, + minimum_step_length, b, disc, previous_f_value = 0.0, + rhs1, rhs2, slope, sum, temp, test, temp_step_length, + MIN_DECREASE = 1.0e-4, TOLERANCE = 1.0e-7; + + sum = la::Dot(*xi, *xi); + sum = sqrt(sum); + if(sum > maximum_step_length) { + la::Scale((maximum_step_length/sum), xi); + } + slope = la::Dot(*grad, *xi); + if(slope >= 0.0){ + return; + } + test = 0.0; + for(i = 0; i < n; i++) { + double fmax; + fmax = (fabs(pold.get(i)) > 1.0 ? fabs(pold.get(i)) : 1.0); + temp = fabs((*xi).get(i)) / fmax; + if(temp > test) test = temp; + } + + minimum_step_length = TOLERANCE/test; + step_length = 1.0; + for(;;) { + + pnew->CopyValues(pold); + la::AddExpert(step_length, *xi, pnew); + + *f_min = (*func_ptr_)((*pnew), data(), grad); + if(step_length < minimum_step_length) { + pnew->CopyValues(pold); + return; + } + else if( *f_min <= fold + MIN_DECREASE*step_length*slope) { + return; + } + else { + if (step_length == 1.0) { + temp_step_length = -slope/(2.0*(*f_min - fold - slope)); + } + else { + rhs1 = *f_min - fold - step_length*slope; + rhs2 = previous_f_value - fold - previous_step_length*slope; + a = (rhs1 / (step_length*step_length) + - rhs2/(previous_step_length*previous_step_length)) + / (step_length-previous_step_length); + b = (-previous_step_length*rhs1/(step_length*step_length) + +step_length*rhs2/(previous_step_length*previous_step_length)) + / (step_length - previous_step_length); + if(a == 0.0) { + temp_step_length = -slope / (2.0*b); + } + else { + disc = b*b - 3.0*a*slope; + if(disc < 0.0) { + temp_step_length = 0.5*step_length; + } + else if (b <= 0.0) { + temp_step_length = (-b+sqrt(disc))/(3.0*a); + } + else { + temp_step_length = -slope / (b+sqrt(disc)); + } + } + if(temp_step_length > 0.5*step_length) { + temp_step_length = 0.5*step_length; + } + } + } + previous_step_length = step_length; + previous_f_value = *f_min; + step_length = (temp_step_length > 0.1*step_length + ? temp_step_length : 0.1*step_length); + } + } +}; + +#endif