Lecture 17

Multidimensional Optimization: The Simplex Method Lecture 17 Biostatistics 615/815 Previous Topic: One-Dimensional Optimization z Bracketing z Golden Search z Quadratic Approximation Bracketing z Find 3 points such that ? a < b < c ? f(b) < f(a) and f(b) < f(c) z Locate minimum by gradually trimming bracketing interval z Bracketing provides additional confidence in result The Golden Ratio A B CX A B C Bracketing Triplet 0.38196 0.38196 New Point Parabolic Interpolation For well behaved functions, faster than Golden Search Today: Multidimensional Optimization z Illustrate the method of Nelder and Mead ? Simplex Method ? Nicknamed "Amoeba" z Simple and, in practice, quite robust ? Counter examples are known z Discuss other standard methods C Utility Functions: Allocating Vectors z Ease allocation of vectors. z Peppered through today's examples double * alloc_vector(int cols) { return (double *) malloc(sizeof(double) * cols); } void free_vector(double * vector, int cols) { free(vector); } C Utility Functions: Allocating Matrices double ** alloc_matrix(int rows, int cols) { int i; double ** matrix = (double **) malloc(sizeof(double *) * rows); for (i = 0; i < rows; i++) matrix[i] = alloc_vector(cols); return matrix; } void free_matrix(double ** matrix, int rows, int cols) { int i; for (i = 0; i < rows; i++) free_vector(matrix[i], cols); free(matrix); } The Simplex Method z Calculate likelihoods at simplex vertices ? Geometric shape with k+1 corners ? E.g. a triangle in k = 2 dimensions z Simplex crawls ? Towards minimum ? Away from maximum z Probably the most widely used optimization method A Simplex in Two Dimensions z Evaluate function at vertices z Note: ? The highest (worst) point ? The next highest point ? The lowest (best) point z Intuition: ? Move away from high point, towards low point x 0 x 1 x 2 C Code: Creating A Simplex double ** make_simplex(double * point, int dim) { int i, j; double ** simplex = alloc_matrix(dim + 1, dim); for (int i = 0; i < dim + 1; i++) for (int j = 0; j < dim; j++) simplex[i][j] = point[j]; for (int i = 0; i < dim; i++) simplex[i][i] += 1.0; return simplex; } C Code: Evaluating Function at Vertices z This function is very simple ? This is a good thing! ? Making each function almost trivial makes debugging easy void evaluate_simplex (double ** simplex, int dim, double * fx, double (* func)(double *, int)) { for (int i = 0; i < dim + 1; i++) fx[i] = (*func)(simplex[i], dim); } C Code: Finding Extremes void simplex_extremes(double *fx, int dim, int & ihi, int & ilo, int & inhi) { int i; if (fx[0] > fx[1]) { ihi = 0; ilo = inhi = 1; } else { ihi = 1; ilo = inhi = 0; } for (i = 2; i < dim + 1; i++) if (fx[i] <= fx[ilo]) ilo = i; else if (fx[i] > fx[ihi]) { inhi = ihi; ihi = i; } else if (fx[i] > fx[inhi]) inhi = i; } Direction for Optimization x 0 x 1 x 2 Average of all points, excluding worst point mid Line through worst point and average of other points C Code: Direction for Optimization void simplex_bearings(double ** simplex, int dim, double * midpoint, double * line, int ihi) { int i, j; for (j = 0; j < dim; j++) midpoint[j] = 0.0; for (i = 0; i < dim + 1; i++) if (i != ihi) for (j = 0; j < dim; j++) midpoint[j] += simplex[i][j]; for (j = 0; j < dim; j++) { midpoint[j] /= dim; line[j] = simplex[ihi][j] - midpoint[j]; } } Reflection x 0 x 1 x 2 mid x' This is the default new trial point Reflection and Expansion x 0 x 1 x 2 mid x' If reflection results in new minimum? x'' Move further along minimization direction Contraction (One Dimension) x 0 x 1 x 2 mid x' If x' is still the worst point? x'' Try a smaller step C Code: Updating The Simplex int update_simplex(double * point, int dim, double & fmax, double * midpoint, double * line, double scale, double (* func)(double *, int)) { int i, update = 0; double * next = alloc_vector(dim), fx; for (i = 0; i < dim; i++) next[i] = midpoint[i] + scale * line[i]; fx = (*func)(next, dim); if (fx < fmax) { for (i = 0; i < dim; i++) point[i] = next[i]; fmax = fx; update = 1; } free_vector(next, dim); return update; } Contraction ? x 0 x 1 x 2 mid If a simple contraction doesn't improve things, then try moving all points towards the current minimum "passing through the eye of a needle" C Code: Contracting The Simplex void contract_simplex(double ** simplex, int dim, double * fx, int ilo, double (*func)(double *, int)) { int i, j; for (int i = 0; i < dim + 1; i++) if (i != ilo) { for (int j = 0; j < dim; j++) simplex[i][j] = (simplex[ilo][j]+simplex[i][j])*0.5; fx[i] = (*func)(simplex[i], dim); } } Summary: The Simplex Method high low Original Simplex reflection reflection and expansion contraction multiple contraction C Code: Minimization Routine (Part I) z Declares local variables and allocates memory double amoeba(double *point, int dim, double (*func)(double *, int), double tol) { int ihi, ilo, inhi, j; double fmin; double * fx = alloc_vector(dim + 1); double * midpoint = alloc_vector(dim); double * line = alloc_vector(dim); double ** simplex = make_simplex(point, dim); evaluate_simplex(simplex, dim, fx, func); C Code: Minimization Route (Part II) while (true) { simplex_extremes(fx, dim, ihi, ilo, inhi); simplex_bearings(simplex, dim, midpoint, line, ihi); if (check_tol(fx[ihi], fx[ilo], tol)) break; update_simplex(simplex[ihi], dim, fx[ihi], midpoint, line, -1.0, func); if (fx[ihi] < fx[ilo]) update_simplex(simplex[ihi], dim, fx[ihi], midpoint, line, -2.0, func); else if (fx[ihi] >= fx[inhi]) if (!update_simplex(simplex[ihi], dim, fx[ihi], midpoint, line, 0.5, func)) contract_simplex(simplex, dim, fx, ilo, func); } C Code: Minimization Routine (Part III) z Store the result and free memory for (j = 0; j < dim; j++) point[j] = simplex[ilo][j]; fmin = fx[ilo]; free_vector(fx, dim); free_vector(midpoint, dim); free_vector(line, dim); free_matrix(simplex, dim + 1, dim); return fmin; } C Code: Checking Convergence #include #define ZEPS 1e-10 int check_tol(double fmax, double fmin, double ftol) { double delta = fabs(fmax - fmin); double accuracy =(fabs(fmax) + fabs(fmin)) * ftol; return (delta <(accuracy + ZEPS)); } amoeba() z A general purpose minimization routine ? Works in multiple dimensions ? Uses only function evaluations ? Does not require derivatives z Typical usage: ?my_func(double * x, int n) { ? } ?amoeba(point, dim, my_func, 1e-7); Example Application Old Faithful Eruptions (n = 272) Old Faithful Eruptions Duration (mins) Fr eq ue nc y 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 0 5 10 15 2 0 Fitting a Normal Distribution z Fit two parameters ? Mean ? Variance z Requires ~165 likelihood evaluations ? Mean = 3.4878 ? Variance = 1.2979 ? Maximum log-likelihood = -421.42 Nice fit, eh? 123456 0 . 0 5 0. 10 0. 15 0. 20 0. 25 0. 30 0 . 3 5 Fitted Distribution Duration (mins) De n s i t y Old Faithful Eruptions Duration (mins) Fr eq ue n c y 123456 0 5 10 15 2 0 A Mixture of Two Normals z Fit 5 parameters ? Proportion in the first component ? Two means ? Two variances z Required about ~700 evaluations ? First component contributes 0.34841 of mixture ? Means are 2.0186 and 4.2734 ? Variances are 0.055517 and 0.19102 ? Maximum log-likelihood = -276.36 Two Components Old Faithful Eruptions Duration (mins) Fr eq ue n c y 123456 0 5 10 15 2 0 123456 0 . 0 0 .1 0 . 2 0 .3 0 . 4 0 . 5 0 . 6 Fitted Distribution Duration (mins) De n s i t y A Mixture of Three Normals z Fit 8 parameters ? Proportion in the first two components ? Three means ? Three variances z Required about ~1400 evaluations ? Did not always converge! z One of the best solutions ? ? Components contributing .339, 0.512 and 0.149 ? Component means are 2.002, 4.401 and 3.727 ? Variances are 0.0455, 0.106, 0.2959 ? Maximum log-likelihood = -267.89 Three Components Old Faithful Eruptions Duration (mins) Fr eq ue n c y 123456 0 5 10 15 2 0 123456 0. 0 0 . 1 0. 2 0 . 3 0. 4 0 . 5 0 . 6 0 . 7 Fitted Distribution Duration (mins) De n s i t y Tricky Minimization Questions z Fitting variables that are constrained ? Proportions vary between 0 and 1 ? Variances must be positive z Selecting the number of components z Checking convergence Improvements to amoeba() z Different scaling along each dimension ? If parameters have different impact on the likelihood z Track total function evaluations ? Avoid getting stuck if function does not cooperate z Rotate simplex ? If the current simplex is leading to slow improvement optim() Function in R z optim(point, function, method) ? Point ? starting point for minimization ? Function that accepts point as argument ? Method can be ? "Nelder-Mead" for simplex method (default) ? "BFGS", "CG" and other options use gradient One parameter at a time z Simple but inefficient approach z Consider ? Parameters ? = (? 1 , ? 2 , ?, ? k ) ? Function f (?) z Maximize ? with respect to each ? i in turn ? Cycle through parameters The Inefficiency? ? 1 ? 2 Steepest Descent z Consider ? Parameters ? = (? 1 , ? 2 , ?, ? k ) ? Function f(?; x) z Score vector ? Find maximum along ? + ?S ? ? ? ? ? ? ? ? == k d fd d fd d fd S ??? ln ...,, lnln 1 Still inefficient? Consecutive steps are still perpendicular! Multidimensional Minimization z Typically, sophisticated methods will? z Use derivatives ? May be calculated numerically. How? z Select a direction for minimization, using: ? Weighted average of previous directions ? Current gradient ? Avoid right angle turns Recommended Reading z Numerical Recipes in C (or C++, or Fortran) ? Press, Teukolsky, Vetterling, Flannery ? Chapter 10.4 z Clear description of Simplex Method ? Other sub-chapters illustrate more sophisticated methods z Online at ? http://www.numerical-recipes.com/ Goncalo Abecasis Microsoft PowerPoint - 615.17 -- Simplex Optimization