Get started today!

Good to have you back!
If you've signed in to StudyBlue with Facebook in the past, please do that again.

- StudyBlue
- Michigan
- University of Michigan - Ann Arbor
- Biostatistics And Medical Informatics
- Biostatistics And Medical Informatics 615
- Abecasis
- Lecture 17

Anonymous

Advertisement

Advertisement

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

"StudyBlue is great for studying. I love the study guides, flashcards and quizzes. So extremely helpful for all of my classes!"

Alice , Arizona State University"I'm a student using StudyBlue, and I can 100% say that it helps me so much. Study materials for almost every subject in school are available in StudyBlue. It is so helpful for my education!"

Tim , University of Florida"StudyBlue provides way more features than other studying apps, and thus allows me to learn very quickly!??I actually feel much more comfortable taking my exams after I study with this app. It's amazing!"

Jennifer , Rutgers University"I love flashcards but carrying around physical flashcards is cumbersome and simply outdated. StudyBlue is exactly what I was looking for!"

Justin , LSU