Package com.imsl.math

Class NelderMead

java.lang.Object
com.imsl.math.NelderMead
All Implemented Interfaces:
Serializable, Cloneable
public class NelderMead extends Object implements Serializable, Cloneable
Minimizes a function of n variables with or without box constraints using a direct search polytope algorithm.

The algorithm uses a polytope method according to the algorithm of Nelder and Mead to find a minimum point of an unconstrained or constrained function of n variables. The constrained problem is of the form $$ \begin{array}{l} \min f(x_1,\ldots,x_n)\\ \text{subject to} \; l_i \le x_i \le u_i, \; i=1,\ldots,n. \end{array} $$ The algorithm is based on function comparison; no smoothness is assumed.

In the unconstrained case, the method iteratively updates a set of n + 1 points \(P = \{p_0,\ldots,p_n\}\) that form a simplex in \(\mathbf{R}^n\). At the start of each iteration, the two "worst" points \(p_{(n)} = {\rm{argmax}}_{p \in P}f(p)\), \(p_{(n-1)} = {\rm{argmax}}_{p \in P \setminus \{p_{(n)}\}}f(p)\) and the "best" point \(p_{(0)} = {\rm{argmin}}_{p \in P}f(p)\) are determined. The centroid of all but the worst point \(p_{(n)}\) $$ c := \frac{1}{n} \sum_{i \ne (n)} p_i $$ is used to compute a reflection point by the formula $$ x_r := c + \alpha \cdot (c-p_{(n)}). $$ Here, \(\alpha > 0\) denotes the reflection coefficient.

Depending on how the new reflection point \(x_r\) compares with the existing simplex vertices \(\{p_0,\ldots,p_n\}\), the iteration proceeds with one of the following three cases:

Case 1 (Reflection step)
If \(f(p_{(0)}) \le f(x_r) \le f(p_{(n-1)})\): Replace \(p_{(n)}\) by \(x_r\), \(P := (P \setminus \{p_{(n)}\}) \cup \{x_r\} \) and start the next iteration.

Case 2 (Expansion step)
If \(f(x_r) < f(p_{(0)})\): Compute $$ x_e := c + \beta \cdot (x_r - c) $$ where \(\beta > 1\) denotes the expansion coefficient.
Replace \(p_{(n)}\) by the better of \(x_r\) and \(x_e\), \(P := (P \setminus \{p_{(n)})\}) \cup \{{\rm{argmin}}\{f(x_r), f(x_e)\}\}\) and start the next iteration.

