You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by er...@apache.org on 2021/08/07 13:36:05 UTC

[commons-math] 08/11: MATH-1622: Hybridation of simplex-based search with simulated annealing (SA).

This is an automated email from the ASF dual-hosted git repository.

erans pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-math.git

commit 50ca244b38e130923b5426092895c587e207f4a2
Author: Gilles Sadowski <gi...@gmail.com>
AuthorDate: Mon Aug 2 00:19:03 2021 +0200

    MATH-1622: Hybridation of simplex-based search with simulated annealing (SA).
    
    Class "SimulatedAnnealing" provides user-defined input ("OptimizationData").
    Behaviour of "SimplexOptimizer" depends on whether SA is activated (if yes,
    each "iteration" is further split in "epoch" stages).
    The "Simplex.TransformFactory" interface is changed to allow passing the SA
    acceptance test ("Metropolis" rule).
    
    The "SimplexOptimizer" is furter refactored (cf. MATH-1614) so that the
    steps (in method "doOptimize") clearly stand out as
     1. Setup and iteration (optionally with simulated annealing).
     2. Simplex transformation rule.
     3. Convergence check.
     4. "Best list" search.
    
    This commit also contains the following modifications:
     * "TestFunction": Override "toString" (for display in unit tests).
     * "Simplex"
       - Self-documenting names for factories.
       - Utility methods: "centroid", "asList", "replaceLast".
       - Unit tests.
     * "MultiDirectionalTransform":
       - Use common utility methods.
       - Fix missing sort in "transform" method.
       - Check validity of constructor parameters.
       - Parameter names to match the equivalent in "NelderMeadMeadTransform".
     * Adapt unit test classes to the updated API of class "Simplex".
     * "PointValuePair":
       - "final".
       - Implements "hashCode".
---
 .../commons/math4/legacy/optim/PointValuePair.java |   8 +-
 .../optim/nonlinear/scalar/SimulatedAnnealing.java | 195 +++++++++++++
 .../scalar/noderiv/HedarFukushimaTransform.java    | 210 ++++++++++++++
 .../scalar/noderiv/MultiDirectionalTransform.java  |  76 +++--
 .../scalar/noderiv/NelderMeadTransform.java        |  28 +-
 .../optim/nonlinear/scalar/noderiv/Simplex.java    | 261 ++++++++++++------
 .../nonlinear/scalar/noderiv/SimplexOptimizer.java | 305 +++++++++++++++++++--
 .../MultivariateFunctionPenaltyAdapterTest.java    |  10 +-
 .../optim/nonlinear/scalar/TestFunction.java       |  38 ++-
 .../SimplexOptimizerHedarFukushimaTest.java        | 161 +++++++++++
 .../SimplexOptimizerMultiDirectionalTest.java      |  21 +-
 .../noderiv/SimplexOptimizerNelderMeadTest.java    |  26 +-
 .../nonlinear/scalar/noderiv/SimplexTest.java      | 199 ++++++++++++++
 13 files changed, 1345 insertions(+), 193 deletions(-)

diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/PointValuePair.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/PointValuePair.java
index 35f883c..2d1a62c 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/PointValuePair.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/PointValuePair.java
@@ -29,7 +29,8 @@ import org.apache.commons.math4.legacy.core.Pair;
  * @see org.apache.commons.math4.legacy.analysis.MultivariateFunction
  * @since 3.0
  */
