From hrz.uni-dortmund.de!a.christmann Wed Mar 05 16:05:31 0800 1997 remote from utstat Received: by utstat; Wed Mar 5 09:58 EST 1997 Received: from nvl1.hrz.uni-dortmund.de by nx1.hrz.uni-dortmund.de with SMTP (PP); Wed, 5 Mar 1997 15:57:40 +0100 Received: from UNIDO_HRZ1/TEMPQ by nvl1.hrz.uni-dortmund.de (Mercury 1.21); 5 Mar 97 16:05:04 +0100 Received: from TEMPQ by UNIDO_HRZ1 (Mercury 1.21); 5 Mar 97 16:04:37 +0100 Received: from dynam142.HRZ.Uni-Dortmund.DE by nvl1.hrz.uni-dortmund.de (Mercury 1.21); 5 Mar 97 16:04:32 +0100 Message-ID: <331E0A4B.629B@hrz.uni-dortmund.de> Date: Wed, 05 Mar 1997 16:05:31 -0800 From: "Dr. A. Christmann" Organization: HRZ X-Mailer: Mozilla 2.02 (Win16; I) MIME-Version: 1.0 To: s-news Subject: Summary : NELDER-MEAD SIMPLEX ALGORITHM Content-Type: text/plain; charset=us-ascii Content-Transfer-Encoding: 7bit Hello, I like to thank all persons who responded on my question. ----------------------------------------------------------------- "Has anyone implemented the Nelder-Mead Simplex minimization procedure in S-Plus ? (This minimization procedure is published as FORTRAN code in ALGORITHM AS 47 APPL. STATIST. (1971) VOL. 20, P. 338 and is also alvailable in StatLib.)" ----------------------------------------------------------------- There are many persons who send me an answer. From my point of view there are two different types of answers. Method a: pure S-Plus code Method b: code in Fortran or C, then use dyn.load or dll.load Of course, using code in Fortran or C is faster in most cases. I work with S-Plus under Windows 3.1 and I have a FORTRAN 77 version published in Applied Statistics and a FORTRAN 90 version from Jeff Breiwick (MINUM) of the Nelder-Mead simplex algorithm. As I do not yet know how to construct a DLL from a NAG FORTRAN 90 compiler, which seems to be necessary to use method b under Windows 3.1, I currently tend to use pure S-Plus code. Or is there another way to use method b under Windows 3.1 ? Therefore, I tried the S-Plus code developed by Bill Clark and it seems to work well for my problem. Thank you, Andreas Christmann Here is a summary of all answers in the sequence I got them. ------------------------------------------------------------ Yingwei Peng: ------------ What do you mean here by "in Splus"? If you mean implementing the procedure purely in S language, it will be too slow for any practical use. If you mean implementing it in other advanced languages and then loading it into Splus, you can easily compile the FORTRAN code and load it in Splus. I have used a C code by @BOOK{70, AUTHOR = "William H. Press and Saul A. Teukolsky and William T. Vetterling and Brian P. Flannery", TITLE = "Numerical recipes in C : the art of scientific computing", PUBLISHER = "Cambridge University Press", YEAR = 1992, ADDRESS = "New York" } It works well. A Fortran code from the same authors is also available. However, I used the C code only as a part of my C programs loaded into S. So I don't have a special Splus interface for the code only. Let me know if you like to have a copy of the code. ------------------------------------------------------------------------------------- Daniel Heitjan: --------------- I have implemented the Nelder-Mead algorithm (fun.amoeba), the Powell algorithm (fun.powell), and the Davidon-Fletcher- Powell algorithm (fun.dfpmin), in Splus. All implementations are based on the _Numerical Recipes_ versions. They are available by anonymous ftp to 156.111.36.191; cd to pub/dheitjan. There you will find three Splus files: the optimization routines (optfcn.q), some test optimization problems (testfcns.q), and code for running the tests (testalg.q). Some of this code is several years old, and I haven't tested it recently, so you will probably have to do a bit of detective work (including changing file names in source() statements) to make it go. But I have used this Nelder-Mead algorithm many times, and it does work well. I can provide some limited help if you experience problems. ---------------------------------------------------------------------------------- Bill Clark: ----------- I wrote and tested the version below. It follows Nelder and Mead's terms and strategy pretty closely. Calling the routine with print.level=2 will track the search. You may want to alter the tolerance parameter (tol) used in the convergence test. nldrmd <- function(fobj, startvec, print.level=0) # Function minimization by the Nelder-Mead simplex method. This is an # inefficient minimizer but is useful for obtaining a decent starting # point and scale information for one of the efficient methods. # Ref.: Press (Numerical Mehtods in C) p. 305. { # Standard tuning parameters. alpha <- 1.0 # alpha >= 1 beta <- 0.5 # 0 < beta < 1 gamma <- 2.0 # gamma > alpha tol <- 1E-6 # Build initial simplex by scaling single elements of starting vector by # a factor exp(+/-lambda). if (print.level==2) cat('Initializing simplex...\n') lambda <- 2 npar <- length(startvec) nverts <- npar + 1 verts <- matrix(startvec,nverts,npar,byrow=T) fvals <- rep(0,nverts) for (iv in 1:npar) repeat { lamb <- lambda { verts[iv,iv] <- exp(lamb) * startvec[iv] if (verts[iv,iv]==0) verts[iv,iv] <- 1 fvals[iv] <- fobj(verts[iv,]) if (is.finite(fvals[iv])) break verts[iv,iv] <- exp(-lamb) * verts[iv,iv] fvals[iv] <- fobj(verts[iv]) if (is.finite(fvals[iv])) break lamb <- lamb/2 } } verts[nverts,] <- startvec fvals[nverts] <- fobj(startvec) scatter.init <- apply(verts,2,sd) centroid <- startvec # Start iteration. iter <- 0 repeat { # Check for convergence. if (print.level==2) ask('Press Return to continue') iter <- iter + 1 scatter <- apply(verts,2,sd) conv.test <- max(scatter/scatter.init) if (conv.test < tol) break if (print.level > 0) { fcen <- fobj(centroid) cat('\nIteration',iter,' ... convergence test =',ff(conv.test), 'fval =',ff(fcen,9),'\n') } # Make a move. First try projecting across centroid from high point. fmax <- max(fvals) iv.fmax <- min((1:nverts)[fvals==fmax]) fnext <- max(fvals[(1:nverts) != iv.fmax]) iv.fnext <- min((1:nverts)[fvals==fnext]) fmin <- min(fvals) iv.fmin <- min((1:nverts)[fvals==fmin]) # The centroid referred to here is the centroid of the face # across from the point where the function is largest. centroid <- if (npar==1) verts[iv.fmin] else apply(verts[(1:nverts) != iv.fmax,],2,mean) proj.vert <- centroid + alpha * (centroid - verts[iv.fmax,]) proj.fval <- fobj(proj.vert) if (print.level == 2) { cat('Vertices:\n') for (iv in 1:nverts) cat(iv,ff(verts[iv,]),'\n') cat('Function values:\n',ff(fvals),'\n') cat('Projected vertex:\n',ff(proj.vert),'\n') cat('Projected function value:',ff(proj.fval),'\n') } # Make sure function values are different. if (all(fvals==fmax)) { cat(paste('nldrmd error: Function values same at all vertices.', 'Change startvec.\n')) return(NA) } # Decide what to do based on projected function value. if (proj.fval >= fmin & proj.fval <= fnext) # Middling result. Take it. { if (print.level==2) cat('Taking projection.\n') verts[iv.fmax,] <- proj.vert fvals[iv.fmax] <- proj.fval next } if (proj.fval < fmin) # Better than the minimum. Perform an expansion. { if (print.level==2) cat('Performing expansion.\n') exp.vert <- centroid + gamma * (centroid - verts[iv.fmax,]) exp.fval <- fobj(exp.vert) verts[iv.fmax,] <- exp.vert fvals[iv.fmax] <- exp.fval next } if (proj.fval > fnext) # Worse than second highest. Try a one-dimensional contraction. { cont.vert <- beta * centroid + (1 - beta) * verts[iv.fmax,] cont.fval <- fobj(cont.vert) if (cont.fval < fmax) { if (print.level==2) cat('Performing 1D contraction.\n') verts[iv.fmax,] <- cont.vert fvals[iv.fmax] <- cont.fval next } # Didn't work. Shrink toward best point. if (print.level==2) cat('Performing wholesale contraction.\n') for (iv in 1:nverts) verts[iv,] <- beta * verts[iv.fmin,] + (1 - beta) * verts[iv,] } # Done with update. } # Convergence test succeeded. return(verts[iv.fmin,]) } ----------------------------------------------------------------------------------- Brian Dawkins: ------------- I implemented this with an interface to C-code several years ago. It is in statlib under the name recipes together with several other bits and pieces. I've had a number of communications over the years from various users and it does seem to work reasonably as far as I know. ------------------------------------------------------------------------------------ Peter Perkins: -------------- i know i hate this kind of answer, but: i've used this algorithm (in C, not from S) for a number of different problems and as far as i can tell its only saving grace is that it seems to be more robust to steep slopes and discontinuities than poorly implemented quasi-newton methods (such as those in Numerical Recipes). the fact that it doesn't require derivatives is a red herring - neither do modern quasi-newton algorithms (but they work faster if you can supply it). the optimizer in S, nlmin(), is really a pretty good one and works much much faster than the simplex search in all of the problems that i have used them. there are better optimizers as well - i use NPSOL outside of S, but nlmin() (based on ATT's PORT routines) is pretty good. ---------------------------------------------------------------------------------- Daniel Heitjan: --------------- Several years ago I implemented a version of the Nelder-Mead simplex algorithm, based on the _Numerical Recipes_ fortran routine amoeba, entirely in Splus. There is no doubt that a Fortran or C implementation would be faster. But there are certain advantages to an all-Splus implementation: It is easier to read and modify, it can be used to optimize any Splus function, it makes it easier to track down problems in the convergence, and it is more portable (e.g., between platforms). If I can do the programming quickly, as I normally can in Splus, I don't mind if I don't get the answer right away. In short, I have found lots of uses for my Splus amoeba. Nevertheless, I always dreamed of an Splus wrapper for the fortran amoeba that would work the same as my Splus amoeba--i.e., optimize any Splus function with any kind of data. But I couldn't figure out how to get the fortran to recognize the function and the data. Has anybody written such a wrapper? ------------------------------------------------------------------------------------- Brian Dawkins: -------------- Several years ago I wrote just such a wrapper. It is in statlib under the rubric of recipes. I've had a little correspondence over the years about this and I'm reasonably confident it works okay. -------------------------------------------------------------------------------------- Martin Maechler: --------------- I would suggest to use the C (instead of fortran) amoeba code. (Numerical recipes in C). Then it is possible to pass S functions as you know from the 'zero' example in the blue book and the "S-plus Programming" manual. --------------------------------------------------------------------------------------- Yingwei Peng : -------------- > Martin Maechler wrote, > > I would suggest to use the C (instead of fortran) amoeba code. > (Numerical recipes in C). Yes, this is exactly what I said in the message to Dr. A. Christmann. > Then it is possible to pass S functions as you know from the 'zero' example > in the blue book and the "S-plus Programming" manual. Yes, it is possible using call_S() function provided by Splus. But it takes even longer time for programs to switch between Splus and C. From a point view of portability, I would prefer either pure S code (if the problem is not large) or C code and try to minimize the times of calling between each other. ---------------------------------------------- A.Christmann@hrz.uni-dortmund.de Andreas Christmann /////// Universitaet Dortmund U N I D O /// Hochschulrechenzentrum ______/////// Wissenschaftl. Anwendungen \_\_\_\///// D-44221 Dortmund \_\_\_\/// Tel. 049-231-755-2763 \_\_\_\/ Fax. -2731