Nelder-Mead Algorithm with custom initial simplex

Configure your own initial simplex and transformation factors in Nelder-Mead algorithm

TL, DR:

Here is a simple implementation:

fun = @(x)100*(x(2) - x(1)^2)^2 + (1 - x(1))^2;

NM = nm(fun,rand(100,2)); % initialize a simplex search project, start from a simplex of 100 points

NM.save_complete_data = false; % do not save data to log unless you want to visualize the iteration

NM.iter_limit = 1000; % limit to 1000 iterations

NM.auto_run();

NM.curr_x % returns simplex vertex of lowest function value

The motivation to write this project is that there MATLAB fminsearch does not allow user to (1) configure custom initial simplex, and (2) configure custom contract factor, etc. This project is just a realization of the algorithm described on this page:

http://www.scholarpedia.org/article/Nelder-Mead_algorithm

In other words, the actual implementation might not be exactly as MATLAB.

Please cite the site below:

Saša Singer and John Nelder (2009), Scholarpedia, 4(7):2928.

A brief introduction:

NM = nm(fun,x0) initializes a nelder-mead simplex search object with function handle fun and initial simplex x0;

How to configure initial simplex x0:

x0 is a matrix with f_len columns and sim_lenrows; In a default Nelder-Mead implementation, sim_len should be f_len + 1. You can define your own custom initial simplex. Here are a couple of things to keep in mind:

sim_len >= f_len + 1 is not required. And even if the criterion is met, it does not necessarily guarantee that input is a sufficient simplex. A function is_simplex which uses built-in function rref is used to determine if an input simplex is sufficient. Otherwise, it shows a warning message. Some may want to perform search on an insufficient simplex.

You can also insert a single vertex of a simplex and get_simplex will be invoked to get a simplex just as the implementation in MATLAB. Make sure it’s a row vector not column vector.

A more useful approach for real world problems in my opinion is starting from a redundant simplex which has (much) more vertices than minimum requirement and make sure that potential optimal solution is located inside the range defined by the vertices. To get a such initial simplex, you can try array2simplexfunction. The function takes n arrays (n = f_len) as initial points and gives a simplex that can be directly sent to nm constructor. Note that each array does not have to be of the same length. In practice in some dimensions the target function always has single global minimum (so you only need 2 values in such dimension) and in some dimensions the function becomes complex. You need to configure based on your understanding of the target function handle.

Optimization of the algorithm: the most time-consuming steps are usually function handle evaluations and sorting of large simplexes. Fortunately, these two can be accelerated. In most cases (excluding initialization and shrinking) only 1~3 function queries are needed in an iteration. If you check the iteration log, shrink is very rare actually. Excluding shrinking, it’s only necessary to replace one vertex and there is no need to sort the entire simplex. This was achieved using binary_insertwhich performs a binary search to find the right place of the new simplex.

Assignment of simplex operation parameters

There are 4 parameters (expand, contract, shrink, reflect) that can be user defined before iteration and the range of these parameters were commented. One parameter that needs to be considered in my opinion is shrink factor as when the number of vertices is very large the simplex tends to converge to one after shrink operations. So you might want to use a more conservative shrink factor (0.85 rather than 0.5, for example).

NM.auto_run() performs automatic iteration until ending criteria is met. It basically just calls NM.single_iteration()repeatedly. If you want to monitor the optimization real-time, you can manually call this function.

Important properties user can check out during iteration:

NM.curr_x: row vector of f_len, which is the vertex that gives best (lowest) function value;

NM.curr_fval: current lowest function value;

NM.simplex: current simplex;

NM.simplex_values: current function values at each vertex;

I have been using the uploaded codes for a while and fixed numerous bugs but there could still be (a lot) more issues and I did not seriously think about access control, etc. There are some sample codes included. Good luck!