Erik Ovegard
11 years ago
committed by
Erik Ovegard
Erik Ovegard
3 changed files with 342 additions and 1 deletions
@ -0,0 +1,339 @@ |
|||
using MathNet.Numerics.LinearAlgebra; |
|||
using System; |
|||
using System.Collections.Generic; |
|||
using System.Linq; |
|||
using System.Text; |
|||
|
|||
namespace MathNet.Numerics.Optimization |
|||
{ |
|||
public sealed class NelderMeadSimplex |
|||
{ |
|||
private static readonly double JITTER = 1e-10d; // a small value used to protect against floating point noise
|
|||
|
|||
public static MinimizationResult Regress(IObjectiveFunction objectiveFunction, Vector<double> initialGuess, |
|||
double convergenceTolerance, int maxEvaluations) |
|||
{ |
|||
SimplexConstant[] simplexConstants = SimplexConstant.CreateFromVector(initialGuess); |
|||
|
|||
// confirm that we are in a position to commence
|
|||
if (objectiveFunction == null) |
|||
throw new InvalidOperationException("ObjectiveFunction must be set to a valid ObjectiveFunctionDelegate"); |
|||
|
|||
if (simplexConstants == null) |
|||
throw new InvalidOperationException("SimplexConstants must be initialized"); |
|||
|
|||
// create the initial simplex
|
|||
int numDimensions = simplexConstants.Length; |
|||
int numVertices = numDimensions + 1; |
|||
Vector<double>[] vertices = _initializeVertices(simplexConstants); |
|||
double[] errorValues = new double[numVertices]; |
|||
|
|||
int evaluationCount = 0; |
|||
MinimizationResult.ExitCondition exitCondition = MinimizationResult.ExitCondition.None; |
|||
ErrorProfile errorProfile; |
|||
|
|||
errorValues = _initializeErrorValues(vertices, objectiveFunction); |
|||
|
|||
// iterate until we converge, or complete our permitted number of iterations
|
|||
while (true) |
|||
{ |
|||
errorProfile = _evaluateSimplex(errorValues); |
|||
|
|||
// see if the range in point heights is small enough to exit
|
|||
if (_hasConverged(convergenceTolerance, errorProfile, errorValues)) |
|||
{ |
|||
exitCondition = MinimizationResult.ExitCondition.Converged; |
|||
break; |
|||
} |
|||
|
|||
// attempt a reflection of the simplex
|
|||
double reflectionPointValue = _tryToScaleSimplex(-1.0, ref errorProfile, vertices, errorValues, objectiveFunction); |
|||
++evaluationCount; |
|||
if (reflectionPointValue <= errorValues[errorProfile.LowestIndex]) |
|||
{ |
|||
// it's better than the best point, so attempt an expansion of the simplex
|
|||
double expansionPointValue = _tryToScaleSimplex(2.0, ref errorProfile, vertices, errorValues, objectiveFunction); |
|||
++evaluationCount; |
|||
} |
|||
else if (reflectionPointValue >= errorValues[errorProfile.NextHighestIndex]) |
|||
{ |
|||
// it would be worse than the second best point, so attempt a contraction to look
|
|||
// for an intermediate point
|
|||
double currentWorst = errorValues[errorProfile.HighestIndex]; |
|||
double contractionPointValue = _tryToScaleSimplex(0.5, ref errorProfile, vertices, errorValues, objectiveFunction); |
|||
++evaluationCount; |
|||
if (contractionPointValue >= currentWorst) |
|||
{ |
|||
// that would be even worse, so let's try to contract uniformly towards the low point;
|
|||
// don't bother to update the error profile, we'll do it at the start of the
|
|||
// next iteration
|
|||
_shrinkSimplex(errorProfile, vertices, errorValues, objectiveFunction); |
|||
evaluationCount += numVertices; // that required one function evaluation for each vertex; keep track
|
|||
} |
|||
} |
|||
// check to see if we have exceeded our alloted number of evaluations
|
|||
if (evaluationCount >= maxEvaluations) |
|||
{ |
|||
exitCondition = MinimizationResult.ExitCondition.LackOfProgress; |
|||
break; |
|||
} |
|||
} |
|||
var regressionResult = new MinimizationResult(null, evaluationCount, exitCondition); |
|||
return regressionResult; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Evaluate the objective function at each vertex to create a corresponding
|
|||
/// list of error values for each vertex
|
|||
/// </summary>
|
|||
/// <param name="vertices"></param>
|
|||
/// <returns></returns>
|
|||
private static double[] _initializeErrorValues(Vector<double>[] vertices, IObjectiveFunction objectiveFunction) |
|||
{ |
|||
double[] errorValues = new double[vertices.Length]; |
|||
for (int i = 0; i < vertices.Length; i++) |
|||
{ |
|||
objectiveFunction.EvaluateAt(vertices[i]); |
|||
errorValues[i] = objectiveFunction.Value; |
|||
} |
|||
return errorValues; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Check whether the points in the error profile have so little range that we
|
|||
/// consider ourselves to have converged
|
|||
/// </summary>
|
|||
/// <param name="errorProfile"></param>
|
|||
/// <param name="errorValues"></param>
|
|||
/// <returns></returns>
|
|||
private static bool _hasConverged(double convergenceTolerance, ErrorProfile errorProfile, double[] errorValues) |
|||
{ |
|||
double range = 2 * Math.Abs(errorValues[errorProfile.HighestIndex] - errorValues[errorProfile.LowestIndex]) / |
|||
(Math.Abs(errorValues[errorProfile.HighestIndex]) + Math.Abs(errorValues[errorProfile.LowestIndex]) + JITTER); |
|||
|
|||
if (range < convergenceTolerance) |
|||
{ |
|||
return true; |
|||
} |
|||
else |
|||
{ |
|||
return false; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Examine all error values to determine the ErrorProfile
|
|||
/// </summary>
|
|||
/// <param name="errorValues"></param>
|
|||
/// <returns></returns>
|
|||
private static ErrorProfile _evaluateSimplex(double[] errorValues) |
|||
{ |
|||
ErrorProfile errorProfile = new ErrorProfile(); |
|||
if (errorValues[0] > errorValues[1]) |
|||
{ |
|||
errorProfile.HighestIndex = 0; |
|||
errorProfile.NextHighestIndex = 1; |
|||
} |
|||
else |
|||
{ |
|||
errorProfile.HighestIndex = 1; |
|||
errorProfile.NextHighestIndex = 0; |
|||
} |
|||
|
|||
for (int index = 0; index < errorValues.Length; index++) |
|||
{ |
|||
double errorValue = errorValues[index]; |
|||
if (errorValue <= errorValues[errorProfile.LowestIndex]) |
|||
{ |
|||
errorProfile.LowestIndex = index; |
|||
} |
|||
if (errorValue > errorValues[errorProfile.HighestIndex]) |
|||
{ |
|||
errorProfile.NextHighestIndex = errorProfile.HighestIndex; // downgrade the current highest to next highest
|
|||
errorProfile.HighestIndex = index; |
|||
} |
|||
else if (errorValue > errorValues[errorProfile.NextHighestIndex] && index != errorProfile.HighestIndex) |
|||
{ |
|||
errorProfile.NextHighestIndex = index; |
|||
} |
|||
} |
|||
|
|||
return errorProfile; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Construct an initial simplex, given starting guesses for the constants, and
|
|||
/// initial step sizes for each dimension
|
|||
/// </summary>
|
|||
/// <param name="simplexConstants"></param>
|
|||
/// <returns></returns>
|
|||
private static Vector<double>[] _initializeVertices(SimplexConstant[] simplexConstants) |
|||
{ |
|||
int numDimensions = simplexConstants.Length; |
|||
Vector<double>[] vertices = new Vector<double>[numDimensions + 1]; |
|||
|
|||
// define one point of the simplex as the given initial guesses
|
|||
var p0 = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(numDimensions); |
|||
for (int i = 0; i < numDimensions; i++) |
|||
{ |
|||
p0[i] = simplexConstants[i].Value; |
|||
} |
|||
|
|||
// now fill in the vertices, creating the additional points as:
|
|||
// P(i) = P(0) + Scale(i) * UnitVector(i)
|
|||
vertices[0] = p0; |
|||
for (int i = 0; i < numDimensions; i++) |
|||
{ |
|||
double scale = simplexConstants[i].InitialPerturbation; |
|||
Vector<double> unitVector = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(numDimensions); |
|||
unitVector[i] = 1; |
|||
vertices[i + 1] = p0.Add(unitVector.Multiply(scale)); |
|||
} |
|||
return vertices; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Test a scaling operation of the high point, and replace it if it is an improvement
|
|||
/// </summary>
|
|||
/// <param name="scaleFactor"></param>
|
|||
/// <param name="errorProfile"></param>
|
|||
/// <param name="vertices"></param>
|
|||
/// <param name="errorValues"></param>
|
|||
/// <returns></returns>
|
|||
private static double _tryToScaleSimplex(double scaleFactor, ref ErrorProfile errorProfile, Vector<double>[] vertices, |
|||
double[] errorValues, IObjectiveFunction objectiveFunction) |
|||
{ |
|||
// find the centroid through which we will reflect
|
|||
Vector<double> centroid = _computeCentroid(vertices, errorProfile); |
|||
|
|||
// define the vector from the centroid to the high point
|
|||
Vector<double> centroidToHighPoint = vertices[errorProfile.HighestIndex].Subtract(centroid); |
|||
|
|||
// scale and position the vector to determine the new trial point
|
|||
Vector<double> newPoint = centroidToHighPoint.Multiply(scaleFactor).Add(centroid); |
|||
|
|||
// evaluate the new point
|
|||
objectiveFunction.EvaluateAt(newPoint); |
|||
double newErrorValue = objectiveFunction.Value; |
|||
|
|||
// if it's better, replace the old high point
|
|||
if (newErrorValue < errorValues[errorProfile.HighestIndex]) |
|||
{ |
|||
vertices[errorProfile.HighestIndex] = newPoint; |
|||
errorValues[errorProfile.HighestIndex] = newErrorValue; |
|||
} |
|||
|
|||
return newErrorValue; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Contract the simplex uniformly around the lowest point
|
|||
/// </summary>
|
|||
/// <param name="errorProfile"></param>
|
|||
/// <param name="vertices"></param>
|
|||
/// <param name="errorValues"></param>
|
|||
private static void _shrinkSimplex(ErrorProfile errorProfile, Vector<double>[] vertices, double[] errorValues, |
|||
IObjectiveFunction objectiveFunction) |
|||
{ |
|||
Vector<double> lowestVertex = vertices[errorProfile.LowestIndex]; |
|||
for (int i = 0; i < vertices.Length; i++) |
|||
{ |
|||
if (i != errorProfile.LowestIndex) |
|||
{ |
|||
vertices[i] = (vertices[i].Add(lowestVertex)).Multiply(0.5); |
|||
objectiveFunction.EvaluateAt(vertices[i]); |
|||
errorValues[i] = objectiveFunction.Value; |
|||
} |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Compute the centroid of all points except the worst
|
|||
/// </summary>
|
|||
/// <param name="vertices"></param>
|
|||
/// <param name="errorProfile"></param>
|
|||
/// <returns></returns>
|
|||
private static Vector<double> _computeCentroid(Vector<double>[] vertices, ErrorProfile errorProfile) |
|||
{ |
|||
int numVertices = vertices.Length; |
|||
// find the centroid of all points except the worst one
|
|||
Vector<double> centroid = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(numVertices - 1); |
|||
for (int i = 0; i < numVertices; i++) |
|||
{ |
|||
if (i != errorProfile.HighestIndex) |
|||
{ |
|||
centroid = centroid.Add(vertices[i]); |
|||
} |
|||
} |
|||
return centroid.Multiply(1.0d / (numVertices - 1)); |
|||
} |
|||
|
|||
private sealed class SimplexConstant |
|||
{ |
|||
private double _value; |
|||
private double _initialPerturbation; |
|||
|
|||
public SimplexConstant(double value, double initialPerturbation) |
|||
{ |
|||
_value = value; |
|||
_initialPerturbation = initialPerturbation; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// The value of the constant
|
|||
/// </summary>
|
|||
public double Value |
|||
{ |
|||
get { return _value; } |
|||
set { _value = value; } |
|||
} |
|||
|
|||
// The size of the initial perturbation
|
|||
public double InitialPerturbation |
|||
{ |
|||
get { return _initialPerturbation; } |
|||
set { _initialPerturbation = value; } |
|||
} |
|||
|
|||
public static SimplexConstant[] CreateFromVector(Vector<double> initialGuess) |
|||
{ |
|||
var constants = new SimplexConstant[initialGuess.Count]; |
|||
|
|||
for (int i = 0; i < constants.Length;i++ ) |
|||
{ |
|||
double pertubation = initialGuess[i]==0.0 ? 1e-5 : initialGuess[i]*1e-5; |
|||
constants[i] = new SimplexConstant(initialGuess[i], pertubation); |
|||
} |
|||
return constants; |
|||
} |
|||
} |
|||
|
|||
private sealed class ErrorProfile |
|||
{ |
|||
private int _highestIndex; |
|||
private int _nextHighestIndex; |
|||
private int _lowestIndex; |
|||
|
|||
public int HighestIndex |
|||
{ |
|||
get { return _highestIndex; } |
|||
set { _highestIndex = value; } |
|||
} |
|||
|
|||
public int NextHighestIndex |
|||
{ |
|||
get { return _nextHighestIndex; } |
|||
set { _nextHighestIndex = value; } |
|||
} |
|||
|
|||
public int LowestIndex |
|||
{ |
|||
get { return _lowestIndex; } |
|||
set { _lowestIndex = value; } |
|||
} |
|||
} |
|||
} |
|||
|
|||
|
|||
|
|||
} |
|||
Loading…
Reference in new issue