Case 3 (Contraction step)
If \(f(x_r) > f(p_{(n-1)})\): Compute the contraction point $$ x_c := \left\{ \begin{array}{ll} c + \gamma \cdot (p_{(n)} - c), & \text{if}\; f(x_r) \ge f(p_{(n)})\\ c + \gamma \cdot (x_r - c), & \text{if}\; f(x_r) < f(p_{(n)}) \end{array} \right. $$ where \(0 < \gamma < 1\) denotes the contraction coefficient.
If \(f(x_c) < \min\{f(p_{(n)}),f(x_r)\}\), replace \(p_{(n)}\) by \(x_c\), \(P := (P \setminus \{p_{(n)}\}) \cup \{x_c\} \).
Otherwise, shrink the polytope by moving the n worst points halfway towards the current best point \(p_{(0)}\): $$ \begin{array}{l} p_i := \frac{1}{2}(p_i + p_{(0)}), \, i = 0, \ldots, n\\ P := \{p_0, \ldots, p_n\} \end{array} $$ Start the next iteration.

The algorithm stops with the best point \(p_{(0)}\) returned as the optimum if at the start of an iteration one of the following stopping criteria is satisfied:
Criterion 1: $$ f(p_{(n)}) - f(p_{(0)}) \le \epsilon_f (1.0+|f(p_{(0)})|) $$ Criterion 2: $$ \frac{1}{n+1}\sum_{i=0}^{n}(f(p_i)-\bar{f})^2 \le \epsilon_f^2 $$ where \(\bar{f} = (\sum_{i=0}^{n}f(p_i))/(n+1)\) and \(\epsilon_f\) is a given tolerance.

The constrained version of the algorithm differs from the unconstrained version in two respects: The polytope is a complex with 2n vertices instead of n + 1, and infeasible points are projected onto the feasible set before their function value is evaluated.

For a complete description of the algorithms, see Nelder and Mead (1965), Box (1965) or Gill et al. (1981).

References
1. Box, M. J. (1965), A new method of constrained optimization and a comparison with other methods, The Computer Journal, Volume 8, Issue 1, 42–52.
2. Gill, P. E., W. Murray, and M. H. Wright (1981), Practical Optimization, Academic Press Inc. Limited, London.
3. Nelder, J. A., and R. Mead (1965), A simplex method for function minimization, The Computer Journal, Volume 7, Issue 4, 308-313.

See Also:

  • Constructor Details

    • NelderMead

      public NelderMead(NelderMead.Function function, int n)
      Constructor for unconstrained NelderMead.
      Parameters:
      function - a Function object, the user-supplied method to evaluate the objective function
      n - an int, the number of variables
      Throws:
      IllegalArgumentException - is thrown if n is less than one

    • NelderMead

      public NelderMead(NelderMead.Function function, double[] lowerBound, double[] upperBound)
      Constructor for constrained NelderMead.
      Parameters:
      function - a Function object, the user-supplied method to evaluate the objective function
      lowerBound - a double array containing the lower bounds on the variables. If variable i has no lower bound, set lowerBound[i] = -Double.MAX_VALUE.
      upperBound - a double array containing the upper bounds on the variables. If variable i has no upper bound, set upperBound[i] = Double.MAX_VALUE.
      Throws:
      IllegalArgumentException - is thrown if the bounds defined by lowerBound and upperBound are inconsistent.

  • Method Details

    • setNumberOfThreads

      public void setNumberOfThreads(int numberOfThreads)
      Sets the number of java.lang.Thread instances to be used for parallel processing.
      Parameters:
      numberOfThreads - an int specifying the number of java.lang.Thread instances to be used for parallel processing. If numberOfThreads is greater than 1, then interface method Function.f is evaluated in parallel and Function.f must be thread-safe. Otherwise, unexpected behavior can occur.

      Default: numberOfThreads = 1.

    • getNumberOfThreads

      public int getNumberOfThreads()
      Returns the number of java.lang.Thread instances used for parallel processing.
      Returns:
      an int containing the number of java.lang.Thread instances used for parallel processing.

    • solve

      public final void solve()
      Solves a minimization problem using a Nelder-Mead type algorithm.

    • setGuess

      public void setGuess(double[] initialGuess)
      Sets an initial guess of the solution.
      Parameters:
      initialGuess - a double array containing an initial guess for the solution of the problem.

      By default, initialGuess = 0.0.

    • setInitialComplex

      public void setInitialComplex(double[][] initialComplex)
      Defines the initial complex.
      Parameters:
      initialComplex - a double array of size numberOfVertices by n containing the numberOfVertices initial points of the complex. In the unconstrained case, numberOfVertices = n + 1, otherwise numberOfVertices = 2n, where n denotes the number of variables.

      By default, initialComplex is computed internally with a random number generator, using the initial guess as the first vertex of the complex.

      Note that after a solve call the internal flag indicating use of a user-defined initial complex is re-set to its default value, i.e. the use of a random complex. If the next solve call requires the initial complex just used, it has to be explicitly set again via method setInitialComplex.

    • getInitialComplex

      public double[][] getInitialComplex()
      Returns the initial complex.
      Returns:
      a double array of size numberOfVertices by n containing the vertices of the initial complex used in method solve. Here, n denotes the number of variables of the problem.

      Note that in the constrained case initial vertices that do not satisfy the box constraints are projected onto the feasible set. Therefore, the initial complex defined via method setInitialComplex needs not necessarily be identical with the complex returned by method getInitialComplex.

    • setTolerance

      public void setTolerance(double tolerance)
      Sets the convergence tolerance.
      Parameters:
      tolerance - a double containing the convergence tolerance

      By default, tolerance = 1.0e-8.

    • setReflectionCoefficient

      public void setReflectionCoefficient(double reflectionCoeff)
      Sets the value for the reflection coefficient.
      Parameters:
      reflectionCoeff - a double containing the reflection coefficient. reflectionCoeff must be greater than 0.

      By default, reflectionCoeff = 1.0.

    • setExpansionCoefficient

      public void setExpansionCoefficient(double expansionCoeff)
      Sets the value for the expansion coefficient.
      Parameters:
      expansionCoeff - a double containing the expansion coefficient. expansionCoeff must be greater than 1.

      By default, expansionCoeff = 2.0.

    • setContractionCoefficient

      public void setContractionCoefficient(double contractionCoeff)
      Sets the value for the contraction coefficient.
      Parameters:
      contractionCoeff - a double containing the contraction coefficient. contractionCoeff must be greater than 0 and less than 1.

      By default, contractionCoeff = 0.5.

    • setMaximumNumberOfFunctionEvaluations

      public void setMaximumNumberOfFunctionEvaluations(int maxEvalNumber)
      Sets the maximum number of allowed function iterations.
      Parameters:
      maxEvalNumber - an int scalar containing the maximum number of iterations.

      By default, maxEvalNumber = 300.

    • getObjectiveValue

      public double getObjectiveValue()
      Returns the value of the objective function at the computed solution.
      Returns:
      a double, the value of the objective function at the computed solution

    • getSolution

      public double[] getSolution()
      Returns the solution.
      Returns:
      a double array containing the computed solution

    • getNumberOfFunctionEvaluations

      public int getNumberOfFunctionEvaluations()
      The number of function evaluations needed.
      Returns:
      an int, the number of function evaluations needed to compute the solution

    • getCentroidDistance

      public double getCentroidDistance()
      Returns the average distance from the final complex vertices to their centroid.
      Returns:
      a double, a measure for the flatness of the function near the computed minimum

    • setRandomObject

      public final void setRandomObject(Random r)
      Sets the random object to be used in the computation of the vertices of the complex.
      Parameters:
      r - a Random object to be used in the computation of the complex vertices.

      Specifying a seed for the Random object can produce repeatable/deterministic output.

    • getRandomObject

      public Random getRandomObject()
      Returns the random object being used in the computation of the vertices of the initial complex.
      Returns:
      a Random object being used for the complex vertices