-public class PointValuePair extends Pair<double[], Double> implements Serializable {
+public final class PointValuePair extends Pair<double[], Double>
+    implements Serializable {
     /** Serializable UID. */
     private static final long serialVersionUID = 20120513L;
 
@@ -96,6 +97,11 @@ public class PointValuePair extends Pair<double[], Double> implements Serializab
     }
 
     @Override
+    public int hashCode() {
+        return Arrays.hashCode(getKey()) + 31 * Double.hashCode(getValue());
+    }
+
+    @Override
     public String toString() {
         return "[" + Arrays.toString(getPointRef()) + ", " + getValue() + "]";
     }
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/SimulatedAnnealing.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/SimulatedAnnealing.java
new file mode 100644
index 0000000..3ce2c7d
--- /dev/null
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/SimulatedAnnealing.java
@@ -0,0 +1,195 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math4.legacy.optim.nonlinear.scalar;
+
+import java.util.function.BiFunction;
+import java.util.function.DoublePredicate;
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.math4.legacy.stat.descriptive.moment.StandardDeviation;
+import org.apache.commons.math4.legacy.optim.OptimizationData;
+import org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv.Simplex;
+
+/**
+ * Simulated annealing setup.
+ *
+ * @since 4.0
+ */
+public class SimulatedAnnealing implements OptimizationData {
+    /** Number of iterations at fixed temperature. */
+    private final int epochDuration;
+    /** Initial acceptance probability. */
+    private final double startProbability;
+    /** Final acceptance probability. */
+    private final double endProbability;
+    /** Cooling function. */
+    private final CoolingSchedule coolingSchedule;
+    /** RNG. */
+    private final UniformRandomProvider rng;
+
+    /**
+     * @param epoch Number of iterations performed at fixed temperature.
+     * @param startProb Initial acceptance probablility.
+     * @param endProb Final acceptance probablility.
+     * @param cooling Computes the temperature as a function of the initial
+     * temperature and the epoch.
+     * It is called for computing a new temperature after each cycle of
+     * {@code epoch} iterations.
+     * Simulated annealing <em>assumes</em> that the function decreases
+     * monotically wrt the epoch (cf. {@link CoolingSchedule#decreasingExponential(double)
+     * provided implementation}).
+     * @param random Random number generator.
+     * @throws IllegalArgumentException if {@code epoch < 1} or
+     * {@code startProb} or {@code endProb} is outside the {@code [0, 1]}
+     * interval.
+     */
+    public SimulatedAnnealing(int epoch,
+                              double startProb,
+                              double endProb,
+                              CoolingSchedule cooling,
+                              UniformRandomProvider random) {
+        if (epoch < 1) {
+            throw new IllegalArgumentException("Epoch out of range: " +
+                                               epoch);
+        }
+        if (startProb < 0 ||
+            startProb > 1) {
+            throw new IllegalArgumentException("Initial acceptance probability out of range: " +
+                                               startProb);
+        }
+        if (endProb < 0 ||
+            endProb > 1) {
+            throw new IllegalArgumentException("Final acceptance probability out of range: " +
+                                               endProb);
+        }
+        if (endProb >= startProb) {
+            throw new IllegalArgumentException("Final probability larger than initial probability");
+        }
+
+        epochDuration = epoch;
+        startProbability = startProb;
+        endProbability = endProb;
+        coolingSchedule = cooling;
+        rng = random;
+    }
+
+    /**
+     * @return the epoch duration.
+     */
+    public int getEpochDuration() {
+        return epochDuration;
+    }
+
+    /**
+     * @return the acceptance probability at the beginning of the SA process.
+     */
+    public double getStartProbability() {
+        return startProbability;
+    }
+
+    /**
+     * @return the acceptance probability at the end of the SA process.
+     */
+    public double getEndProbability() {
+        return endProbability;
+    }
+
+    /**
+     * @return the cooling schedule.
+     */
+    public CoolingSchedule getCoolingSchedule() {
+        return coolingSchedule;
+    }
+
+    /**
+     * Specifies the cooling schedule.
+     * It computes the current temperature as a function of two arguments:
+     * <ol>
+     *  <li>the previous temperature,</li>
+     *  <li>the current simplex.</li>
+     * </ol>
+     */
+    public interface CoolingSchedule extends BiFunction<Double, Simplex, Double> {
+        /**
+         * Power-law cooling scheme:
+         * \[
+         *   T_i = T_0 * f^i
+         * \], where \( i \) is the current iteration.
+         * <p>
+         * Note: Simplex argument (of the returned function) is not used.
+         *
+         * @param f Factor by which the temperature is decreased.
+         * @return the cooling schedule.
+         */
+        static CoolingSchedule decreasingExponential(final double f) {
+            if (f <= 0 ||
+                f >= 1) {
+                throw new IllegalArgumentException("Factor out of range: " + f);
+            }
+
+            return (previousTemperature, simplex) -> f * previousTemperature;
+        }
+
+        /**
+         * Aarst and van Laarhoven (1985) scheme:
+         * \[
+         *   T_{i + 1} = \frac{T_{i}}{1 + \frac{T_i \ln(1 + \delta)}{3 \sigma}}
+         * \]
+         * <p>
+         * The simplex argument is used to compute the standard deviation
+         * (\(\sigma\)) of all the vertices' objective function value.
+         *
+         * @param delta Trajectory parameter. Values smaller than 1 entail slow
+         * convergence; values larger than 1 entail convergence to local optimum.
+         * @return the cooling schedule.
+         */
+        static CoolingSchedule aarstAndVanLaarhoven(final double delta) {
+            if (delta <= 0) {
+                throw new IllegalArgumentException("Trajectory parameter out of range: " +
+                                                   delta);
+            }
+
+            return (previousTemperature, simplex) -> {
+                // Standard deviation of the values of the objective function.
+                final StandardDeviation stddev = new StandardDeviation();
+                for (int i = 0; i < simplex.getSize(); i++) {
+                    stddev.increment(simplex.get(i).getValue());
+                }
+                final double sigma = stddev.getResult();
+
+                final double a = previousTemperature * Math.log(1 + delta);
+                final double b = 3 * sigma;
+                return previousTemperature / (1 + a / b);
+            };
+        }
+    }
+
+    /**
+     * Factory for the Metropolis check for accepting a worse state.
+     * \( e^{-|\Delta E| / T} \geq \mathrm{rand}(0, 1) \).
+     *
+     * <p>
+     * It is assumed that this check is performed <em>after</em> ensuring
+     * that the alternate state is <em>worse</em> than the current best state.
+     *
+     * @param temperature Current temperature.
+     * @return the acceptance test.
+     */
+    public DoublePredicate metropolis(final double temperature) {
+        // Absolute value takes care of both minimization and maximization cases.
+        return deltaE -> Math.exp(Math.abs(deltaE) / temperature) >= rng.nextDouble();
+    }
+}
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/HedarFukushimaTransform.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/HedarFukushimaTransform.java
new file mode 100644
index 0000000..2328b72
--- /dev/null
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/HedarFukushimaTransform.java
@@ -0,0 +1,210 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv;
+
+import java.util.Comparator;
+import java.util.List;
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.function.UnaryOperator;
+import java.util.function.DoublePredicate;
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.simple.RandomSource;
+import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
+import org.apache.commons.rng.sampling.distribution.ContinuousUniformSampler;
+import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
+import org.apache.commons.math4.legacy.optim.PointValuePair;
+import org.apache.commons.math4.legacy.optim.OptimizationData;
+
+/**
+ * DSSA algorithm.
+ *
+ * Described in
+ * <blockquote>
+ *  <em>Abdel-Rahman Hedar and Masao Fukushima (2002)</em>,
+ *  <b>
+ *   Hybrid simulated annealing and direct search method
+ *   for nonlinear unconstrained global optimization
+ *  </b>,
+ *  Optimization Methods and Software, 17:5, 891-912,
+ *  DOI: 10.1080/1055678021000030084
+ * </blockquote>
+ *
+ * <p>
+ * A note about the {@link #HedarFukushimaTransform(double) "shrink" factor}:
+ * Per DSSA's description, the simplex must keep its size during the simulated
+ * annealing (SA) phase to avoid premature convergence.  This assumes that the
+ * best candidates from the SA phase will each subsequently serve as starting
+ * point for another optimization to hone in on the local optimum.
+ * Values lower than 1 and no subsequent "best list" search correspond to the
+ * "SSA" algorithm in the above paper.
+ */
+public class HedarFukushimaTransform
+    implements Simplex.TransformFactory,
+               OptimizationData {
+    /** Shrinkage coefficient. */
+    private final double sigma;
+    /** Sampler for reflection coefficient. */
+    private final ContinuousSampler alphaSampler;
+    /** No shrink indicator. */
+    private final boolean noShrink;
+
+    /**
+     * @param sigma Shrink factor.
+     * @param rng Random generator.
+     * @throws IllegalArgumentException if {@code sigma <= 0} or
+     * {@code sigma > 1}.
+     */
+    public HedarFukushimaTransform(double sigma,
+                                   UniformRandomProvider rng) {
+        if (sigma <= 0 ||
+            sigma > 1) {
+            throw new IllegalArgumentException("Shrink factor out of range: " +
+                                               sigma);
+        }
+
+        this.sigma = sigma;
+        alphaSampler = ContinuousUniformSampler.of(rng, 0.9, 1.1);
+        noShrink = sigma == 1d;
+    }
+
+    /**
+     * @param sigma Shrink factor.
+     * @throws IllegalArgumentException if {@code sigma <= 0} or
+     * {@code sigma > 1}.
+     */
+    public HedarFukushimaTransform(double sigma) {
+        this(sigma, RandomSource.KISS.create());
+    }
+
+    /**
+     * Disable shrinking of the simplex (as mandated by DSSA).
+     */
+    public HedarFukushimaTransform() {
+        this(1d);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public UnaryOperator<Simplex> create(final MultivariateFunction evaluationFunction,
+                                         final Comparator<PointValuePair> comparator,
+                                         final DoublePredicate saAcceptance) {
+        if (saAcceptance == null) {
+            throw new IllegalArgumentException("Missing SA acceptance test");
+        }
+
+        return original -> transform(original,
+                                     saAcceptance,
+                                     evaluationFunction,
+                                     comparator);
+    }
+
+    /**
+     * Simulated annealing step (at fixed temperature).
+     *
+     * @param original Current simplex.  Points must be sorted from best to worst.
+     * @param sa Simulated annealing acceptance test.
+     * @param eval Evaluation function.
+     * @param comp Objective function comparator.
+     * @return a new instance.
+     */
+    private Simplex transform(Simplex original,
+                              DoublePredicate sa,
+                              MultivariateFunction eval,
+                              Comparator<PointValuePair> comp) {
+        final int size = original.getSize();
+        // Current best point.
+        final PointValuePair best = original.get(0);
+        final double bestValue = best.getValue();
+
+        for (int k = 1; k < size; k++) {
+            // Perform reflections of the "k" worst points.
+            final List<PointValuePair> reflected = reflectPoints(original, k, eval);
+            Collections.sort(reflected, comp);
+
+            // Check whether the best of the reflected points is better than the
+            // current overall best.
+            final PointValuePair candidate = reflected.get(0);
+            final boolean candidateIsBetter = comp.compare(candidate, best) < 0;
+            final boolean candidateIsAccepted = candidateIsBetter ||
+                sa.test(candidate.getValue() - bestValue);
+
+            if (candidateIsAccepted) {
+                // Replace worst points with the reflected points.
+                return original.replaceLast(reflected);
+            }
+        }
+
+        // No direction provided a better point.
+        return noShrink ?
+            original :
+            original.shrink(sigma, eval);
+    }
+
+    /**
+     * @param simplex Current simplex (whose points must be sorted, from best
+     * to worst).
+     * @param nPoints Number of points to reflect.
+     * The {@code nPoints} worst points will be reflected through the centroid
+     * of the {@code n + 1 - nPoints} best points.
+     * @param eval Evaluation function.
+     * @return the (unsorted) list of reflected points.
+     * @throws IllegalArgumentException if {@code nPoints < 1} or
+     * {@code nPoints > n}.
+     */
+    private List<PointValuePair> reflectPoints(Simplex simplex,
+                                               int nPoints,
+                                               MultivariateFunction eval) {
+        final int size = simplex.getSize();
+        if (nPoints < 1 ||
+            nPoints >= size) {
+            throw new IllegalArgumentException("Out of range: " + nPoints);
+        }
+
+        final int nCentroid = size - nPoints;
+        final List<PointValuePair> centroidList = simplex.asList(0, nCentroid);
+        final List<PointValuePair> reflectList = simplex.asList(nCentroid, size);
+
+        final double[] centroid = Simplex.centroid(centroidList);
+
+        final List<PointValuePair> reflected = new ArrayList<>(nPoints);
+        for (int i = 0; i < reflectList.size(); i++) {
+            reflected.add(newReflectedPoint(reflectList.get(i),
+                                            centroid,
+                                            eval));
+        }
+
+        return reflected;
+    }
+
+    /**
+     * @param point Current point.
+     * @param centroid Coordinates through which reflection must be performed.
+     * @param eval Evaluation function.
+     * @return a new point with Cartesian coordinates set to the reflection
+     * of {@code point} through {@code centroid}.
+     */
+    private PointValuePair newReflectedPoint(PointValuePair point,
+                                             double[] centroid,
+                                             MultivariateFunction eval) {
+        final double alpha = alphaSampler.sample();
+        return Simplex.newPoint(centroid,
+                                -alpha,
+                                point.getPoint(),
+                                eval);
+    }
+}
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/MultiDirectionalTransform.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/MultiDirectionalTransform.java
index e7ae9e5..6b4e1f8 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/MultiDirectionalTransform.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/MultiDirectionalTransform.java
@@ -16,8 +16,11 @@
  */
 package org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv;
 
+import java.util.List;
+import java.util.ArrayList;
 import java.util.Comparator;
 import java.util.function.UnaryOperator;
+import java.util.function.DoublePredicate;
 
 import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
 import org.apache.commons.math4.legacy.optim.PointValuePair;
@@ -29,23 +32,33 @@ import org.apache.commons.math4.legacy.optim.OptimizationData;
 public class MultiDirectionalTransform
     implements Simplex.TransformFactory,
                OptimizationData {
+    /** Reflection coefficient. */
+    private static final double ALPHA = 1;
     /** Default value for {@link #gamma}: {@value}. */
     private static final double DEFAULT_GAMMA = 2;
-    /** Default value for {@link #rho}: {@value}. */
-    private static final double DEFAULT_RHO = 0.5;
+    /** Default value for {@link #sigma}: {@value}. */
+    private static final double DEFAULT_SIGMA = 0.5;
     /** Expansion coefficient. */
     private final double gamma;
     /** Contraction coefficient. */
-    private final double rho;
+    private final double sigma;
 
     /**
      * @param gamma Expansion coefficient.
-     * @param rho Contraction coefficient.
+     * @param sigma Shrinkage coefficient.
      */
     public MultiDirectionalTransform(double gamma,
-                                     double rho) {
+                                     double sigma) {
+        if (gamma < 1) {
+            throw new IllegalArgumentException("gamma: " + gamma);
+        }
+        if (sigma < 0 ||
+            sigma > 1) {
+            throw new IllegalArgumentException("sigma: " + sigma);
+        }
+
         this.gamma = gamma;
-        this.rho = rho;
+        this.sigma = sigma;
     }
 
     /**
@@ -53,31 +66,38 @@ public class MultiDirectionalTransform
      */
     public MultiDirectionalTransform() {
         this(DEFAULT_GAMMA,
-             DEFAULT_RHO);
+             DEFAULT_SIGMA);
     }
 
     /** {@inheritDoc} */
     @Override
-    public UnaryOperator<Simplex> apply(final MultivariateFunction evaluationFunction,
-                                        final Comparator<PointValuePair> comparator) {
+    public UnaryOperator<Simplex> create(final MultivariateFunction evaluationFunction,
+                                         final Comparator<PointValuePair> comparator,
+                                         final DoublePredicate unused) {
         return original -> {
             final PointValuePair best = original.get(0);
 
             // Perform a reflection step.
-            final Simplex reflectedSimplex = transform(original, 1);
+            final Simplex reflectedSimplex = transform(original,
+                                                       ALPHA,
+                                                       comparator,
+                                                       evaluationFunction);
             final PointValuePair reflectedBest = reflectedSimplex.get(0);
 
             if (comparator.compare(reflectedBest, best) < 0) {
                 // Compute the expanded simplex.
-                final Simplex expandedSimplex = transform(original, gamma);
+                final Simplex expandedSimplex = transform(original,
+                                                          gamma,
+                                                          comparator,
+                                                          evaluationFunction);
                 final PointValuePair expandedBest = expandedSimplex.get(0);
 
-                return comparator.compare(reflectedBest, expandedBest) <= 0 ?
+                return comparator.compare(reflectedBest, expandedBest) < 0 ?
                     reflectedSimplex :
                     expandedSimplex;
             } else {
                 // Compute the contracted simplex.
-                return transform(original, rho);
+                return original.shrink(sigma, evaluationFunction);
             }
         };
     }
@@ -87,30 +107,28 @@ public class MultiDirectionalTransform
      *
      * @param original Original simplex.
      * @param coeff Linear coefficient.
+     * @param comp Fitness comparator.
+     * @param evalFunc Objective function.
      * @return the transformed simplex.
      * @throws org.apache.commons.math4.legacy.exception.TooManyEvaluationsException
      * if the maximal number of evaluations is exceeded.
      */
     private Simplex transform(Simplex original,
-                              double coeff) {
+                              double coeff,
+                              Comparator<PointValuePair> comp,
+                              MultivariateFunction evalFunc) {
         // Transformed simplex is the result a linear transformation on all
         // points except the first one.
-        final int dim = original.getDimension();
-        final int len = original.getSize();
-        final double[][] simplex = new double[len][];
-
-        simplex[0] = original.get(0).getPoint();
-        final double[] xSmallest = simplex[0];
-
-        for (int i = 1; i < len; i++) {
-            final double[] xOriginal = original.get(i).getPoint();
-            final double[] xTransformed = new double[dim];
-            for (int j = 0; j < dim; j++) {
-                xTransformed[j] = xSmallest[j] + coeff * (xSmallest[j] - xOriginal[j]);
-            }
-            simplex[i] = xTransformed;
+        final int replSize = original.getSize() - 1;
+        final List<PointValuePair> replacement = new ArrayList<>();
+        final double[] bestPoint = original.get(0).getPoint();
+        for (int i = 0; i < replSize; i++) {
+            replacement.add(Simplex.newPoint(bestPoint,
+                                             -coeff,
+                                             original.get(i + 1).getPoint(),
+                                             evalFunc));
         }
 
-        return Simplex.of(simplex);
+        return original.replaceLast(replacement).evaluate(evalFunc, comp);
     }
 }
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/NelderMeadTransform.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/NelderMeadTransform.java
index 6933c6f..775d5a3 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/NelderMeadTransform.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/NelderMeadTransform.java
@@ -18,6 +18,7 @@ package org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv;
 
 import java.util.Comparator;
 import java.util.function.UnaryOperator;
+import java.util.function.DoublePredicate;
 
 import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
 import org.apache.commons.math4.legacy.optim.PointValuePair;
@@ -74,8 +75,9 @@ public class NelderMeadTransform
 
     /** {@inheritDoc} */
     @Override
-    public UnaryOperator<Simplex> apply(final MultivariateFunction evaluationFunction,
-                                        final Comparator<PointValuePair> comparator) {
+    public UnaryOperator<Simplex> create(final MultivariateFunction evaluationFunction,
+                                         final Comparator<PointValuePair> comparator,
+                                         final DoublePredicate unused) {
         return original -> {
             Simplex newSimplex = original;
 
@@ -89,17 +91,7 @@ public class NelderMeadTransform
             final double[] xWorst = worst.getPoint();
 
             // Centroid of the best vertices, dismissing the worst point (at index n).
-            final double[] centroid = new double[n];
-            for (int i = 0; i < n; i++) {
-                final double[] x = newSimplex.get(i).getPoint();
-                for (int j = 0; j < n; j++) {
-                    centroid[j] += x[j];
-                }
-            }
-            final double scaling = 1d / n;
-            for (int j = 0; j < n; j++) {
-                centroid[j] *= scaling;
-            }
+            final double[] centroid = Simplex.centroid(newSimplex.asList().subList(0, n));
 
             // Reflection.
             final PointValuePair reflected = Simplex.newPoint(centroid,
@@ -108,7 +100,7 @@ public class NelderMeadTransform
                                                               evaluationFunction);
             if (comparator.compare(reflected, secondWorst) < 0 &&
                 comparator.compare(best, reflected) <= 0) {
-                return newSimplex.withReplacement(n, reflected);
+                return newSimplex.replaceLast(reflected);
             }
 
             if (comparator.compare(reflected, best) < 0) {
@@ -118,9 +110,9 @@ public class NelderMeadTransform
                                                                  xWorst,
                                                                  evaluationFunction);
                 if (comparator.compare(expanded, reflected) < 0) {
-                    return newSimplex.withReplacement(n, expanded);
+                    return newSimplex.replaceLast(expanded);
                 } else {
-                    return newSimplex.withReplacement(n, reflected);
+                    return newSimplex.replaceLast(reflected);
                 }
             }
 
@@ -131,7 +123,7 @@ public class NelderMeadTransform
                                                                    reflected.getPoint(),
                                                                    evaluationFunction);
                 if (comparator.compare(contracted, reflected) < 0) {
-                    return newSimplex.withReplacement(n, contracted); // Accept contracted point.
+                    return newSimplex.replaceLast(contracted); // Accept contracted point.
                 }
             } else {
                 // Inside contraction.
@@ -140,7 +132,7 @@ public class NelderMeadTransform
                                                                    xWorst,
                                                                    evaluationFunction);
                 if (comparator.compare(contracted, worst) < 0) {
-                    return newSimplex.withReplacement(n, contracted); // Accept contracted point.
+                    return newSimplex.replaceLast(contracted); // Accept contracted point.
                 }
             }
 
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/Simplex.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/Simplex.java
index a2d4d4b..bcc94aa 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/Simplex.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/Simplex.java
@@ -22,8 +22,8 @@ import java.util.ArrayList;
 import java.util.Arrays;
 import java.util.Comparator;
 import java.util.Collections;
-import java.util.function.BiFunction;
 import java.util.function.UnaryOperator;
+import java.util.function.DoublePredicate;
 
 import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
@@ -35,11 +35,9 @@ import org.apache.commons.math4.legacy.optim.OptimizationData;
 import org.apache.commons.math4.legacy.optim.PointValuePair;
 
 /**
- * Represents a simplex.
+ * Represents a <a href="https://en.wikipedia.org/wiki/Simplex">simplex</a>.
  *
  * @see SimplexOptimizer
- * @see MultiDirectionalTransform
- * @see NelderMeadTransform
  */
 public final class Simplex implements OptimizationData {
     /** Coordinates. */
@@ -105,7 +103,7 @@ public final class Simplex implements OptimizationData {
     /**
      * Builds from a given set of coordinates.
      *
-     * @param simplex Reference simplex.
+     * @param simplex Simplex coordinates.
      * @return a new instance.
      * @throws NotStrictlyPositiveException if the reference simplex does not
      * contain at least one point.
@@ -118,27 +116,17 @@ public final class Simplex implements OptimizationData {
     }
 
     /**
-     * Builds a unit hypercube simplex.
+     * Builds simplex with the given side length.
      *
-     * @param n Dimension of the simplex.
-     * @return a new instance.
-     */
-    public static Simplex of(int n) {
-        return of(n, 1d);
-    }
-
-    /**
-     * Builds a hypercube simplex with the given side length.
-     *
-     * @param n Dimension of the simplex.
+     * @param dim Space dimensions.
      * @param sideLength Length of the sides of the hypercube.
      * @return a new instance.
      */
-    public static Simplex of(int n,
-                             double sideLength) {
-        final double[] steps = new double[n];
+    public static Simplex equalSidesAlongAxes(int dim,
+                                              double sideLength) {
+        final double[] steps = new double[dim];
         Arrays.fill(steps, sideLength);
-        return of(steps);
+        return alongAxes(steps);
     }
 
     /**
@@ -147,19 +135,28 @@ public final class Simplex implements OptimizationData {
      * of a box parallel to the canonical axes. It is built as the path followed
      * while traveling from one vertex of the box to the diagonally opposite
      * vertex moving only along the box edges. The first vertex of the box will
-     * be located at the start point of the optimization.
-     * As an example, in dimension 3 a simplex has 4 vertices. Setting the
+     * be located at the origin of the coordinate system.
+     *
+     * To be used for simplex-based optimization, the simplex must be
+     * {@link #translate(double[]) translated} so that its first vertex will be
+     * the {@link org.apache.commons.math4.legacy.optim.InitialGuess initial guess}.
+     *
+     * For example, in dimension 3 a simplex has 4 vertices. Setting the
      * steps to (1, 10, 2) and the start point to (1, 1, 1) would imply the
-     * start simplex would be: { (1, 1, 1), (2, 1, 1), (2, 11, 1), (2, 11, 3) }.
-     * The first vertex would be set to the start point at (1, 1, 1) and the
-     * last vertex would be set to the diagonally opposite vertex at (2, 11, 3).
+     * initial simplex would be:
+     * <ol>
+     *  <li>(1, 1, 1),</li>
+     *  <li>(2, 1, 1),</li>
+     *  <li>(2, 11, 1),</li>
+     *  <li>(2, 11, 3).</li>
+     * </ol>
      *
      * @param steps Steps along the canonical axes representing box edges.
      * They may be negative but not zero.
      * @throws ZeroException if one of the steps is zero.
      * @return a new instance.
      */
-    public static Simplex of(double[] steps) {
+    public static Simplex alongAxes(double[] steps) {
         if (steps.length == 0) {
             throw new ZeroException();
         }
@@ -237,70 +234,64 @@ public final class Simplex implements OptimizationData {
     }
 
     /**
-     * Translates the simplex such that its first point is at the given {@code point}.
+     * Creates a (deep) copy of the simplex points.
      *
-     * @param point Coordinates of the new simplex's first point.
-     * @return a new instance.
-     * @throws DimensionMismatchException if the start point does not match
-     * simplex dimension.
+     * @return the points.
      */
-    /* package private */ Simplex translate(final double[] point) {
-        final int dim = getDimension();
-        if (dim != point.length) {
-            throw new DimensionMismatchException(dim, point.length);
-        }
-
-        final int len = points.size();
-        final double[][] coordinates = new double[len][dim];
-        final double[] current0 = points.get(0).getPoint(); // Current first point.
-
-        // Set new vertices.
-        for (int i = 0; i < len; i++) {
-            final double[] currentI = points.get(i).getPoint();
-            final double[] newI = coordinates[i];
-            for (int k = 0; k < dim; k++) {
-                newI[k] = point[k] + currentI[k] - current0[k];
-            }
-        }
+    public List<PointValuePair> asList() {
+        return asList(0, points.size());
+    }
 
-        return of(coordinates);
+    /**
+     * Generator of simplex transform.
+     *
+     * @see MultiDirectionalTransform
+     * @see NelderMeadTransform
+     * @see HedarFukushimaTransform
+     */
+    public interface TransformFactory {
+        /**
+         * Creates a simplex transformation.
+         *
+         * @param evaluationFunction Evaluation function.
+         * @param comparator Vertex fitness comparator.
+         * @param saAcceptance Simulated annealing acceptance test.
+         * @return the simplex transform operator.
+         */
+        UnaryOperator<Simplex> create(MultivariateFunction evaluationFunction,
+                                      Comparator<PointValuePair> comparator,
+                                      DoublePredicate saAcceptance);
     }
 
     /**
-     * Builds a new simplex where the given {@code point} replaces the
-     * one at {@code index} in this instance.
+     * Creates a (deep) copy of the simplex points within slots
+     * {@code from} (included) and {@code to} (excluded).
      *
-     * @param index Index of the point to replace.
-     * @param point Replacement for the point currently at {@code index}.
-     * @return a new instance.
-     * @throws IndexOutOfBoundsException if {@code index} is out of bounds.
+     * @param from Index of the first point to retrieve.
+     * @param to One past the index of the last point to retrieve.
+     * @return the points.
+     * @throws IllegalArgumentException if {@code from} and {@code to} are
+     * not within the {@code [0, n + 1]} interval (where {@code n} is the
+     * space dimension) or {@code from > to}.
      */
-    /* package private */ Simplex withReplacement(int index,
-                                                  PointValuePair point) {
-        final int len = points.size();
-        if (index < 0 ||
-            index >= len) {
-            throw new IndexOutOfBoundsException("index: " + index);
+    /* package private */ List<PointValuePair> asList(int from,
+                                                      int to) {
+        if (from < 0 ||
+            to > points.size() ||
+            from > to) {
+            throw new IllegalArgumentException("Index");
         }
 
-        final List<PointValuePair> newPoints = new ArrayList<>(len);
-        for (int i = 0; i < len; i++) {
-            final PointValuePair pv = i == index ?
-                point :
-                points.get(i);
-            newPoints.add(new PointValuePair(pv.getPoint(), pv.getValue(), false));
+        final int len = to - from;
+        final List<PointValuePair> copy = new ArrayList<>(len);
+        for (int i = from; i < to; i++) {
+            copy.add(get(i));
         }
 
-        return new Simplex(newPoints);
+        return copy;
     }
 
     /**
-     * Generator of simplex transform.
-     */
-    public interface TransformFactory
-        extends BiFunction<MultivariateFunction, Comparator<PointValuePair>, UnaryOperator<Simplex>> {}
-
-    /**
      * Utility for evaluating a point with coordinates \( a_i + s (b_i - a_i) \).
      *
      * @param a Cartesian coordinates.
@@ -324,8 +315,8 @@ public final class Simplex implements OptimizationData {
     }
 
     /**
-     * Utility for the "shrinking" a simplex: All the points will be transformed
-     * except the one at index 0.
+     * Utility for the "shrinking" a simplex: All the points will be
+     * transformed except the one at index 0.
      *
      * @param sigma Shrink factor.
      * @param function Evaluation function.
@@ -333,17 +324,111 @@ public final class Simplex implements OptimizationData {
      */
     /* package private */ Simplex shrink(double sigma,
                                          MultivariateFunction function) {
-        final int size = getSize();
-        final double[] xBest = get(0).getPoint();
-        Simplex newSimplex = this;
-        for (int i = 1; i < size; i++) {
-            final PointValuePair p = newPoint(xBest,
-                                              sigma,
-                                              get(i).getPoint(),
-                                              function);
-            newSimplex = newSimplex.withReplacement(i, p);
+        final int replSize = getSize() - 1;
+        final List<PointValuePair> replacement = new ArrayList<>();
+        final double[] bestPoint = get(0).getPoint();
+        for (int i = 0; i < replSize; i++) {
+            replacement.add(Simplex.newPoint(bestPoint,
+                                             sigma,
+                                             get(i + 1).getPoint(),
+                                             function));
+        }
+
+        return replaceLast(replacement);
+    }
+
+    /**
+     * Translates the simplex such that the first point's new coordinates
+     * will be at the given {@code point}.
+     *
+     * @param point Coordinates of the new simplex's first point.
+     * @return the translated points.
+     * @throws DimensionMismatchException if the dimensions do not match.
+     */
+    /* package private */ Simplex translate(double[] point) {
+        final int dim = point.length;
+        if (getDimension() != dim) {
+            throw new DimensionMismatchException(getDimension(), dim);
+        }
+        final int len = points.size();
+        final double[][] coordinates = new double[len][dim];
+        final double[] current0 = points.get(0).getPoint(); // Current first point.
+
+        // Set new vertices.
+        for (int i = 0; i < len; i++) {
+            final double[] currentI = points.get(i).getPoint();
+
+            final double[] newI = coordinates[i];
+            for (int k = 0; k < dim; k++) {
+                newI[k] = point[k] + currentI[k] - current0[k];
+            }
+        }
+
+        return new Simplex(coordinates);
+    }
+
+    /**
+     * Creates a new simplex where the given {@code point} replaces the one at the
+     * last position.
+     * Caveat: No check is done that the resulting set of points forms is a simplex.
+     *
+     * @param point Point.
+     * @return a new instance.
+     */
+    /* package private */ Simplex replaceLast(PointValuePair point) {
+        final List<PointValuePair> newPoints = asList(0, getDimension()); // Deep copy.
+        newPoints.add(new PointValuePair(point.getPoint(), // Deep copy.
+                                         point.getValue(),
+                                         false));
+
+        return new Simplex(newPoints);
+    }
+
+    /**
+     * Replace the last points of the simplex with the points from the given
+     * {@code replacement} list.
+     * Caveat: No check is done that the resulting set of points is a simplex.
+     *
+     * @param replacement List of points that will replace the last points of
+     * the {@code simplex}.
+     * @return a new instance.
+     */
+    /* package private */ Simplex replaceLast(List<PointValuePair> replacement) {
+        final int nPoints = replacement.size();
+        final int from = points.size() - nPoints;
+        final List<PointValuePair> newPoints = asList(0, from); // Deep copy.
+
+
+        for (int i = 0; i < nPoints; i++) {
+            final PointValuePair p = replacement.get(i);
+            newPoints.add(new PointValuePair(p.getPoint(), // Deep copy.
+                                             p.getValue(),
+                                             false));
+        }
+
+        return new Simplex(newPoints);
+    }
+
+    /**
+     * @param list List of simplex points.
+     * @return the centroid of the points in the given {@code list}.
+     */
+    /* package private */ static double[] centroid(List<PointValuePair> list) {
+        final double[] centroid = list.get(0).getPoint();
+
+        final int nPoints = list.size();
+        final int dim = centroid.length;
+        for (int i = 1; i < nPoints; i++) {
+            final double[] p = list.get(i).getPoint();
+            for (int k = 0; k < dim; k++) {
+                centroid[k] += p[k];
+            }
+        }
+
+        for (int k = 0; k < dim; k++) {
+            centroid[k] /= nPoints;
         }
 
-        return newSimplex;
+        return centroid;
     }
 }
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizer.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizer.java
index 60af50d..ac367af 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizer.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizer.java
@@ -16,18 +16,29 @@
  */
 package org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv;
 
+import java.util.Arrays;
+import java.util.List;
+import java.util.ArrayList;
 import java.util.Comparator;
+import java.util.Collections;
 import java.util.function.UnaryOperator;
 
+import org.apache.commons.math4.legacy.core.MathArrays;
 import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
 import org.apache.commons.math4.legacy.exception.MathUnsupportedOperationException;
+import org.apache.commons.math4.legacy.exception.ConvergenceException;
 import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
 import org.apache.commons.math4.legacy.optim.ConvergenceChecker;
 import org.apache.commons.math4.legacy.optim.OptimizationData;
 import org.apache.commons.math4.legacy.optim.PointValuePair;
 import org.apache.commons.math4.legacy.optim.SimpleValueChecker;
+import org.apache.commons.math4.legacy.optim.InitialGuess;
+import org.apache.commons.math4.legacy.optim.MaxEval;
 import org.apache.commons.math4.legacy.optim.nonlinear.scalar.GoalType;
 import org.apache.commons.math4.legacy.optim.nonlinear.scalar.MultivariateOptimizer;
+import org.apache.commons.math4.legacy.optim.nonlinear.scalar.SimulatedAnnealing;
+import org.apache.commons.math4.legacy.optim.nonlinear.scalar.PopulationSize;
+import org.apache.commons.math4.legacy.optim.nonlinear.scalar.ObjectiveFunction;
 
 /**
  * This class implements simplex-based direct search optimization.
@@ -57,8 +68,9 @@ import org.apache.commons.math4.legacy.optim.nonlinear.scalar.MultivariateOptimi
  * an instance of this class will register the following data:
  * <ul>
  *  <li>{@link Simplex}</li>
- *  <li>{@link Simplex.TransformFactory} (either {@link NelderMeadTransform}
- *   or {@link MultiDirectionalTransform})</li>
+ *  <li>{@link Simplex.TransformFactory}</li>
+ *  <li>{@link SimulatedAnnealing}</li>
+ *  <li>{@link PopulationSize}</li>
  * </ul>
  *
  * <p>
@@ -72,18 +84,36 @@ import org.apache.commons.math4.legacy.optim.nonlinear.scalar.MultivariateOptimi
  * <p>
  * Convergence is considered achieved when <em>all</em> the simplex points
  * have converged.
+ * <p>
+ * Whenever {@link SimulatedAnnealing simulated annealing (SA)} is activated,
+ * the SA phase will terminate when the temperature reaches \( 10^{-2} \).
+ * Indeed the chance of accepting the alternative state is vanishingly small.
+ * It does not mean however that the optimization has converged.
+ * In such a case, the {@link PopulationSize} argument to method
+ * {@link #optimize(OptimizationData[]) optimize} will trigger an additional
+ * "best list" search.
  *
  * <p>
  * This implementation does not directly support constrained optimization
  * with simple bounds.
  * The call to {@link #optimize(OptimizationData[]) optimize} will throw
  * {@link MathUnsupportedOperationException} if bounds are passed to it.
+ *
+ * @see NelderMeadTransform
+ * @see MultiDirectionalTransform
+ * @see HedarFukushimaTransform
  */
 public class SimplexOptimizer extends MultivariateOptimizer {
+    /** Final temperature for simulated annealing. */
+    private static final double SA_FINAL_TEMPERATURE = 1e-2;
     /** Simplex update function factory. */
     private Simplex.TransformFactory updateRule;
-    /** Current simplex. */
-    private Simplex simplex;
+    /** Initial simplex. */
+    private Simplex initialSimplex;
+    /** Simulated annealing setup (optional). */
+    private SimulatedAnnealing simulatedAnnealing = null;
+    /** Number of additional optimizations (optional). */
+    private int bestListSize = 0;
 
     /**
      * @param checker Convergence checker.
@@ -117,36 +147,79 @@ public class SimplexOptimizer extends MultivariateOptimizer {
             return isMinim ? Double.compare(v1, v2) : Double.compare(v2, v1);
         };
 
-        final UnaryOperator<Simplex> update = updateRule.apply(evalFunc, comparator);
+        // Start points for additional search.
+        final List<PointValuePair> bestList = new ArrayList<>();
 
-        // Initialize search.
-        simplex = simplex.translate(getStartPoint()).evaluate(evalFunc, comparator);
+        Simplex currentSimplex = initialSimplex.translate(getStartPoint()).evaluate(evalFunc, comparator);
+        double temperature = Double.NaN; // Only used with simulated annealing.
+        Simplex previousSimplex = null;
+
+        if (simulatedAnnealing != null) {
+            temperature =
+                temperature(currentSimplex.get(0),
+                            currentSimplex.get(currentSimplex.getDimension()),
+                            simulatedAnnealing.getStartProbability());
+        }
 
-        Simplex previous = null;
-        final ConvergenceChecker<PointValuePair> checker = getConvergenceChecker();
         while (true) {
-            if (previous != null) { // Skip check at first iteration.
-                boolean converged = true;
-                for (int i = 0; i < simplex.getSize(); i++) {
-                    if (!checker.converged(getIterations(),
-                                           previous.get(i),
-                                           simplex.get(i))) {
-                        converged = false;
-                        break;
-                    }
-                }
-                if (converged) {
-                    // We have found an optimum.
-                    return simplex.get(0);
+            if (previousSimplex != null) { // Skip check at first iteration.
+                if (hasConverged(previousSimplex, currentSimplex)) {
+                    return currentSimplex.get(0);
                 }
             }
 
             // We still need to search.
-            previous = simplex;
-            simplex = update.apply(simplex).evaluate(evalFunc, comparator);
+            previousSimplex = currentSimplex;
+
+            if (simulatedAnnealing != null) {
+                // Update current temperature.
+                temperature =
+                    simulatedAnnealing.getCoolingSchedule().apply(temperature,
+                                                                  currentSimplex);
+
+                final double endTemperature =
+                    temperature(currentSimplex.get(0),
+                                currentSimplex.get(currentSimplex.getDimension()),
+                                simulatedAnnealing.getEndProbability());
+
+                if (temperature < endTemperature) {
+                    break;
+                }
+
+                final UnaryOperator<Simplex> update =
+                    updateRule.create(evalFunc,
+                                      comparator,
+                                      simulatedAnnealing.metropolis(temperature));
+
+                for (int i = 0; i < simulatedAnnealing.getEpochDuration(); i++) {
+                    currentSimplex = update.apply(currentSimplex).evaluate(evalFunc, comparator);
+                }
+            } else {
+                // No simulated annealing.
+                final UnaryOperator<Simplex> update =
+                    updateRule.create(evalFunc, comparator, null);
+
+                currentSimplex = update.apply(currentSimplex).evaluate(evalFunc, comparator);
+            }
+
+            if (bestListSize != 0) {
+                // Store best points.
+                for (int i = 0; i < currentSimplex.getSize(); i++) {
+                    keepIfBetter(currentSimplex.get(i), comparator, bestList, bestListSize);
+                }
+            }
 
             incrementIterationCount();
         }
+
+        // No convergence.
+        if (bestList.isEmpty()) {
+            // Bail out.
+            throw new ConvergenceException();
+        } else {
+            // Additional optimizations.
+            return bestListSearch(evalFunc, comparator, bestList);
+        }
     }
 
     /**
@@ -158,6 +231,8 @@ public class SimplexOptimizer extends MultivariateOptimizer {
      * <ul>
      *  <li>{@link Simplex}</li>
      *  <li>{@link Simplex.TransformFactory}</li>
+     *  <li>{@link SimulatedAnnealing}</li>
+     *  <li>{@link PopulationSize}</li>
      * </ul>
      */
     @Override
@@ -169,11 +244,39 @@ public class SimplexOptimizer extends MultivariateOptimizer {
         // if not provided in the argument list.
         for (OptimizationData data : optData) {
             if (data instanceof Simplex) {
-                simplex = (Simplex) data;
+                initialSimplex = (Simplex) data;
             } else if (data instanceof Simplex.TransformFactory) {
                 updateRule = (Simplex.TransformFactory) data;
+            } else if (data instanceof SimulatedAnnealing) {
+                simulatedAnnealing = (SimulatedAnnealing) data;
+            } else if (data instanceof PopulationSize) {
+                bestListSize = ((PopulationSize) data).getPopulationSize();
+            }
+        }
+    }
+
+    /**
+     * Detects whether the simplex has shrunk below the user-defined
+     * tolerance.
+     *
+     * @param previous Simplex at previous iteration.
+     * @param current Simplex at current iteration.
+     * @return {@code true} if convergence is considered achieved.
+     */
+    private boolean hasConverged(Simplex previous,
+                                 Simplex current) {
+        final ConvergenceChecker<PointValuePair> checker = getConvergenceChecker();
+
+        for (int i = 0; i < current.getSize(); i++) {
+            final PointValuePair prev = previous.get(i);
+            final PointValuePair curr = current.get(i);
+
+            if (!checker.converged(getIterations(), prev, curr)) {
+                return false;
             }
         }
+
+        return true;
     }
 
     /**
@@ -186,7 +289,7 @@ public class SimplexOptimizer extends MultivariateOptimizer {
         if (updateRule == null) {
             throw new NullPointerException("No update rule");
         }
-        if (simplex == null) {
+        if (initialSimplex == null) {
             throw new NullPointerException("No initial simplex");
         }
         if (getLowerBound() != null ||
@@ -194,4 +297,154 @@ public class SimplexOptimizer extends MultivariateOptimizer {
             throw new MathUnsupportedOperationException(LocalizedFormats.CONSTRAINT);
         }
     }
+
+    /**
+     * Computes the temperature as a function of the acceptance probability
+     * and the fitness difference between two of the simplex vertices (usually
+     * the best and worst points).
+     *
+     * @param p1 Simplex point.
+     * @param p2 Simplex point.
+     * @param prob Acceptance probability.
+     * @return the temperature.
+     */
+    private double temperature(PointValuePair p1,
+                               PointValuePair p2,
+                               double prob) {
+        return -Math.abs(p1.getValue() - p2.getValue()) / Math.log(prob);
+    }
+
+    /**
+     * Stores the given {@code candidate} if its fitness is better than
+     * that of the last (assumed to be the worst) point in {@code list}.
+     *
+     * @param candidate Point to be stored.
+     * @param comp Fitness comparator.
+     * @param list Starting points (modified in-place).
+     * @param max Maximum size of the {@code list}.
+     */
+    private static void keepIfBetter(PointValuePair candidate,
+                                     Comparator<PointValuePair> comp,
+                                     List<PointValuePair> list,
+                                     int max) {
+        final int listSize = list.size();
+        final double[] candidatePoint = candidate.getPoint();
+        if (listSize == 0) {
+            list.add(candidate);
+        } else if (listSize < max) {
+            // List is not fully populated yet.
+            for (PointValuePair p : list) {
+                final double[] pPoint = p.getPoint();
+                if (Arrays.equals(pPoint, candidatePoint)) {
+                    // Point was already stored.
+                    return;
+                } else {
+                    // Store candidate.
+                    list.add(candidate);
+                    return;
+                }
+            }
+        } else {
+            final int last = max - 1;
+            if (comp.compare(candidate, list.get(last)) < 0) {
+                for (PointValuePair p : list) {
+                    final double[] pPoint = p.getPoint();
+                    if (Arrays.equals(pPoint, candidatePoint)) {
+                        // Point was already stored.
+                        return;
+                    }
+                }
+
+                // Store better candidate and reorder the list.
+                list.set(last, candidate);
+                Collections.sort(list, comp);
+            }
+        }
+    }
+
+    /**
+     * Computes the smallest distance between the given {@code point}
+     * and any of the other points in the {@code list}.
+     *
+     * @param point Point.
+     * @param list List.
+     * @return the smallest distance.
+     */
+    private static double shortestDistance(PointValuePair point,
+                                           List<PointValuePair> list) {
+        double minDist = Double.POSITIVE_INFINITY;
+
+        final double[] p = point.getPoint();
+        for (PointValuePair other : list) {
+            final double[] pOther = other.getPoint();
+            if (!Arrays.equals(p, pOther)) {
+                final double dist = MathArrays.distance(p, pOther);
+                if (dist < minDist) {
+                    minDist = dist;
+                }
+            }
+        }
+
+        return minDist;
+    }
+
+    /**
+     * Perform additional optimizations.
+     *
+     * @param evalFunc Objective function.
+     * @param comp Fitness comparator.
+     * @param starts Starting points.
+     * @return the optimum.
+     */
+    private PointValuePair bestListSearch(MultivariateFunction evalFunc,
+                                          Comparator<PointValuePair> comp,
+                                          List<PointValuePair> starts) {
+        PointValuePair best = starts.get(0); // Overall best result.
+
+        // Additional local optimizations using each of the best
+        // points visited during the main search.
+        for (PointValuePair p : starts) {
+            // Find shortest distance to the other points.
+            final double dist = shortestDistance(p, starts);
+            final double[] init = p.getPoint();
+            // Create smaller initial simplex.
+            final Simplex simplex = Simplex.equalSidesAlongAxes(init.length,
+                                                                0.5 * dist);
+
+            final PointValuePair r = directSearch(init,
+                                                  simplex,
+                                                  evalFunc,
+                                                  getConvergenceChecker(),
+                                                  getGoalType());
+            if (comp.compare(r, best) < 0) {
+                best = r; // New overall best.
+            }
+        }
+
+        return best;
+    }
+
+    /**
+     * @param init Start point.
+     * @param simplex Initial simplex.
+     * @param eval Objective function.
+     * Note: It is assumed that evaluations of this function are
+     * incrementing the main counter.
+     * @param checker Convergence checker.
+     * @param goalType Whether to minimize or maximize the objective function.
+     * @return the optimum.
+     */
+    private static PointValuePair directSearch(double[] init,
+                                               Simplex simplex,
+                                               MultivariateFunction eval,
+                                               ConvergenceChecker<PointValuePair> checker,
+                                               GoalType goalType) {
+        final SimplexOptimizer optim = new SimplexOptimizer(checker);
+        return optim.optimize(MaxEval.unlimited(),
+                              new ObjectiveFunction(eval),
+                              goalType,
+                              new InitialGuess(init),
+                              simplex,
+                              new NelderMeadTransform());
+    }
 }
diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/MultivariateFunctionPenaltyAdapterTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/MultivariateFunctionPenaltyAdapterTest.java
index 01ac256..197350f 100644
--- a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/MultivariateFunctionPenaltyAdapterTest.java
+++ b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/MultivariateFunctionPenaltyAdapterTest.java
@@ -42,7 +42,7 @@ public class MultivariateFunctionPenaltyAdapterTest {
         final PointValuePair optimum
             = optimizer.optimize(new MaxEval(300),
                                  new ObjectiveFunction(wrapped),
-                                 Simplex.of(new double[] { 1.0, 0.5 }),
+                                 Simplex.alongAxes(new double[] { 1.0, 0.5 }),
                                  new NelderMeadTransform(),
                                  GoalType.MINIMIZE,
                                  new InitialGuess(new double[] { 1.5, 2.25 }));
@@ -65,7 +65,7 @@ public class MultivariateFunctionPenaltyAdapterTest {
         final PointValuePair optimum
             = optimizer.optimize(new MaxEval(300),
                                  new ObjectiveFunction(wrapped),
-                                 Simplex.of(new double[] { 1.0, 0.5 }),
+                                 Simplex.alongAxes(new double[] { 1.0, 0.5 }),
                                  new NelderMeadTransform(),
                                  GoalType.MINIMIZE,
                                  new InitialGuess(new double[] { -1.5, 4.0 }));
@@ -88,7 +88,7 @@ public class MultivariateFunctionPenaltyAdapterTest {
         final PointValuePair optimum
             = optimizer.optimize(new MaxEval(600),
                                  new ObjectiveFunction(wrapped),
-                                 Simplex.of(new double[] { 1.0, 0.5 }),
+                                 Simplex.alongAxes(new double[] { 1.0, 0.5 }),
                                  new NelderMeadTransform(),
                                  GoalType.MINIMIZE,
                                  new InitialGuess(new double[] { -1.5, 4.0 }));
@@ -113,7 +113,7 @@ public class MultivariateFunctionPenaltyAdapterTest {
         final PointValuePair optimum
             = optimizer.optimize(new MaxEval(300),
                                  new ObjectiveFunction(wrapped),
-                                 Simplex.of(new double[] { 1.0, 0.5 }),
+                                 Simplex.alongAxes(new double[] { 1.0, 0.5 }),
                                  new NelderMeadTransform(),
                                  GoalType.MINIMIZE,
                                  new InitialGuess(new double[] { -1.5, 4.0 }));
@@ -138,7 +138,7 @@ public class MultivariateFunctionPenaltyAdapterTest {
         final PointValuePair optimum
             = optimizer.optimize(new MaxEval(400),
                                  new ObjectiveFunction(wrapped),
-                                 Simplex.of(new double[] { 1.0, 0.5 }),
+                                 Simplex.alongAxes(new double[] { 1.0, 0.5 }),
                                  new NelderMeadTransform(),
                                  GoalType.MINIMIZE,
                                  new InitialGuess(new double[] { -1.5, 4.0 }));
diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/TestFunction.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/TestFunction.java
index fe457d4..2dd7048 100644
--- a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/TestFunction.java
+++ b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/TestFunction.java
@@ -180,6 +180,16 @@ public enum TestFunction {
                 }
                 return f;
             };
+        }),
+    PARABOLA(dim -> {
+            return x -> {
+                double f = 0;
+                for (int i = 0; i < dim; i++) {
+                    final double xi = x[i];
+                    f += xi * xi;
+                }
+                return f;
+            };
         });
 
     /** Template for variable dimension. */
@@ -197,13 +207,29 @@ public enum TestFunction {
      * @return the function for the given dimension.
      */
     public MultivariateFunction withDimension(final int dim) {
-        final MultivariateFunction f = generator.apply(dim);
-        return x -> {
-            if (x.length != dim) {
-                throw new IllegalArgumentException("Dimension mismatch: " + x.length +
-                                                   "(expected: " + dim + ")");
+        return new MultivariateFunction() {
+            /** Delegate. */
+            private final MultivariateFunction f = generator.apply(dim);
+
+            @Override
+            public double value(double[] x) {
+                if (x.length != dim) {
+                    throw new IllegalArgumentException("Dimension mismatch: " + x.length +
+                                                       "(expected: " + dim + ")");
+                }
+                return f.value(x);
+            };
+
+            @Override
+            public String toString() {
+                final StringBuilder sb = new StringBuilder();
+                sb.append("[")
+                    .append(TestFunction.this.toString())
+                    .append(" dim=")
+                    .append(dim)
+                    .append("]");
+                return sb.toString();
             }
-            return f.value(x);
         };
     }
 }
diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizerHedarFukushimaTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizerHedarFukushimaTest.java
new file mode 100644
index 0000000..9b266f4
--- /dev/null
+++ b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizerHedarFukushimaTest.java
@@ -0,0 +1,161 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv;
+
+import java.util.Arrays;
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.simple.RandomSource;
+import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
+import org.apache.commons.math4.legacy.optim.InitialGuess;
+import org.apache.commons.math4.legacy.optim.MaxEval;
+import org.apache.commons.math4.legacy.optim.PointValuePair;
+import org.apache.commons.math4.legacy.optim.nonlinear.scalar.GoalType;
+import org.apache.commons.math4.legacy.optim.nonlinear.scalar.ObjectiveFunction;
+import org.apache.commons.math4.legacy.optim.nonlinear.scalar.SimulatedAnnealing;
+import org.apache.commons.math4.legacy.optim.nonlinear.scalar.PopulationSize;
+import org.apache.commons.math4.legacy.optim.nonlinear.scalar.TestFunction;
+import org.apache.commons.math4.legacy.core.MathArrays;
+import org.junit.Assert;
+import org.junit.Test;
+
+/**
+ * Tests for {@link HedarFukushimaTransform}.
+ */
+public class SimplexOptimizerHedarFukushimaTest {
+    @Test
+    public void testParabola2D() {
+        final int dim = 2;
+        doTest(TestFunction.PARABOLA.withDimension(dim),
+               OptimTestUtils.point(new double[] { -4, 8 }, 5e-1),
+               GoalType.MINIMIZE,
+               250,
+               new PointValuePair(OptimTestUtils.point(dim, 0.0), 0.0),
+               1e-4);
+    }
+
+    @Test
+    public void testRosenbrock2D() {
+        final int dim = 2;
+        doTest(TestFunction.ROSENBROCK.withDimension(dim),
+               OptimTestUtils.point(new double[] { -4, 8 }, 5e-1),
+               GoalType.MINIMIZE,
+               600,
+               new PointValuePair(OptimTestUtils.point(dim, 1.0), 0.0),
+               1e-4);
+    }
+
+    @Test
+    public void testRosenbrock5D() {
+        final int dim = 5;
+        doTest(TestFunction.ROSENBROCK.withDimension(dim),
+               OptimTestUtils.point(new double[] { 2, -1, -4, 5, 9 }, 5e-1),
+               GoalType.MINIMIZE,
+               4650,
+               new PointValuePair(OptimTestUtils.point(dim, 1.0), 0.0),
+               5e-4);
+    }
+
+    @Test
+    public void testPowell() {
+        final int dim = 4;
+        doTest(TestFunction.POWELL.withDimension(dim),
+               OptimTestUtils.point(new double[] { 3, -1, 0, 1 }, 1e-1),
+               GoalType.MINIMIZE,
+               1000,
+               new PointValuePair(OptimTestUtils.point(dim, 0.0), 0.0),
+               1e-2);
+    }
+
+    /**
+     * @param func Function to optimize.
+     * @param startPoint Starting point.
+     * @param goal Minimization or maximization.
+     * @param maxEvaluations Maximum number of evaluations.
+     * @param expected Expected optimum.
+     * @param tol Tolerance for checking that the optimum is correct.
+     */
+    private void doTest(MultivariateFunction func,
+                        double[] startPoint,
+                        GoalType goal,
+                        int maxEvaluations,
+                        PointValuePair expected,
+                        double tol) {
+        doTest(func,
+               startPoint,
+               goal,
+               maxEvaluations,
+               Simplex.equalSidesAlongAxes(startPoint.length, 0.2),
+               expected,
+               tol);
+    }
+
+    /**
+     * @param func Function to optimize.
+     * @param startPoint Starting point.
+     * @param goal Minimization or maximization.
+     * @param maxEvaluations Maximum number of evaluations.
+     * @param simplexSteps Initial simplex.
+     * @param expected Expected optimum.
+     * @param tol Tolerance for checking that the optimum is correct.
+     */
+    private void doTest(MultivariateFunction func,
+                        double[] startPoint,
+                        GoalType goal,
+                        int maxEvaluations,
+                        Simplex simplex,
+                        PointValuePair expected,
+                        double tol) {
+        final String name = func.toString();
+        final int dim = startPoint.length;
+
+        // Simulated annealing setup.
+        final SimulatedAnnealing.CoolingSchedule coolSched =
+            SimulatedAnnealing.CoolingSchedule.decreasingExponential(0.7);
+        final UniformRandomProvider rng = RandomSource.KISS.create();
+        final SimulatedAnnealing sa = new SimulatedAnnealing(dim,
+                                                             0.9,
+                                                             1e-13,
+                                                             coolSched,
+                                                             rng);
+        final PopulationSize popSize = new PopulationSize(dim);
+
+        final int maxEval = Math.max(maxEvaluations, 25000);
+        final SimplexOptimizer optim = new SimplexOptimizer(-1, 1e-8);
+        final PointValuePair result = optim.optimize(new MaxEval(maxEval),
+                                                     new ObjectiveFunction(func),
+                                                     goal,
+                                                     new InitialGuess(startPoint),
+                                                     simplex,
+                                                     new HedarFukushimaTransform(),
+                                                     popSize,
+                                                     sa);
+        final double[] endPoint = result.getPoint();
+        final double funcValue = result.getValue();
+        Assert.assertEquals(name + ": value at " + Arrays.toString(endPoint),
+                            expected.getValue(),
+                            funcValue, 1e-2);
+
+        final double dist = MathArrays.distance(expected.getPoint(),
+                                                endPoint);
+        Assert.assertEquals(name + ": distance to optimum", 0d, dist, tol);
+
+        final int nEval = optim.getEvaluations();
+        Assert.assertTrue(name + ": nEval=" + nEval,
+                          nEval < maxEvaluations);
+    }
+}
diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizerMultiDirectionalTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizerMultiDirectionalTest.java
index 5bde7d8..a652cfb 100644
--- a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizerMultiDirectionalTest.java
+++ b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizerMultiDirectionalTest.java
@@ -33,6 +33,9 @@ import org.junit.Assert;
 import org.junit.Test;
 import org.junit.Ignore;
 
+/**
+ * Tests for {@link MultiDirectionalTransform}.
+ */
 public class SimplexOptimizerMultiDirectionalTest {
     private static final int DIM = 13;
 
@@ -45,7 +48,7 @@ public class SimplexOptimizerMultiDirectionalTest {
                            new ObjectiveFunction(fourExtrema),
                            GoalType.MINIMIZE,
                            new InitialGuess(new double[] { -3, 0 }),
-                           Simplex.of(new double[] { 0.2, 0.2 }),
+                           Simplex.alongAxes(new double[] { 0.2, 0.2 }),
                            new MultiDirectionalTransform(),
                            new SimpleBounds(new double[] { -5, -1 },
                                             new double[] { 5, 1 }));
@@ -59,7 +62,7 @@ public class SimplexOptimizerMultiDirectionalTest {
                                                      new ObjectiveFunction(function),
                                                      GoalType.MAXIMIZE,
                                                      new InitialGuess(function.getMaximumPosition()),
-                                                     Simplex.of(2),
+                                                     Simplex.equalSidesAlongAxes(2, 1d),
                                                      new MultiDirectionalTransform());
         final double EPSILON = 1e-5;
         final double expectedMaximum = function.getMaximum();
@@ -79,7 +82,7 @@ public class SimplexOptimizerMultiDirectionalTest {
                OptimTestUtils.point(new double[] {-3, 0}, 1e-1),
                GoalType.MINIMIZE,
                105,
-               Simplex.of(OptimTestUtils.point(2, 0.2, 1e-2)),
+               Simplex.alongAxes(OptimTestUtils.point(2, 0.2, 1e-2)),
                new PointValuePair(new double[] {f.xM, f.yP}, f.valueXmYp),
                1e-6);
     }
@@ -90,7 +93,7 @@ public class SimplexOptimizerMultiDirectionalTest {
                OptimTestUtils.point(new double[] {-3, 0}, 1e-1),
                GoalType.MAXIMIZE,
                100,
-               Simplex.of(OptimTestUtils.point(2, 0.2, 1e-2)),
+               Simplex.alongAxes(OptimTestUtils.point(2, 0.2, 1e-2)),
                new PointValuePair(new double[] {f.xM, f.yM}, f.valueXmYm),
                1e-6);
     }
@@ -101,7 +104,7 @@ public class SimplexOptimizerMultiDirectionalTest {
                OptimTestUtils.point(new double[] {1, 0}, 1e-1),
                GoalType.MINIMIZE,
                100,
-               Simplex.of(OptimTestUtils.point(2, 0.2, 1e-2)),
+               Simplex.alongAxes(OptimTestUtils.point(2, 0.2, 1e-2)),
                new PointValuePair(new double[] {f.xP, f.yM}, f.valueXpYm),
                1e-6);
     }
@@ -112,7 +115,7 @@ public class SimplexOptimizerMultiDirectionalTest {
                OptimTestUtils.point(new double[] {1, 0}, 1e-1),
                GoalType.MAXIMIZE,
                110,
-               Simplex.of(OptimTestUtils.point(2, 0.2, 1e-2)),
+               Simplex.alongAxes(OptimTestUtils.point(2, 0.2, 1e-2)),
                new PointValuePair(new double[] {f.xP, f.yP}, f.valueXpYp),
                1e-6);
     }
@@ -130,7 +133,7 @@ public class SimplexOptimizerMultiDirectionalTest {
 
     @Test
     public void testPowell() {
-        final int dim =4;
+        final int dim = 4;
         doTest(TestFunction.POWELL.withDimension(dim),
                OptimTestUtils.point(new double[] { 3, -1, 0, 1 }, 1e-1),
                GoalType.MINIMIZE,
@@ -278,7 +281,7 @@ public class SimplexOptimizerMultiDirectionalTest {
                startPoint,
                goal,
                maxEvaluations,
-               Simplex.of(startPoint.length, 1),
+               Simplex.equalSidesAlongAxes(startPoint.length, 1),
                expected,
                tol);
     }
@@ -307,7 +310,7 @@ public class SimplexOptimizerMultiDirectionalTest {
                                                      new InitialGuess(startPoint),
                                                      simplex,
                                                      new NelderMeadTransform());
-        final String name = func.getClass().getSimpleName();
+        final String name = func.toString();
 
         final double[] endPoint = result.getPoint();
         final double funcValue = result.getValue();
diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizerNelderMeadTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizerNelderMeadTest.java
index acc1ede..e2ed741 100644
--- a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizerNelderMeadTest.java
+++ b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizerNelderMeadTest.java
@@ -37,6 +37,9 @@ import org.junit.Assert;
 import org.junit.Test;
 import org.junit.Ignore;
 
+/**
+ * Tests for {@link NelderMeadTransform}.
+ */
 public class SimplexOptimizerNelderMeadTest {
     private static final int DIM = 13;
 
@@ -49,7 +52,7 @@ public class SimplexOptimizerNelderMeadTest {
                            new ObjectiveFunction(fourExtrema),
                            GoalType.MINIMIZE,
                            new InitialGuess(new double[] { -3, 0 }),
-                           Simplex.of(new double[] { 0.2, 0.2 }),
+                           Simplex.alongAxes(new double[] { 0.2, 0.2 }),
                            new NelderMeadTransform(),
                            new SimpleBounds(new double[] { -5, -1 },
                                             new double[] { 5, 1 }));
@@ -74,7 +77,7 @@ public class SimplexOptimizerNelderMeadTest {
                                new ObjectiveFunction(ls),
                                GoalType.MINIMIZE,
                                new InitialGuess(new double[] { 10, 10 }),
-                               Simplex.of(2),
+                               Simplex.equalSidesAlongAxes(2, 1d),
                                new NelderMeadTransform());
         Assert.assertEquals( 2, optimum.getPointRef()[0], 1e-3);
         Assert.assertEquals(-3, optimum.getPointRef()[1], 4e-4);
@@ -103,7 +106,7 @@ public class SimplexOptimizerNelderMeadTest {
                                new ObjectiveFunction(ls),
                                GoalType.MINIMIZE,
                                new InitialGuess(new double[] { 10, 10 }),
-                               Simplex.of(2),
+                               Simplex.equalSidesAlongAxes(2, 1d),
                                new NelderMeadTransform());
         Assert.assertEquals( 2, optimum.getPointRef()[0], 1e-4);
         Assert.assertEquals(-3, optimum.getPointRef()[1], 8e-4);
@@ -134,7 +137,7 @@ public class SimplexOptimizerNelderMeadTest {
                                  new ObjectiveFunction(ls),
                                  GoalType.MINIMIZE,
                                  new InitialGuess(new double[] { 10, 10 }),
-                                 Simplex.of(2),
+                                 Simplex.equalSidesAlongAxes(2, 1d),
                                  new NelderMeadTransform());
         Assert.assertEquals( 2, optimum.getPointRef()[0], 1e-2);
         Assert.assertEquals(-3, optimum.getPointRef()[1], 1e-2);
@@ -152,7 +155,7 @@ public class SimplexOptimizerNelderMeadTest {
                            new ObjectiveFunction(TestFunction.POWELL.withDimension(dim)),
                            GoalType.MINIMIZE,
                            new InitialGuess(new double[] { 3, -1, 0, 1 }),
-                           Simplex.of(dim),
+                           Simplex.equalSidesAlongAxes(dim, 1d),
                            new NelderMeadTransform());
     }
 
@@ -163,7 +166,7 @@ public class SimplexOptimizerNelderMeadTest {
                OptimTestUtils.point(new double[] {-3, 0}, 1e-1),
                GoalType.MINIMIZE,
                105,
-               Simplex.of(OptimTestUtils.point(2, 0.2, 1e-2)),
+               Simplex.alongAxes(OptimTestUtils.point(2, 0.2, 1e-2)),
                new PointValuePair(new double[] {f.xM, f.yP}, f.valueXmYp),
                1e-6);
     }
@@ -174,7 +177,7 @@ public class SimplexOptimizerNelderMeadTest {
                OptimTestUtils.point(new double[] {-3, 0}, 1e-1),
                GoalType.MAXIMIZE,
                100,
-               Simplex.of(OptimTestUtils.point(2, 0.2, 1e-2)),
+               Simplex.alongAxes(OptimTestUtils.point(2, 0.2, 1e-2)),
                new PointValuePair(new double[] {f.xM, f.yM}, f.valueXmYm),
                1e-6);
     }
@@ -185,7 +188,7 @@ public class SimplexOptimizerNelderMeadTest {
                OptimTestUtils.point(new double[] {1, 0}, 1e-1),
                GoalType.MINIMIZE,
                100,
-               Simplex.of(OptimTestUtils.point(2, 0.2, 1e-2)),
+               Simplex.alongAxes(OptimTestUtils.point(2, 0.2, 1e-2)),
                new PointValuePair(new double[] {f.xP, f.yM}, f.valueXpYm),
                1e-6);
     }
@@ -196,7 +199,7 @@ public class SimplexOptimizerNelderMeadTest {
                OptimTestUtils.point(new double[] {1, 0}, 1e-1),
                GoalType.MAXIMIZE,
                110,
-               Simplex.of(OptimTestUtils.point(2, 0.2, 1e-2)),
+               Simplex.alongAxes(OptimTestUtils.point(2, 0.2, 1e-2)),
                new PointValuePair(new double[] {f.xP, f.yP}, f.valueXpYp),
                1e-6);
     }
@@ -362,7 +365,7 @@ public class SimplexOptimizerNelderMeadTest {
                startPoint,
                goal,
                maxEvaluations,
-               Simplex.of(startPoint.length, 1),
+               Simplex.equalSidesAlongAxes(startPoint.length, 1d),
                expected,
                tol);
     }
@@ -383,6 +386,8 @@ public class SimplexOptimizerNelderMeadTest {
                         Simplex simplex,
                         PointValuePair expected,
                         double tol) {
+        final String name = func.toString();
+
         final int maxEval = Math.max(maxEvaluations, 12000);
         final SimplexOptimizer optim = new SimplexOptimizer(1e-13, 1e-14);
         final PointValuePair result = optim.optimize(new MaxEval(maxEval),
@@ -391,7 +396,6 @@ public class SimplexOptimizerNelderMeadTest {
                                                      new InitialGuess(startPoint),
                                                      simplex,
                                                      new NelderMeadTransform());
-        final String name = func.getClass().getSimpleName();
 
         final double[] endPoint = result.getPoint();
         final double funcValue = result.getValue();
diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexTest.java
new file mode 100644
index 0000000..260d251
--- /dev/null
+++ b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexTest.java
@@ -0,0 +1,199 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv;
+
+import java.util.Arrays;
+import java.util.List;
+import java.util.ArrayList;
+import org.junit.Assert;
+import org.junit.Test;
+import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
+import org.apache.commons.math4.legacy.optim.PointValuePair;
+
+/**
+ * Test for {@link Simplex}.
+ */
+public class SimplexTest {
+    @Test
+    public void testTriangle() {
+        final double[] o = new double[] { 1, 2 };
+        final double s = 3.45;
+        final int dim = o.length;
+        final int size = dim + 1;
+
+        final double[][] expected = new double[][] {
+            o,
+            new double[] { o[0] + s, o[1] },
+            new double[] { o[0] + s, o[1] + s }
+        };
+
+        final Simplex simplex = Simplex.equalSidesAlongAxes(dim, s).translate(o);
+        Assert.assertEquals(dim, simplex.getDimension());
+        Assert.assertEquals(size, simplex.getSize());
+
+        for (int i = 0; i < size; i++) {
+            final double[] e = expected[i];
+            final double[] a = simplex.get(i).getPoint();
+            Assert.assertTrue("i=" + i + ": e=" + Arrays.toString(e) + " a=" + Arrays.toString(a),
+                              Arrays.equals(e, a));
+        }
+    }
+
+    @Test
+    public void testJavadocExample() {
+        final double[] o = new double[] { 1, 10, 2 };
+        final double[] start = new double[] { 1, 1, 1 };
+        final int dim = o.length;
+        final int size = dim + 1;
+
+        final double[][] expected = new double[][] {
+            new double[] { 1, 1, 1 },
+            new double[] { 2, 1, 1 },
+            new double[] { 2, 11, 1 },
+            new double[] { 2, 11, 3 },
+        };
+
+        final Simplex simplex = Simplex.alongAxes(o).translate(start);
+        Assert.assertEquals(dim, simplex.getDimension());
+        Assert.assertEquals(size, simplex.getSize());
+
+        for (int i = 0; i < size; i++) {
+            final double[] e = expected[i];
+            final double[] a = simplex.get(i).getPoint();
+
+            final String msg = "i=" + i +
+                ": e=" + Arrays.toString(e) +
+                " a=" + Arrays.toString(a);
+            Assert.assertTrue(msg, Arrays.equals(e, a));
+        }
+    }
+
+    @Test
+    public void testCentroid() {
+        final List<PointValuePair> list = new ArrayList<>();
+
+        double[] centroid = null;
+
+        // Interval.
+        list.clear();
+        list.add(new PointValuePair(new double[] { 1, 1 }, 0d));
+        list.add(new PointValuePair(new double[] { 2, 2 }, 0d));
+        centroid = Simplex.centroid(list);
+        Assert.assertEquals(1.5, centroid[0], 0d);
+        Assert.assertEquals(1.5, centroid[1], 0d);
+
+        // Square.
+        list.clear();
+        list.add(new PointValuePair(new double[] { 1, 0 }, 0d));
+        list.add(new PointValuePair(new double[] { 2, 0 }, 0d));
+        list.add(new PointValuePair(new double[] { 2, 1 }, 0d));
+        list.add(new PointValuePair(new double[] { 1, 1 }, 0d));
+        centroid = Simplex.centroid(list);
+        Assert.assertEquals(1.5, centroid[0], 0d);
+        Assert.assertEquals(0.5, centroid[1], 0d);
+    }
+
+    @Test
+    public void testAsList() {
+        final int dim = 10;
+        final Simplex simplex = Simplex.equalSidesAlongAxes(dim, 1d);
+        final List<PointValuePair> list = simplex.asList();
+        Assert.assertEquals(dim + 1, list.size());
+
+        for (int i = 0; i < dim + 1; i++) {
+            final PointValuePair p = list.get(i);
+            final double[] a = simplex.get(i).getPointRef();
+            final double[] b = p.getPointRef();
+            Assert.assertFalse(a == b);
+            Assert.assertTrue(Arrays.equals(a, b));
+            Assert.assertTrue(Double.isNaN(p.getValue()));
+        }
+    }
+
+    @Test
+    public void testReplaceLast1() {
+        final int dim = 7;
+        final int nPoints = dim + 1;
+        final Simplex initSimplex = Simplex.equalSidesAlongAxes(dim, 1d);
+        for (PointValuePair pv : initSimplex.asList()) {
+            Assert.assertTrue(Double.isNaN(pv.getValue()));
+        }
+
+        final int nRepl = 3;
+        final double value = 12.345;
+        final List<PointValuePair> replace = new ArrayList<>();
+        for (int i = 0; i < nRepl; i++) {
+            replace.add(new PointValuePair(new double[dim], value));
+        }
+
+        final Simplex newSimplex = initSimplex.replaceLast(replace);
+        Assert.assertEquals(nPoints, newSimplex.getSize());
+
+        final int from = nPoints - nRepl;
+        for (int i = 0; i < from; i++) {
+            Assert.assertTrue(Double.isNaN(newSimplex.get(i).getValue()));
+        }
+        for (int i = from; i < nPoints; i++) {
+            Assert.assertEquals(value, newSimplex.get(i).getValue(), 0d);
+        }
+    }
+
+    @Test
+    public void testReplaceLast2() {
+        final int dim = 7;
+        final int nPoints = dim + 1;
+        final Simplex initSimplex = Simplex.equalSidesAlongAxes(dim, 1d);
+        for (PointValuePair pv : initSimplex.asList()) {
+            Assert.assertTrue(Double.isNaN(pv.getValue()));
+        }
+
+        final double value = 12.345;
+        final PointValuePair replace = new PointValuePair(new double[dim], value);
+        final Simplex newSimplex = initSimplex.replaceLast(replace);
+        Assert.assertEquals(nPoints, newSimplex.getSize());
+
+        final int from = nPoints - 1;
+        for (int i = 0; i < from; i++) {
+            Assert.assertTrue(Double.isNaN(newSimplex.get(i).getValue()));
+        }
+        for (int i = from; i < nPoints; i++) {
+            Assert.assertEquals(value, newSimplex.get(i).getValue(), 0d);
+        }
+    }
+
+    @Test
+    public void testNewPoint() {
+        final double[] a = new double[] { 1, 2 };
+        final double[] b = new double[] { 3, -1 };
+        final double s = 7;
+
+        final MultivariateFunction f = x -> {
+            double v = 0;
+            for (int i = 0; i < x.length; i++) {
+                v += x[i] * x[i];
+            }
+            return v;
+        };
+
+        final double[] eP = new double[] { 15, -19 };
+        final double eV = 586;
+
+        final PointValuePair pv = Simplex.newPoint(a, s, b, f);
+        Assert.assertTrue(Arrays.equals(eP, pv.getPoint()));
+        Assert.assertEquals(eV, pv.getValue(), 0d);
+    }
+}