diff --git a/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs b/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs index 9f403286..793113a2 100644 --- a/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs +++ b/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs @@ -1720,7 +1720,8 @@ namespace MathNet.Numerics.LinearAlgebra /// matrix and a given other matrix being the 'x' of atan2 and the /// 'this' matrix being the 'y' /// - /// + /// The other matrix 'y' + /// The matrix with the result and 'x' /// public void PointwiseAtan2(Matrix other, Matrix result) { diff --git a/src/Numerics/LinearAlgebra/Vector.Arithmetic.cs b/src/Numerics/LinearAlgebra/Vector.Arithmetic.cs index cbb9bc03..edc9cb7e 100644 --- a/src/Numerics/LinearAlgebra/Vector.Arithmetic.cs +++ b/src/Numerics/LinearAlgebra/Vector.Arithmetic.cs @@ -1028,6 +1028,7 @@ namespace MathNet.Numerics.LinearAlgebra /// /// Function which takes a scalar and a vector, modifies the vector in place and returns void /// The scalar to be passed to the function + /// The vector where the result will be placed /// If this vector and are not the same size. protected void PointwiseBinary(Action> f, T x, Vector result) { diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 9f218718..2eb58dfa 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -44,7 +44,7 @@ 4 AllRules.ruleset 1591 - 5 + 6 ..\..\out\lib-debug\Net40\ @@ -58,7 +58,7 @@ prompt 4 AllRules.ruleset - 5 + 6 ..\..\out\lib-signed\Net40\ @@ -75,7 +75,7 @@ 4 AllRules.ruleset 1591 - 5 + 6 @@ -115,6 +115,24 @@ + + + + + + + + + + + + + + + + + + @@ -233,6 +251,7 @@ + @@ -242,6 +261,17 @@ + + + + + + + + + + + diff --git a/src/Numerics/Optimization/BfgsBMinimizer.cs b/src/Numerics/Optimization/BfgsBMinimizer.cs new file mode 100644 index 00000000..9281a589 --- /dev/null +++ b/src/Numerics/Optimization/BfgsBMinimizer.cs @@ -0,0 +1,305 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2016 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; +using System.Collections.Generic; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.LinearAlgebra.Double; +using MathNet.Numerics.Optimization.LineSearch; + +namespace MathNet.Numerics.Optimization +{ + /// + /// Broyden–Fletcher–Goldfarb–Shanno Bounded (BFGS-B) algorithm is an iterative method for solving box-constrained nonlinear optimization problems + /// http://www.ece.northwestern.edu/~nocedal/PSfiles/limited.ps.gz + /// + public class BfgsBMinimizer : BfgsMinimizerBase + { + public BfgsBMinimizer(double gradientTolerance, double parameterTolerance, double functionProgressTolerance, int maximumIterations = 1000) + : base(gradientTolerance,parameterTolerance,functionProgressTolerance,maximumIterations) + { + } + + /// + /// Find the minimum of the objective function given lower and upper bounds + /// + /// The objective function, must support a gradient + /// The lower bound + /// The upper bound + /// The initial guess + /// The MinimizationResult which contains the minimum and the ExitCondition + public MinimizationResult FindMinimum(IObjectiveFunction objective, Vector lowerBound, Vector upperBound, Vector initialGuess) + { + _lowerBound = lowerBound; + _upperBound = upperBound; + if (!objective.IsGradientSupported) + throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for BFGS minimization."); + + // Check that dimensions match + if (lowerBound.Count != upperBound.Count || lowerBound.Count != initialGuess.Count) + throw new ArgumentException("Dimensions of bounds and/or initial guess do not match."); + + // Check that initial guess is feasible + for (int ii = 0; ii < initialGuess.Count; ++ii) + if (initialGuess[ii] < lowerBound[ii] || initialGuess[ii] > upperBound[ii]) + throw new ArgumentException("Initial guess is not in the feasible region"); + + objective.EvaluateAt(initialGuess); + ValidateGradientAndObjective(objective); + + // Check that we're not already done + MinimizationResult.ExitCondition currentExitCondition = ExitCriteriaSatisfied(objective, null, 0); + if (currentExitCondition != MinimizationResult.ExitCondition.None) + return new MinimizationResult(objective, 0, currentExitCondition); + + // Set up line search algorithm + var lineSearcher = new StrongWolfeLineSearch(1e-4, 0.9, Math.Max(ParameterTolerance, 1e-5), maxIterations: 1000); + + // Declare state variables + Vector reducedSolution1, reducedGradient, reducedInitialPoint, reducedCauchyPoint, solution1; + Matrix reducedHessian; + List reducedMap; + + // First step + var pseudoHessian = CreateMatrix.DiagonalIdentity(initialGuess.Count); + + // Determine active set + var gradientProjectionResult = QuadraticGradientProjectionSearch.Search(objective.Point, objective.Gradient, pseudoHessian, lowerBound, upperBound); + var cauchyPoint = gradientProjectionResult.CauchyPoint; + var fixedCount = gradientProjectionResult.FixedCount; + var isFixed = gradientProjectionResult.IsFixed; + var freeCount = lowerBound.Count - fixedCount; + + if (freeCount > 0) + { + reducedGradient = new DenseVector(freeCount); + reducedHessian = new DenseMatrix(freeCount, freeCount); + reducedMap = new List(freeCount); + reducedInitialPoint = new DenseVector(freeCount); + reducedCauchyPoint = new DenseVector(freeCount); + + CreateReducedData(objective.Point, cauchyPoint, isFixed, lowerBound, upperBound, objective.Gradient, pseudoHessian, reducedInitialPoint, reducedCauchyPoint, reducedGradient, reducedHessian, reducedMap); + + // Determine search direction and maximum step size + reducedSolution1 = reducedInitialPoint + reducedHessian.Cholesky().Solve(-reducedGradient); + + solution1 = ReducedToFull(reducedMap, reducedSolution1, cauchyPoint); + } + else + { + solution1 = cauchyPoint; + } + + var directionFromCauchy = solution1 - cauchyPoint; + var maxStepFromCauchyPoint = FindMaxStep(cauchyPoint, directionFromCauchy, lowerBound, upperBound); + + var solution2 = cauchyPoint + Math.Min(maxStepFromCauchyPoint, 1.0)*directionFromCauchy; + + var lineSearchDirection = solution2 - objective.Point; + var maxLineSearchStep = FindMaxStep(objective.Point, lineSearchDirection, lowerBound, upperBound); + var estStepSize = -objective.Gradient*lineSearchDirection/(lineSearchDirection*pseudoHessian*lineSearchDirection); + + var startingStepSize = Math.Min(Math.Max(estStepSize, 1.0), maxLineSearchStep); + + // Line search + LineSearchResult lineSearchResult; + try + { + lineSearchResult = lineSearcher.FindConformingStep(objective, lineSearchDirection, startingStepSize, upperBound: maxLineSearchStep); + } + catch (Exception e) + { + throw new InnerOptimizationException("Line search failed.", e); + } + + var previousPoint = objective.Fork(); + var candidatePoint = lineSearchResult.FunctionInfoAtMinimum; + var gradient = candidatePoint.Gradient; + var step = candidatePoint.Point - initialGuess; + + // Subsequent steps + int totalLineSearchSteps = lineSearchResult.Iterations; + int iterationsWithNontrivialLineSearch = lineSearchResult.Iterations > 0 ? 0 : 1; + + int iterations = DoBfgsUpdate(ref currentExitCondition, lineSearcher, ref pseudoHessian, ref lineSearchDirection, ref previousPoint, ref lineSearchResult, ref candidatePoint, ref step, ref totalLineSearchSteps, ref iterationsWithNontrivialLineSearch); + + if (iterations == MaximumIterations && currentExitCondition == MinimizationResult.ExitCondition.None) + throw new MaximumIterationsException(string.Format("Maximum iterations ({0}) reached.", MaximumIterations)); + + return new MinimizationWithLineSearchResult(candidatePoint, iterations, currentExitCondition, totalLineSearchSteps, iterationsWithNontrivialLineSearch); + } + + protected override Vector CalculateSearchDirection(ref Matrix pseudoHessian, + out double maxLineSearchStep, + out double startingStepSize, + IObjectiveFunction previousPoint, + IObjectiveFunction candidatePoint, + Vector step) + { + Vector lineSearchDirection; + var y = candidatePoint.Gradient - previousPoint.Gradient; + + double sy = step * y; + if (sy > 0.0) // only do update if it will create a positive definite matrix + { + double sts = step * step; + + var Hs = pseudoHessian * step; + var sHs = step * pseudoHessian * step; + pseudoHessian = pseudoHessian + y.OuterProduct(y) * (1.0 / sy) - Hs.OuterProduct(Hs) * (1.0 / sHs); + } + else + { + //pseudo_hessian = LinearAlgebra.Double.DiagonalMatrix.Identity(initial_guess.Count); + } + + // Determine active set + var gradientProjectionResult = QuadraticGradientProjectionSearch.Search(candidatePoint.Point, candidatePoint.Gradient, pseudoHessian, _lowerBound, _upperBound); + var cauchyPoint = gradientProjectionResult.CauchyPoint; + var fixedCount = gradientProjectionResult.FixedCount; + var isFixed = gradientProjectionResult.IsFixed; + var freeCount = _lowerBound.Count - fixedCount; + Vector solution1; + if (freeCount > 0) + { + var reducedGradient = new DenseVector(freeCount); + var reducedHessian = new DenseMatrix(freeCount, freeCount); + var reducedMap = new List(freeCount); + var reducedInitialPoint = new DenseVector(freeCount); + var reducedCauchyPoint = new DenseVector(freeCount); + + CreateReducedData(candidatePoint.Point, cauchyPoint, isFixed, _lowerBound, _upperBound, candidatePoint.Gradient, pseudoHessian, reducedInitialPoint, reducedCauchyPoint, reducedGradient, reducedHessian, reducedMap); + + // Determine search direction and maximum step size + Vector reducedSolution1 = reducedInitialPoint + reducedHessian.Cholesky().Solve(-reducedGradient); + + solution1 = ReducedToFull(reducedMap, reducedSolution1, cauchyPoint); + } + else + { + solution1 = cauchyPoint; + } + + var directionFromCauchy = solution1 - cauchyPoint; + var maxStepFromCauchyPoint = FindMaxStep(cauchyPoint, directionFromCauchy, _lowerBound, _upperBound); + + var solution2 = cauchyPoint + Math.Min(maxStepFromCauchyPoint, 1.0) * directionFromCauchy; + + lineSearchDirection = solution2 - candidatePoint.Point; + maxLineSearchStep = FindMaxStep(candidatePoint.Point, lineSearchDirection, _lowerBound, _upperBound); + + if (maxLineSearchStep == 0.0) + { + lineSearchDirection = cauchyPoint - candidatePoint.Point; + maxLineSearchStep = FindMaxStep(candidatePoint.Point, lineSearchDirection, _lowerBound, _upperBound); + } + + double estStepSize = -candidatePoint.Gradient * lineSearchDirection / (lineSearchDirection * pseudoHessian * lineSearchDirection); + + startingStepSize = Math.Min(Math.Max(estStepSize, 1.0), maxLineSearchStep); + return lineSearchDirection; + } + + private static Vector ReducedToFull(List reducedMap, Vector reducedVector, Vector fullVector) + { + var output = fullVector.Clone(); + for (int ii = 0; ii < reducedMap.Count; ++ii) + output[reducedMap[ii]] = reducedVector[ii]; + return output; + } + + private Vector _lowerBound; + private Vector _upperBound; + + private static double FindMaxStep(Vector startingPoint, Vector searchDirection, Vector lowerBound, Vector upperBound) + { + double maxStep = Double.PositiveInfinity; + for (int ii = 0; ii < startingPoint.Count; ++ii) + { + double paramMaxStep; + if (searchDirection[ii] > 0) + paramMaxStep = (upperBound[ii] - startingPoint[ii])/searchDirection[ii]; + else if (searchDirection[ii] < 0) + paramMaxStep = (startingPoint[ii] - lowerBound[ii])/-searchDirection[ii]; + else + paramMaxStep = Double.PositiveInfinity; + + if (paramMaxStep < maxStep) + maxStep = paramMaxStep; + } + return maxStep; + } + + private static void CreateReducedData(Vector initialPoint, Vector cauchyPoint, List isFixed, Vector lowerBound, Vector upperBound, Vector gradient, Matrix pseudoHessian, Vector reducedInitialPoint, Vector reducedCauchyPoint, Vector reducedGradient, Matrix reducedHessian, List reducedMap) + { + int ll = 0; + for (int ii = 0; ii < lowerBound.Count; ++ii) + { + if (!isFixed[ii]) + { + // hessian + int mm = 0; + for (int jj = 0; jj < lowerBound.Count; ++jj) + { + if (!isFixed[jj]) + { + reducedHessian[ll, mm++] = pseudoHessian[ii, jj]; + } + } + + // gradient + reducedInitialPoint[ll] = initialPoint[ii]; + reducedCauchyPoint[ll] = cauchyPoint[ii]; + reducedGradient[ll] = gradient[ii]; + ll += 1; + reducedMap.Add(ii); + + } + } + } + + protected override double GetProjectedGradient(IObjectiveFunction candidatePoint, int ii) + { + double projectedGradient; + bool atLowerBound = candidatePoint.Point[ii] - _lowerBound[ii] < VerySmall; + bool atUpperBound = _upperBound[ii] - candidatePoint.Point[ii] < VerySmall; + + if (atLowerBound && atUpperBound) + projectedGradient = 0.0; + else if (atLowerBound) + projectedGradient = Math.Min(candidatePoint.Gradient[ii], 0.0); + else if (atUpperBound) + projectedGradient = Math.Max(candidatePoint.Gradient[ii], 0.0); + else + projectedGradient = base.GetProjectedGradient(candidatePoint, ii); + return projectedGradient; + } + } +} diff --git a/src/Numerics/Optimization/BfgsMinimizer.cs b/src/Numerics/Optimization/BfgsMinimizer.cs new file mode 100644 index 00000000..afac00f7 --- /dev/null +++ b/src/Numerics/Optimization/BfgsMinimizer.cs @@ -0,0 +1,142 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2016 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.Optimization.LineSearch; + +namespace MathNet.Numerics.Optimization +{ + /// + /// Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm is an iterative method for solving unconstrained nonlinear optimization problems + /// + public class BfgsMinimizer : BfgsMinimizerBase + { + /// + /// Creates BFGS minimizer + /// + /// The gradient tolerance + /// The parameter tolerance + /// The funciton progress tolerance + /// The maximum number of iterations + public BfgsMinimizer(double gradientTolerance, double parameterTolerance, double functionProgressTolerance, int maximumIterations=1000) + :base(gradientTolerance,parameterTolerance,functionProgressTolerance,maximumIterations) + { + } + + /// + /// Find the minimum of the objective function given lower and upper bounds + /// + /// The objective function, must support a gradient + /// The initial guess + /// The MinimizationResult which contains the minimum and the ExitCondition + public MinimizationResult FindMinimum(IObjectiveFunction objective, Vector initialGuess) + { + if (!objective.IsGradientSupported) + throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for BFGS minimization."); + + objective.EvaluateAt(initialGuess); + ValidateGradientAndObjective(objective); + + // Check that we're not already done + MinimizationResult.ExitCondition currentExitCondition = ExitCriteriaSatisfied(objective, null, 0); + if (currentExitCondition != MinimizationResult.ExitCondition.None) + return new MinimizationResult(objective, 0, currentExitCondition); + + // Set up line search algorithm + var lineSearcher = new WeakWolfeLineSearch(1e-4, 0.9, Math.Max(ParameterTolerance, 1e-10), 1000); + + // First step + var inversePseudoHessian = CreateMatrix.DenseIdentity(initialGuess.Count); + var lineSearchDirection = -objective.Gradient; + var stepSize = 100 * GradientTolerance / (lineSearchDirection * lineSearchDirection); + + var previousPoint = objective; + + LineSearchResult lineSearchResult; + try + { + lineSearchResult = lineSearcher.FindConformingStep(objective, lineSearchDirection, stepSize); + } + catch (OptimizationException e) + { + throw new InnerOptimizationException("Line search failed.", e); + } + catch (ArgumentException e) + { + throw new InnerOptimizationException("Line search failed.", e); + } + + var candidate = lineSearchResult.FunctionInfoAtMinimum; + ValidateGradientAndObjective(candidate); + + var gradient = candidate.Gradient; + var step = candidate.Point - initialGuess; + + // Subsequent steps + Matrix I = CreateMatrix.DiagonalIdentity(initialGuess.Count); + int iterations; + int totalLineSearchSteps = lineSearchResult.Iterations; + int iterationsWithNontrivialLineSearch = lineSearchResult.Iterations > 0 ? 0 : 1; + iterations = DoBfgsUpdate(ref currentExitCondition, lineSearcher, ref inversePseudoHessian, ref lineSearchDirection, ref previousPoint, ref lineSearchResult, ref candidate, ref step, ref totalLineSearchSteps, ref iterationsWithNontrivialLineSearch); + + if (iterations == MaximumIterations && currentExitCondition == MinimizationResult.ExitCondition.None) + throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations)); + + return new MinimizationWithLineSearchResult(candidate, iterations, MinimizationResult.ExitCondition.AbsoluteGradient, totalLineSearchSteps, iterationsWithNontrivialLineSearch); + } + + protected override Vector CalculateSearchDirection(ref Matrix inversePseudoHessian, + out double maxLineSearchStep, + out double startingStepSize, + IObjectiveFunction previousPoint, + IObjectiveFunction candidate, + Vector step) + { + startingStepSize = 1.0; + maxLineSearchStep = double.PositiveInfinity; + + Vector lineSearchDirection; + var y = candidate.Gradient - previousPoint.Gradient; + + double sy = step * y; + inversePseudoHessian = inversePseudoHessian + ((sy + y * inversePseudoHessian * y) / Math.Pow(sy, 2.0)) * step.OuterProduct(step) - ((inversePseudoHessian * y.ToColumnMatrix()) * step.ToRowMatrix() + step.ToColumnMatrix() * (y.ToRowMatrix() * inversePseudoHessian)) * (1.0 / sy); + lineSearchDirection = -inversePseudoHessian * candidate.Gradient; + + if (lineSearchDirection * candidate.Gradient >= 0.0) + { + lineSearchDirection = -candidate.Gradient; + inversePseudoHessian = CreateMatrix.DenseIdentity(candidate.Point.Count); + } + + return lineSearchDirection; + } + } +} diff --git a/src/Numerics/Optimization/BfgsMinimizerBase.cs b/src/Numerics/Optimization/BfgsMinimizerBase.cs new file mode 100644 index 00000000..16b06490 --- /dev/null +++ b/src/Numerics/Optimization/BfgsMinimizerBase.cs @@ -0,0 +1,158 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2016 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.LinearAlgebra.Double; +using MathNet.Numerics.Optimization.LineSearch; +using System; + +namespace MathNet.Numerics.Optimization +{ + public abstract class BfgsMinimizerBase + { + public double GradientTolerance { get; set; } + public double ParameterTolerance { get; set; } + public double FunctionProgressTolerance { get; set; } + public int MaximumIterations { get; set; } + + protected const double VerySmall = 1e-15; + + /// + /// Creates a base class for BFGS minimization + /// + /// The gradient tolerance + /// The parameter tolerance + /// The funciton progress tolerance + /// The maximum number of iterations + public BfgsMinimizerBase(double gradientTolerance, double parameterTolerance, double functionProgressTolerance, int maximumIterations) + { + GradientTolerance = gradientTolerance; + ParameterTolerance = parameterTolerance; + FunctionProgressTolerance = functionProgressTolerance; + MaximumIterations = maximumIterations; + } + + protected MinimizationResult.ExitCondition ExitCriteriaSatisfied(IObjectiveFunction candidatePoint, IObjectiveFunction lastPoint, int iterations) + { + Vector relGrad = new DenseVector(candidatePoint.Point.Count); + double relativeGradient = 0.0; + double normalizer = Math.Max(Math.Abs(candidatePoint.Value), 1.0); + for (int ii = 0; ii < relGrad.Count; ++ii) + { + double projectedGradient = GetProjectedGradient(candidatePoint, ii); + + double tmp = projectedGradient * + Math.Max(Math.Abs(candidatePoint.Point[ii]), 1.0) / normalizer; + relativeGradient = Math.Max(relativeGradient, Math.Abs(tmp)); + } + if (relativeGradient < GradientTolerance) + { + return MinimizationResult.ExitCondition.RelativeGradient; + } + + if (lastPoint != null) + { + double mostProgress = 0.0; + for (int ii = 0; ii < candidatePoint.Point.Count; ++ii) + { + var tmp = Math.Abs(candidatePoint.Point[ii] - lastPoint.Point[ii]) / + Math.Max(Math.Abs(lastPoint.Point[ii]), 1.0); + mostProgress = Math.Max(mostProgress, tmp); + } + if (mostProgress < ParameterTolerance) + { + return MinimizationResult.ExitCondition.LackOfProgress; + } + + double functionChange = candidatePoint.Value - lastPoint.Value; + if (iterations > 500 && functionChange < 0 && Math.Abs(functionChange) < FunctionProgressTolerance) + return MinimizationResult.ExitCondition.LackOfProgress; + } + + return MinimizationResult.ExitCondition.None; + } + + protected virtual double GetProjectedGradient(IObjectiveFunction candidatePoint, int ii) + { + return candidatePoint.Gradient[ii]; + } + + protected void ValidateGradientAndObjective(IObjectiveFunction eval) + { + foreach (var x in eval.Gradient) + { + if (Double.IsNaN(x) || Double.IsInfinity(x)) + throw new EvaluationException("Non-finite gradient returned.", eval); + } + if (Double.IsNaN(eval.Value) || Double.IsInfinity(eval.Value)) + throw new EvaluationException("Non-finite objective function returned.", eval); + } + + protected int DoBfgsUpdate(ref MinimizationResult.ExitCondition currentExitCondition, WolfeLineSearch lineSearcher, ref Matrix inversePseudoHessian, ref Vector lineSearchDirection, ref IObjectiveFunction previousPoint, ref LineSearchResult lineSearchResult, ref IObjectiveFunction candidate, ref Vector step, ref int totalLineSearchSteps, ref int iterationsWithNontrivialLineSearch) + { + int iterations; + for (iterations = 1; iterations < MaximumIterations; ++iterations) + { + double startingStepSize; + double maxLineSearchStep; + lineSearchDirection = CalculateSearchDirection(ref inversePseudoHessian, out maxLineSearchStep, out startingStepSize, previousPoint, candidate, step); + + try + { + lineSearchResult = lineSearcher.FindConformingStep(candidate, lineSearchDirection, startingStepSize, maxLineSearchStep); + } + catch (Exception e) + { + throw new InnerOptimizationException("Line search failed.", e); + } + + iterationsWithNontrivialLineSearch += lineSearchResult.Iterations > 0 ? 1 : 0; + totalLineSearchSteps += lineSearchResult.Iterations; + + step = lineSearchResult.FunctionInfoAtMinimum.Point - candidate.Point; + previousPoint = candidate; + candidate = lineSearchResult.FunctionInfoAtMinimum; + + currentExitCondition = ExitCriteriaSatisfied(candidate, previousPoint, iterations); + if (currentExitCondition != MinimizationResult.ExitCondition.None) + break; + } + + return iterations; + } + + protected abstract Vector CalculateSearchDirection(ref Matrix inversePseudoHessian, + out double maxLineSearchStep, + out double startingStepSize, + IObjectiveFunction previousPoint, + IObjectiveFunction candidate, + Vector step); + } +} diff --git a/src/Numerics/Optimization/BfgsSolver.cs b/src/Numerics/Optimization/BfgsSolver.cs new file mode 100644 index 00000000..d02c7583 --- /dev/null +++ b/src/Numerics/Optimization/BfgsSolver.cs @@ -0,0 +1,110 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2015 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.LinearAlgebra.Double; +using MathNet.Numerics.Optimization.LineSearch; + +namespace MathNet.Numerics.Optimization +{ + /// + /// Broyden-Fletcher-Goldfarb-Shanno solver for finding function minima + /// See http://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm + /// Inspired by implementation: https://github.com/PatWie/CppNumericalSolvers/blob/master/src/BfgsSolver.cpp + /// + public static class BfgsSolver + { + private const double GradientTolerance = 1e-5; + private const int MaxIterations = 100000; + + /// + /// Finds a minimum of a function by the BFGS quasi-Newton method + /// This uses the function and it's gradient (partial derivatives in each direction) and approximates the Hessian + /// + /// An initial guess + /// Evaluates the function at a point + /// Evaluates the gradient of the function at a point + /// The minimum found + public static Vector Solve(Vector initialGuess, Func, double> functionValue, Func, Vector> functionGradient) + { + var objectiveFunction = ObjectiveFunction.Gradient(functionValue, functionGradient); + objectiveFunction.EvaluateAt(initialGuess); + + int dim = initialGuess.Count; + int iter = 0; + // H represents the approximation of the inverse hessian matrix + // it is updated via the Sherman–Morrison formula (http://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula) + Matrix H = DenseMatrix.CreateIdentity(dim); + + Vector x = initialGuess; + Vector x_old = x; + Vector grad; + WolfeLineSearch wolfeLineSearch = new WeakWolfeLineSearch(1e-4, 0.9, 1e-5, 200); + do + { + // search along the direction of the gradient + grad = objectiveFunction.Gradient; + Vector p = -1 * H * grad; + var lineSearchResult = wolfeLineSearch.FindConformingStep(objectiveFunction, p, 1.0); + double rate = lineSearchResult.FinalStep; + x = x + rate * p; + Vector grad_old = grad; + + // update the gradient + objectiveFunction.EvaluateAt(x); + grad = objectiveFunction.Gradient;// functionGradient(x); + + Vector s = x - x_old; + Vector y = grad - grad_old; + + double rho = 1.0 / (y * s); + if (iter == 0) + { + // set up an initial hessian + H = (y * s) / (y * y) * DenseMatrix.CreateIdentity(dim); + } + + var sM = s.ToColumnMatrix(); + var yM = y.ToColumnMatrix(); + + // Update the estimate of the hessian + H = H + - rho * (sM * (yM.TransposeThisAndMultiply(H)) + (H * yM).TransposeAndMultiply(sM)) + + rho * rho * (y.DotProduct(H * y) + 1.0 / rho) * (sM.TransposeAndMultiply(sM)); + x_old = x; + iter++; + } + while ((grad.InfinityNorm() > GradientTolerance) && (iter < MaxIterations)); + + return x; + } + } +} diff --git a/src/Numerics/Optimization/ConjugateGradientMinimizer.cs b/src/Numerics/Optimization/ConjugateGradientMinimizer.cs new file mode 100644 index 00000000..f33a9043 --- /dev/null +++ b/src/Numerics/Optimization/ConjugateGradientMinimizer.cs @@ -0,0 +1,115 @@ +using System; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.Optimization.LineSearch; + +namespace MathNet.Numerics.Optimization +{ + public class ConjugateGradientMinimizer + { + public double GradientTolerance { get; set; } + public int MaximumIterations { get; set; } + + public ConjugateGradientMinimizer(double gradientTolerance, int maximumIterations) + { + GradientTolerance = gradientTolerance; + MaximumIterations = maximumIterations; + } + + public MinimizationResult FindMinimum(IObjectiveFunction objective, Vector initialGuess) + { + if (!objective.IsGradientSupported) + throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for ConjugateGradient minimization."); + + objective.EvaluateAt(initialGuess); + var gradient = objective.Gradient; + ValidateGradient(objective); + + // Check that we're not already done + if (ExitCriteriaSatisfied(initialGuess, gradient)) + return new MinimizationResult(objective, 0, MinimizationResult.ExitCondition.AbsoluteGradient); + + // Set up line search algorithm + var lineSearcher = new WeakWolfeLineSearch(1e-4, 0.1, 1e-4, 1000); + + // First step + var steepestDirection = -gradient; + var searchDirection = steepestDirection; + double initialStepSize = 100 * GradientTolerance / (gradient * gradient); + + LineSearchResult result; + try + { + result = lineSearcher.FindConformingStep(objective, searchDirection, initialStepSize); + } + catch (Exception e) + { + throw new InnerOptimizationException("Line search failed.", e); + } + + objective = result.FunctionInfoAtMinimum; + ValidateGradient(objective); + + double stepSize = result.FinalStep; + + // Subsequent steps + int iterations = 1; + int totalLineSearchSteps = result.Iterations; + int iterationsWithNontrivialLineSearch = result.Iterations > 0 ? 0 : 1; + int steepestDescentResets = 0; + while (!ExitCriteriaSatisfied(objective.Point, objective.Gradient) && iterations < MaximumIterations) + { + var previousSteepestDirection = steepestDirection; + steepestDirection = -objective.Gradient; + var searchDirectionAdjuster = Math.Max(0, steepestDirection*(steepestDirection - previousSteepestDirection)/(previousSteepestDirection*previousSteepestDirection)); + searchDirection = steepestDirection + searchDirectionAdjuster * searchDirection; + if (searchDirection * objective.Gradient >= 0) + { + searchDirection = steepestDirection; + steepestDescentResets += 1; + } + + try + { + result = lineSearcher.FindConformingStep(objective, searchDirection, stepSize); + } + catch (Exception e) + { + throw new InnerOptimizationException("Line search failed.", e); + } + + iterationsWithNontrivialLineSearch += result.Iterations == 0 ? 1 : 0; + totalLineSearchSteps += result.Iterations; + stepSize = result.FinalStep; + objective = result.FunctionInfoAtMinimum; + iterations += 1; + } + + if (iterations == MaximumIterations) + { + throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations)); + } + + return new MinimizationWithLineSearchResult(objective, iterations, MinimizationResult.ExitCondition.AbsoluteGradient, totalLineSearchSteps, iterationsWithNontrivialLineSearch); + } + + bool ExitCriteriaSatisfied(Vector candidatePoint, Vector gradient) + { + return gradient.Norm(2.0) < GradientTolerance; + } + + void ValidateGradient(IObjectiveFunction objective) + { + foreach (var x in objective.Gradient) + { + if (Double.IsNaN(x) || Double.IsInfinity(x)) + throw new EvaluationException("Non-finite gradient returned.", objective); + } + } + + void ValidateObjective(IObjectiveFunction objective) + { + if (Double.IsNaN(objective.Value) || Double.IsInfinity(objective.Value)) + throw new EvaluationException("Non-finite objective function returned.", objective); + } + } +} diff --git a/src/Numerics/Optimization/Exceptions.cs b/src/Numerics/Optimization/Exceptions.cs new file mode 100644 index 00000000..4b95d8fa --- /dev/null +++ b/src/Numerics/Optimization/Exceptions.cs @@ -0,0 +1,52 @@ +using System; + +namespace MathNet.Numerics.Optimization +{ + public class OptimizationException : Exception + { + public OptimizationException(string message) + : base(message) { } + + public OptimizationException(string message, Exception innerException) + : base(message, innerException) { } + } + + public class MaximumIterationsException : OptimizationException + { + public MaximumIterationsException(string message) + : base(message) { } + } + + public class EvaluationException : OptimizationException + { + public IObjectiveFunction ObjectiveFunction { get; private set; } + + public EvaluationException(string message, IObjectiveFunction eval) + : base(message) + { + ObjectiveFunction = eval; + } + + public EvaluationException(string message, IObjectiveFunction eval, Exception innerException) + : base(message, innerException) + { + ObjectiveFunction = eval; + } + + } + + public class InnerOptimizationException : OptimizationException + { + public InnerOptimizationException(string message) + : base(message) { } + + public InnerOptimizationException(string message, Exception inner_exception) + : base(message, inner_exception) { } + } + + public class IncompatibleObjectiveException : OptimizationException + { + public IncompatibleObjectiveException(string message) + : base(message) { } + } +} diff --git a/src/Numerics/Optimization/GoldenSectionMinimizer.cs b/src/Numerics/Optimization/GoldenSectionMinimizer.cs new file mode 100644 index 00000000..80b09778 --- /dev/null +++ b/src/Numerics/Optimization/GoldenSectionMinimizer.cs @@ -0,0 +1,107 @@ +using System; + +namespace MathNet.Numerics.Optimization +{ + public class GoldenSectionMinimizer + { + public double XTolerance { get; set; } + public int MaximumIterations { get; set; } + public int MaximumExpansionSteps { get; set; } + public double LowerExpansionFactor { get; set; } + public double UpperExpansionFactor { get; set; } + + public GoldenSectionMinimizer(double xTolerance = 1e-5, int maxIterations = 1000, int maxExpansionSteps = 10, double lowerExpansionFactor = 2.0, double upperExpansionFactor = 2.0) + { + XTolerance = xTolerance; + MaximumIterations = maxIterations; + MaximumExpansionSteps = maxExpansionSteps; + LowerExpansionFactor = lowerExpansionFactor; + UpperExpansionFactor = upperExpansionFactor; + } + + public MinimizationResult1D FindMinimum(IObjectiveFunction1D objective, double lowerBound, double upperBound) + { + if (upperBound <= lowerBound) + throw new OptimizationException("Lower bound must be lower than upper bound."); + + double middlePointX = lowerBound + (upperBound - lowerBound)/(1 + Constants.GoldenRatio); + IEvaluation1D lower = objective.Evaluate(lowerBound); + IEvaluation1D middle = objective.Evaluate(middlePointX); + IEvaluation1D upper = objective.Evaluate(upperBound); + + ValueChecker(lower.Value, lowerBound); + ValueChecker(middle.Value, middlePointX); + ValueChecker(upper.Value, upperBound); + + int expansion_steps = 0; + while ((expansion_steps < this.MaximumExpansionSteps) && (upper.Value < middle.Value || lower.Value < middle.Value)) + { + if (lower.Value < middle.Value) + { + lowerBound = 0.5*(upperBound + lowerBound) - this.LowerExpansionFactor*0.5*(upperBound - lowerBound); + lower = objective.Evaluate(lowerBound); + } + + if (upper.Value < middle.Value) + { + upperBound = 0.5*(upperBound + lowerBound) + this.UpperExpansionFactor*0.5*(upperBound - lowerBound); + upper = objective.Evaluate(upperBound); + } + + middlePointX = lowerBound + (upperBound - lowerBound)/(1 + Constants.GoldenRatio); + middle = objective.Evaluate(middlePointX); + + expansion_steps += 1; + } + + if (upper.Value < middle.Value || lower.Value < middle.Value) + throw new OptimizationException("Lower and upper bounds do not necessarily bound a minimum."); + + int iterations = 0; + while (Math.Abs(upper.Point - lower.Point) > XTolerance && iterations < MaximumIterations) + { + double testX = lower.Point + (upper.Point - middle.Point); + var test = objective.Evaluate(testX); + ValueChecker(test.Value, testX); + + if (test.Point < middle.Point) + { + if (test.Value > middle.Value) + { + lower = test; + } + else + { + upper = middle; + middle = test; + } + } + else + { + if (test.Value > middle.Value) + { + upper = test; + } + else + { + lower = middle; + middle = test; + } + } + + iterations += 1; + } + + if (iterations == MaximumIterations) + throw new MaximumIterationsException("Max iterations reached."); + + return new MinimizationResult1D(middle, iterations, MinimizationResult.ExitCondition.BoundTolerance); + } + + void ValueChecker(double value, double point) + { + if (Double.IsNaN(value) || Double.IsInfinity(value)) + throw new Exception("Objective function returned non-finite value."); + } + } +} diff --git a/src/Numerics/Optimization/IObjectiveFunction.cs b/src/Numerics/Optimization/IObjectiveFunction.cs new file mode 100644 index 00000000..7c54093c --- /dev/null +++ b/src/Numerics/Optimization/IObjectiveFunction.cs @@ -0,0 +1,33 @@ +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization +{ + /// + /// Objective function with a frozen evaluation that must not be changed from the outside. + /// + public interface IObjectiveFunctionEvaluation + { + /// Create a new unevaluated and independent copy of this objective function + IObjectiveFunction CreateNew(); + + /// Create a new independent copy of this objective function, evaluated at the same point. + IObjectiveFunction Fork(); + + Vector Point { get; } + double Value { get; } + + bool IsGradientSupported { get; } + Vector Gradient { get; } + + bool IsHessianSupported { get; } + Matrix Hessian { get; } + } + + /// + /// Objective function with a mutable evaluation. + /// + public interface IObjectiveFunction : IObjectiveFunctionEvaluation + { + void EvaluateAt(Vector point); + } +} diff --git a/src/Numerics/Optimization/IUnconstrainedMinimizer.cs b/src/Numerics/Optimization/IUnconstrainedMinimizer.cs new file mode 100644 index 00000000..3c88e8ca --- /dev/null +++ b/src/Numerics/Optimization/IUnconstrainedMinimizer.cs @@ -0,0 +1,9 @@ +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization +{ + public interface IUnconstrainedMinimizer + { + MinimizationResult FindMinimum(IObjectiveFunction objective, Vector initialGuess); + } +} diff --git a/src/Numerics/Optimization/LineSearch/LineSearchResult.cs b/src/Numerics/Optimization/LineSearch/LineSearchResult.cs new file mode 100644 index 00000000..2473e891 --- /dev/null +++ b/src/Numerics/Optimization/LineSearch/LineSearchResult.cs @@ -0,0 +1,43 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2016 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.Optimization.LineSearch +{ + public class LineSearchResult : MinimizationResult + { + public double FinalStep { get; private set; } + + public LineSearchResult(IObjectiveFunction functionInfo, int iterations, double finalStep, ExitCondition reasonForExit) + : base(functionInfo, iterations, reasonForExit) + { + FinalStep = finalStep; + } + } +} diff --git a/src/Numerics/Optimization/LineSearch/StrongWolfeLineSearch.cs b/src/Numerics/Optimization/LineSearch/StrongWolfeLineSearch.cs new file mode 100644 index 00000000..69a8675a --- /dev/null +++ b/src/Numerics/Optimization/LineSearch/StrongWolfeLineSearch.cs @@ -0,0 +1,50 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2016 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; + +namespace MathNet.Numerics.Optimization.LineSearch +{ + public class StrongWolfeLineSearch : WolfeLineSearch + { + public StrongWolfeLineSearch(double c1, double c2, double parameterTolerance, int maxIterations = 10) + : base(c1, c2, parameterTolerance, maxIterations) + { + // Argument validation in base class + } + + protected override MinimizationResult.ExitCondition WolfeExitCondition { get { return MinimizationResult.ExitCondition.StrongWolfeCriteria; } } + + protected override bool WolfeCondition(double stepDd, double initialDd) + { + return Math.Abs(stepDd) > C2 * Math.Abs(initialDd); + } + } +} diff --git a/src/Numerics/Optimization/LineSearch/WeakWolfeLineSearch.cs b/src/Numerics/Optimization/LineSearch/WeakWolfeLineSearch.cs new file mode 100644 index 00000000..f3a3baf4 --- /dev/null +++ b/src/Numerics/Optimization/LineSearch/WeakWolfeLineSearch.cs @@ -0,0 +1,97 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2016 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization.LineSearch +{ + /// + /// Search for a step size alpha that satisfies the weak wolfe conditions. The weak Wolfe + /// Conditions are + /// i) Armijo Rule: f(x_k + alpha_k p_k) <= f(x_k) + c1 alpha_k p_k^T g(x_k) + /// ii) Curvature Condition: p_k^T g(x_k + alpha_k p_k) >= c2 p_k^T g(x_k) + /// where g(x) is the gradient of f(x), 0 < c1 < c2 < 1. + /// + /// Implementation is based on http://www.math.washington.edu/~burke/crs/408/lectures/L9-weak-Wolfe.pdf + /// + /// references: + /// http://en.wikipedia.org/wiki/Wolfe_conditions + /// http://www.math.washington.edu/~burke/crs/408/lectures/L9-weak-Wolfe.pdf + /// + public class WeakWolfeLineSearch : WolfeLineSearch + { + public WeakWolfeLineSearch(double c1, double c2, double parameterTolerance, int maxIterations = 10) + : base(c1,c2,parameterTolerance,maxIterations) + { + // Validation in base class + } + + protected override MinimizationResult.ExitCondition WolfeExitCondition + { + get { return MinimizationResult.ExitCondition.WeakWolfeCriteria; } + } + + protected override bool WolfeCondition(double stepDd, double initialDd) + { + return stepDd < C2 * initialDd; + } + + protected override void ValidateValue(IObjectiveFunction eval) + { + if (!IsFinite(eval.Value)) + { + throw new EvaluationException(String.Format("Non-finite value returned by objective function: {0}", eval.Value), eval); + } + } + + protected override void ValidateInputArguments(IObjectiveFunctionEvaluation startingPoint, Vector searchDirection, double initialStep, double upperBound) + { + if (!startingPoint.IsGradientSupported) + throw new ArgumentException("objective function does not support gradient"); + } + + protected override void ValidateGradient(IObjectiveFunction eval) + { + foreach (double x in eval.Gradient) + { + if (!IsFinite(x)) + { + throw new EvaluationException(string.Format("Non-finite value returned by gradient: {0}", x), eval); + } + } + } + + static bool IsFinite(double x) + { + return !(double.IsNaN(x) || double.IsInfinity(x)); + } + } +} diff --git a/src/Numerics/Optimization/LineSearch/WolfeLineSearch.cs b/src/Numerics/Optimization/LineSearch/WolfeLineSearch.cs new file mode 100644 index 00000000..73ad86cb --- /dev/null +++ b/src/Numerics/Optimization/LineSearch/WolfeLineSearch.cs @@ -0,0 +1,158 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2016 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using MathNet.Numerics.LinearAlgebra; +using System; + +namespace MathNet.Numerics.Optimization.LineSearch +{ + public abstract class WolfeLineSearch + { + protected double C1 { get; } + protected double C2 { get; } + protected double ParameterTolerance { get; } + protected int MaximumIterations { get; } + + public WolfeLineSearch(double c1, double c2, double parameterTolerance, int maxIterations = 10) + { + if (c1 <= 0) + throw new ArgumentException(string.Format("c1 {0} should be greater than 0", c1)); + if (c2 <= c1) + throw new ArgumentException(string.Format("c1 {0} should be less than c2 {1}", c1, c2)); + if (c2 >= 1) + throw new ArgumentException(string.Format("c2 {0} should be less than 1", c2)); + + C1 = c1; + C2 = c2; + ParameterTolerance = parameterTolerance; + MaximumIterations = maxIterations; + } + + /// Implemented following http://www.math.washington.edu/~burke/crs/408/lectures/L9-weak-Wolfe.pdf + /// The objective function being optimized, evaluated at the starting point of the search + /// Search direction + /// Initial size of the step in the search direction + public LineSearchResult FindConformingStep(IObjectiveFunctionEvaluation startingPoint, Vector searchDirection, double initialStep) + { + return FindConformingStep(startingPoint, searchDirection, initialStep, double.PositiveInfinity); + } + + /// + /// The objective function being optimized, evaluated at the starting point of the search + /// Search direction + /// Initial size of the step in the search direction + /// The upper bound + public LineSearchResult FindConformingStep(IObjectiveFunctionEvaluation startingPoint, Vector searchDirection, double initialStep, double upperBound) + { + ValidateInputArguments(startingPoint, searchDirection, initialStep, upperBound); + + double lowerBound = 0.0; + double step = initialStep; + + double initialValue = startingPoint.Value; + Vector initialGradient = startingPoint.Gradient; + + double initialDd = searchDirection * initialGradient; + + IObjectiveFunction objective = startingPoint.CreateNew(); + int ii; + MinimizationResult.ExitCondition reasonForExit = MinimizationResult.ExitCondition.None; + for (ii = 0; ii < MaximumIterations; ++ii) + { + objective.EvaluateAt(startingPoint.Point + searchDirection * step); + ValidateGradient(objective); + ValidateValue(objective); + + double stepDd = searchDirection * objective.Gradient; + + if (objective.Value > initialValue + C1 * step * initialDd) + { + upperBound = step; + step = 0.5 * (lowerBound + upperBound); + } + else if (WolfeCondition(stepDd,initialDd)) + { + lowerBound = step; + step = double.IsPositiveInfinity(upperBound) ? 2 * lowerBound : 0.5 * (lowerBound + upperBound); + } + else + { + reasonForExit = WolfeExitCondition; + break; + } + + if (!double.IsInfinity(upperBound)) + { + double maxRelChange = 0.0; + for (int jj = 0; jj < objective.Point.Count; ++jj) + { + double tmp = Math.Abs(searchDirection[jj] * (upperBound - lowerBound)) / Math.Max(Math.Abs(objective.Point[jj]), 1.0); + maxRelChange = Math.Max(maxRelChange, tmp); + } + if (maxRelChange < ParameterTolerance) + { + reasonForExit = MinimizationResult.ExitCondition.LackOfProgress; + break; + } + } + } + + if (ii == MaximumIterations && Double.IsPositiveInfinity(upperBound)) + { + throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached. Function appears to be unbounded in search direction.", MaximumIterations)); + } + + if (ii == MaximumIterations) + { + throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations)); + } + + return new LineSearchResult(objective, ii, step, reasonForExit); + } + protected abstract MinimizationResult.ExitCondition WolfeExitCondition { get; } + + protected abstract bool WolfeCondition(double stepDd, double initialDd); + + protected virtual void ValidateGradient(IObjectiveFunction objective) + { + } + protected virtual void ValidateValue(IObjectiveFunction objective) + { + } + + protected virtual void ValidateInputArguments(IObjectiveFunctionEvaluation startingPoint, Vector searchDirection, double initialStep, double upperBound) + { + + } + } + + + +} diff --git a/src/Numerics/Optimization/MinimizationResult.cs b/src/Numerics/Optimization/MinimizationResult.cs new file mode 100644 index 00000000..2a9d8fcd --- /dev/null +++ b/src/Numerics/Optimization/MinimizationResult.cs @@ -0,0 +1,32 @@ +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization +{ + public class MinimizationResult + { + public enum ExitCondition + { + None, + RelativeGradient, + LackOfProgress, + AbsoluteGradient, + WeakWolfeCriteria, + BoundTolerance, + StrongWolfeCriteria, + LackOfFunctionImprovement, + Converged + } + + public Vector MinimizingPoint { get { return FunctionInfoAtMinimum.Point; } } + public IObjectiveFunction FunctionInfoAtMinimum { get; private set; } + public int Iterations { get; private set; } + public ExitCondition ReasonForExit { get; private set; } + + public MinimizationResult(IObjectiveFunction functionInfo, int iterations, ExitCondition reasonForExit) + { + FunctionInfoAtMinimum = functionInfo; + Iterations = iterations; + ReasonForExit = reasonForExit; + } + } +} diff --git a/src/Numerics/Optimization/MinimizationResult1D.cs b/src/Numerics/Optimization/MinimizationResult1D.cs new file mode 100644 index 00000000..c1e6f65a --- /dev/null +++ b/src/Numerics/Optimization/MinimizationResult1D.cs @@ -0,0 +1,17 @@ +namespace MathNet.Numerics.Optimization +{ + public class MinimizationResult1D + { + public double MinimizingPoint { get { return FunctionInfoAtMinimum.Point; } } + public IEvaluation1D FunctionInfoAtMinimum { get; private set; } + public int Iterations { get; private set; } + public MinimizationResult.ExitCondition ReasonForExit { get; private set; } + + public MinimizationResult1D(IEvaluation1D functionInfo, int iterations, MinimizationResult.ExitCondition reasonForExit) + { + FunctionInfoAtMinimum = functionInfo; + Iterations = iterations; + ReasonForExit = reasonForExit; + } + } +} diff --git a/src/Numerics/Optimization/MinimizationWithLineSearchResult.cs b/src/Numerics/Optimization/MinimizationWithLineSearchResult.cs new file mode 100644 index 00000000..72002bee --- /dev/null +++ b/src/Numerics/Optimization/MinimizationWithLineSearchResult.cs @@ -0,0 +1,15 @@ +namespace MathNet.Numerics.Optimization +{ + public class MinimizationWithLineSearchResult : MinimizationResult + { + public int TotalLineSearchIterations { get; private set; } + public int IterationsWithNonTrivialLineSearch { get; private set; } + + public MinimizationWithLineSearchResult(IObjectiveFunction functionInfo, int iterations, ExitCondition reasonForExit, int totalLineSearchIterations, int iterationsWithNonTrivialLineSearch) + : base(functionInfo, iterations, reasonForExit) + { + TotalLineSearchIterations = totalLineSearchIterations; + IterationsWithNonTrivialLineSearch = iterationsWithNonTrivialLineSearch; + } + } +} diff --git a/src/Numerics/Optimization/NelderMeadSimplex.cs b/src/Numerics/Optimization/NelderMeadSimplex.cs new file mode 100644 index 00000000..0024e96f --- /dev/null +++ b/src/Numerics/Optimization/NelderMeadSimplex.cs @@ -0,0 +1,415 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2015 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +// Converted from code relased with a MIT liscense available at https://code.google.com/p/nelder-mead-simplex/ + +using MathNet.Numerics.LinearAlgebra; +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace MathNet.Numerics.Optimization +{ + /// + /// Class implementing the Nelder-Mead simplex algorithm, used to find a minima when no gradient is available. + /// Called fminsearch() in Matlab. A description of the algorithm can be found at + /// http://se.mathworks.com/help/matlab/math/optimizing-nonlinear-functions.html#bsgpq6p-11 + /// or + /// https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method + /// + public sealed class NelderMeadSimplex + { + private static readonly double JITTER = 1e-10d; // a small value used to protect against floating point noise + + public double ConvergenceTolerance { get; set; } + public int MaximumIterations { get; set; } + + public NelderMeadSimplex(double convergenceTolerance, int maximumIterations) + { + ConvergenceTolerance = convergenceTolerance; + MaximumIterations = maximumIterations; + } + + /// + /// Finds the minimum of the objective function without an intial pertubation, the default values used + /// by fminsearch() in Matlab are used instead + /// http://se.mathworks.com/help/matlab/math/optimizing-nonlinear-functions.html#bsgpq6p-11 + /// + /// The objective function, no gradient or hessian needed + /// The intial guess + /// The minimum point + public MinimizationResult FindMinimum(IObjectiveFunction objectiveFunction, Vector initialGuess) + { + var initalPertubation = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(initialGuess.Count); + for (int i = 0; i < initialGuess.Count; i++) + { + initalPertubation[i] = initialGuess[i] == 0.0 ? 0.00025 : initialGuess[i] * 0.05; + } + return FindMinimum(objectiveFunction, initialGuess, initalPertubation); + } + + /// + /// Finds the minimum of the objective function with an intial pertubation + /// + /// The objective function, no gradient or hessian needed + /// The intial guess + /// The inital pertubation + /// The minimum point + public MinimizationResult FindMinimum(IObjectiveFunction objectiveFunction, Vector initialGuess, Vector initalPertubation) + { + // confirm that we are in a position to commence + if (objectiveFunction == null) + throw new ArgumentNullException("objectiveFunction","ObjectiveFunction must be set to a valid ObjectiveFunctionDelegate"); + + if (initialGuess == null) + throw new ArgumentNullException("initialGuess", "initialGuess must be initialized"); + + if (initialGuess == null) + throw new ArgumentNullException("initalPertubation", "initalPertubation must be initialized, if unknown use overloaded version of FindMinimum()"); + + SimplexConstant[] simplexConstants = SimplexConstant.CreateSimplexConstantsFromVectors(initialGuess,initalPertubation); + + // create the initial simplex + int numDimensions = simplexConstants.Length; + int numVertices = numDimensions + 1; + Vector[] 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 >= MaximumIterations) + { + throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations)); + } + } + var regressionResult = new MinimizationResult(objectiveFunction, evaluationCount, exitCondition); + return regressionResult; + } + + /// + /// Evaluate the objective function at each vertex to create a corresponding + /// list of error values for each vertex + /// + /// + /// + /// + private static double[] InitializeErrorValues(Vector[] 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; + } + + /// + /// Check whether the points in the error profile have so little range that we + /// consider ourselves to have converged + /// + /// + /// + /// + /// + 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; + } + } + + /// + /// Examine all error values to determine the ErrorProfile + /// + /// + /// + 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; + } + + /// + /// Construct an initial simplex, given starting guesses for the constants, and + /// initial step sizes for each dimension + /// + /// + /// + private static Vector[] InitializeVertices(SimplexConstant[] simplexConstants) + { + int numDimensions = simplexConstants.Length; + Vector[] vertices = new Vector[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 unitVector = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(numDimensions); + unitVector[i] = 1; + vertices[i + 1] = p0.Add(unitVector.Multiply(scale)); + } + return vertices; + } + + /// + /// Test a scaling operation of the high point, and replace it if it is an improvement + /// + /// + /// + /// + /// + /// + /// + private static double TryToScaleSimplex(double scaleFactor, ref ErrorProfile errorProfile, Vector[] vertices, + double[] errorValues, IObjectiveFunction objectiveFunction) + { + // find the centroid through which we will reflect + Vector centroid = ComputeCentroid(vertices, errorProfile); + + // define the vector from the centroid to the high point + Vector centroidToHighPoint = vertices[errorProfile.HighestIndex].Subtract(centroid); + + // scale and position the vector to determine the new trial point + Vector 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; + } + + /// + /// Contract the simplex uniformly around the lowest point + /// + /// + /// + /// + /// + private static void ShrinkSimplex(ErrorProfile errorProfile, Vector[] vertices, double[] errorValues, + IObjectiveFunction objectiveFunction) + { + Vector 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; + } + } + } + + /// + /// Compute the centroid of all points except the worst + /// + /// + /// + /// + private static Vector ComputeCentroid(Vector[] vertices, ErrorProfile errorProfile) + { + int numVertices = vertices.Length; + // find the centroid of all points except the worst one + Vector 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; + } + + /// + /// The value of the constant + /// + 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[] CreateSimplexConstantsFromVectors(Vector initialGuess, Vector initialPertubation) + { + var constants = new SimplexConstant[initialGuess.Count]; + for (int i = 0; i < constants.Length;i++ ) + { + constants[i] = new SimplexConstant(initialGuess[i], initialPertubation[i]); + } + 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; } + } + } + } + + + +} diff --git a/src/Numerics/Optimization/NewtonMinimizer.cs b/src/Numerics/Optimization/NewtonMinimizer.cs new file mode 100644 index 00000000..2a0419bd --- /dev/null +++ b/src/Numerics/Optimization/NewtonMinimizer.cs @@ -0,0 +1,130 @@ +using System; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.Optimization.LineSearch; + +namespace MathNet.Numerics.Optimization +{ + public class NewtonMinimizer + { + public double GradientTolerance { get; set; } + public int MaximumIterations { get; set; } + public bool UseLineSearch { get; set; } + + public NewtonMinimizer(double gradientTolerance, int maximumIterations, bool useLineSearch = false) + { + GradientTolerance = gradientTolerance; + MaximumIterations = maximumIterations; + UseLineSearch = useLineSearch; + } + + public MinimizationResult FindMinimum(IObjectiveFunction objective, Vector initialGuess) + { + if (!objective.IsGradientSupported) + { + throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for Newton minimization."); + } + + if (!objective.IsHessianSupported) + { + throw new IncompatibleObjectiveException("Hessian not supported in objective function, but required for Newton minimization."); + } + + // Check that we're not already done + objective.EvaluateAt(initialGuess); + ValidateGradient(objective); + if (ExitCriteriaSatisfied(objective.Gradient)) + { + return new MinimizationResult(objective, 0, MinimizationResult.ExitCondition.AbsoluteGradient); + } + + // Set up line search algorithm + var lineSearcher = new WeakWolfeLineSearch(1e-4, 0.9, 1e-4, maxIterations: 1000); + + // Subsequent steps + int iterations = 0; + int totalLineSearchSteps = 0; + int iterationsWithNontrivialLineSearch = 0; + bool tmpLineSearch = false; + while (!ExitCriteriaSatisfied(objective.Gradient) && iterations < MaximumIterations) + { + ValidateHessian(objective); + + var searchDirection = objective.Hessian.LU().Solve(-objective.Gradient); + if (searchDirection * objective.Gradient >= 0) + { + searchDirection = -objective.Gradient; + tmpLineSearch = true; + } + + if (UseLineSearch || tmpLineSearch) + { + LineSearchResult result; + try + { + result = lineSearcher.FindConformingStep(objective, searchDirection, 1.0); + } + catch (Exception e) + { + throw new InnerOptimizationException("Line search failed.", e); + } + + iterationsWithNontrivialLineSearch += result.Iterations > 0 ? 1 : 0; + totalLineSearchSteps += result.Iterations; + objective = result.FunctionInfoAtMinimum; + } + else + { + objective.EvaluateAt(objective.Point + searchDirection); + } + + ValidateGradient(objective); + + tmpLineSearch = false; + iterations += 1; + } + + if (iterations == MaximumIterations) + { + throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations)); + } + + return new MinimizationWithLineSearchResult(objective, iterations, MinimizationResult.ExitCondition.AbsoluteGradient, totalLineSearchSteps, iterationsWithNontrivialLineSearch); + } + + bool ExitCriteriaSatisfied(Vector gradient) + { + return gradient.Norm(2.0) < GradientTolerance; + } + + static void ValidateGradient(IObjectiveFunction eval) + { + foreach (var x in eval.Gradient) + { + if (Double.IsNaN(x) || Double.IsInfinity(x)) + { + throw new EvaluationException("Non-finite gradient returned.", eval); + } + } + } + + private void ValidateObjective(IObjectiveFunction eval) + { + if (Double.IsNaN(eval.Value) || Double.IsInfinity(eval.Value)) + throw new EvaluationException("Non-finite objective function returned.", eval); + } + + private void ValidateHessian(IObjectiveFunction eval) + { + for (int ii = 0; ii < eval.Hessian.RowCount; ++ii) + { + for (int jj = 0; jj < eval.Hessian.ColumnCount; ++jj) + { + if (Double.IsNaN(eval.Hessian[ii, jj]) || Double.IsInfinity(eval.Hessian[ii, jj])) + { + throw new EvaluationException("Non-finite Hessian returned.", eval); + } + } + } + } + } +} diff --git a/src/Numerics/Optimization/ObjectiveFunction.cs b/src/Numerics/Optimization/ObjectiveFunction.cs new file mode 100644 index 00000000..e3a1d577 --- /dev/null +++ b/src/Numerics/Optimization/ObjectiveFunction.cs @@ -0,0 +1,65 @@ +using System; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.Optimization.ObjectiveFunctions; + +namespace MathNet.Numerics.Optimization +{ + public static class ObjectiveFunction + { + /// + /// Objective function where neither Gradient nor Hessian is available. + /// + public static IObjectiveFunction Value(Func, double> function) + { + return new ValueObjectiveFunction(function); + } + + /// + /// Objective function where the Gradient is available. Greedy evaluation. + /// + public static IObjectiveFunction Gradient(Func, Tuple>> function) + { + return new GradientObjectiveFunction(function); + } + + /// + /// Objective function where the Gradient is available. Lazy evaluation. + /// + public static IObjectiveFunction Gradient(Func, double> function, Func, Vector> gradient) + { + return new LazyObjectiveFunction(function, gradient: gradient); + } + + /// + /// Objective function where the Hessian is available. Greedy evaluation. + /// + public static IObjectiveFunction Hessian(Func, Tuple>> function) + { + return new HessianObjectiveFunction(function); + } + + /// + /// Objective function where the Hessian is available. Lazy evaluation. + /// + public static IObjectiveFunction Hessian(Func, double> function, Func, Matrix> hessian) + { + return new LazyObjectiveFunction(function, hessian: hessian); + } + + /// + /// Objective function where both Gradient and Hessian are available. Greedy evaluation. + /// + public static IObjectiveFunction GradientHessian(Func, Tuple, Matrix>> function) + { + return new GradientHessianObjectiveFunction(function); + } + + /// + /// Objective function where both Gradient and Hessian are available. Lazy evaluation. + /// + public static IObjectiveFunction GradientHessian(Func, double> function, Func, Vector> gradient, Func, Matrix> hessian) + { + return new LazyObjectiveFunction(function, gradient: gradient, hessian: hessian); + } + } +} diff --git a/src/Numerics/Optimization/ObjectiveFunction1D.cs b/src/Numerics/Optimization/ObjectiveFunction1D.cs new file mode 100644 index 00000000..00a35d65 --- /dev/null +++ b/src/Numerics/Optimization/ObjectiveFunction1D.cs @@ -0,0 +1,98 @@ +using System; + +namespace MathNet.Numerics.Optimization +{ + public interface IEvaluation1D + { + double Point { get; } + double Value { get; } + double Derivative { get; } + double SecondDerivative { get; } + } + + public interface IObjectiveFunction1D + { + bool DerivativeSupported { get; } + bool SecondDerivativeSupported { get; } + IEvaluation1D Evaluate(double point); + } + + public class CachedEvaluation1D : IEvaluation1D + { + private double? _value; + private double? _derivative; + private double? _secondDerivative; + private readonly SimpleObjectiveFunction1D _objectiveObject; + private readonly double _point; + + public CachedEvaluation1D(SimpleObjectiveFunction1D f, double point) + { + _objectiveObject = f; + _point = point; + } + private double SetValue() + { + _value = _objectiveObject.Objective(_point); + return _value.Value; + } + private double SetDerivative() + { + _derivative = _objectiveObject.Derivative(_point); + return _derivative.Value; + } + private double SetSecondDerivative() + { + _secondDerivative = _objectiveObject.SecondDerivative(_point); + return _secondDerivative.Value; + } + + public double Point { get { return _point; } } + public double Value { get { return _value ?? SetValue(); } } + public double Derivative { get { return _derivative ?? SetDerivative(); } } + public double SecondDerivative { get { return _secondDerivative ?? SetSecondDerivative(); } } + + } + + public class SimpleObjectiveFunction1D : IObjectiveFunction1D + { + public Func Objective { get; private set; } + public Func Derivative { get; private set; } + public Func SecondDerivative { get; private set; } + + public SimpleObjectiveFunction1D(Func objective) + { + Objective = objective; + Derivative = null; + SecondDerivative = null; + } + + public SimpleObjectiveFunction1D(Func objective, Func derivative) + { + Objective = objective; + Derivative = derivative; + SecondDerivative = null; + } + + public SimpleObjectiveFunction1D(Func objective, Func derivative, Func secondDerivative) + { + Objective = objective; + Derivative = derivative; + SecondDerivative = secondDerivative; + } + + public bool DerivativeSupported + { + get { return Derivative != null; } + } + + public bool SecondDerivativeSupported + { + get { return SecondDerivative != null; } + } + + public IEvaluation1D Evaluate(double point) + { + return new CachedEvaluation1D(this, point); + } + } +} diff --git a/src/Numerics/Optimization/ObjectiveFunctions/ForwardDifferenceGradientObjectiveFunction.cs b/src/Numerics/Optimization/ObjectiveFunctions/ForwardDifferenceGradientObjectiveFunction.cs new file mode 100644 index 00000000..21b3d029 --- /dev/null +++ b/src/Numerics/Optimization/ObjectiveFunctions/ForwardDifferenceGradientObjectiveFunction.cs @@ -0,0 +1,145 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization.ObjectiveFunctions +{ + /// + /// Adapts an objective function with only value implemented + /// to provide a gradient as well. Gradient calculation is + /// done using the finite difference method, specifically + /// forward differences. + /// + /// For each gradient computed, the algorithm requires an + /// additional number of function evaluations equal to the + /// functions's number of input parameters. + /// + public class ForwardDifferenceGradientObjectiveFunction : IObjectiveFunction + { + public IObjectiveFunction InnerObjectiveFunction { get; protected set; } + protected Vector LowerBound { get; set; } + protected Vector UpperBound { get; set; } + + protected bool ValueEvaluated { get; set; } = false; + protected bool GradientEvaluated { get; set; } = false; + private Vector _gradient; + + public double MinimumIncrement { get; set; } + public double RelativeIncrement { get; set; } + + public ForwardDifferenceGradientObjectiveFunction(IObjectiveFunction valueOnlyObj, Vector lowerBound, Vector upperBound, double relativeIncrement=1e-5, double minimumIncrement=1e-8) + { + InnerObjectiveFunction = valueOnlyObj; + LowerBound = lowerBound; + UpperBound = upperBound; + _gradient = new LinearAlgebra.Double.DenseVector(LowerBound.Count); + RelativeIncrement = relativeIncrement; + MinimumIncrement = minimumIncrement; + } + + protected void EvaluateValue() + { + ValueEvaluated = true; + } + + protected void EvaluateGradient() + { + if (!ValueEvaluated) + EvaluateValue(); + + var tmp_point = Point.Clone(); + var tmp_obj = InnerObjectiveFunction.CreateNew(); + for (int ii = 0; ii < _gradient.Count; ++ii) + { + var orig_point = tmp_point[ii]; + var rel_incr = orig_point * RelativeIncrement; + var h = Math.Max(rel_incr, MinimumIncrement); + var mult = 1; + if (orig_point + h > UpperBound[ii]) + mult = -1; + + tmp_point[ii] = orig_point + mult*h; + tmp_obj.EvaluateAt(tmp_point); + double bumped_value = tmp_obj.Value; + _gradient[ii] = (mult * bumped_value - mult * InnerObjectiveFunction.Value) / h; + + tmp_point[ii] = orig_point; + } + GradientEvaluated = true; + } + + public Vector Gradient + { + get + { + if (!GradientEvaluated) + EvaluateGradient(); + return _gradient; + } + protected set { _gradient = value; } + } + + public Matrix Hessian + { + get + { + throw new NotImplementedException(); + } + } + + public bool IsGradientSupported + { + get + { + return true; + } + } + + public bool IsHessianSupported + { + get + { + return false; + } + } + + public Vector Point { get; protected set; } + + public double Value + { + get + { + if (!ValueEvaluated) + EvaluateValue(); + return this.InnerObjectiveFunction.Value; + } + } + + public IObjectiveFunction CreateNew() + { + var tmp = new ForwardDifferenceGradientObjectiveFunction(this.InnerObjectiveFunction.CreateNew(), LowerBound, UpperBound, this.RelativeIncrement, this.MinimumIncrement); + return tmp; + } + + public void EvaluateAt(Vector point) + { + Point = point; + ValueEvaluated = false; + GradientEvaluated = false; + InnerObjectiveFunction.EvaluateAt(point); + } + + public IObjectiveFunction Fork() + { + return new ForwardDifferenceGradientObjectiveFunction(this.InnerObjectiveFunction.Fork(), LowerBound, UpperBound, this.RelativeIncrement, this.MinimumIncrement) + { + Point = Point?.Clone(), + GradientEvaluated = GradientEvaluated, + ValueEvaluated = ValueEvaluated, + _gradient = _gradient?.Clone() + }; + } + } +} diff --git a/src/Numerics/Optimization/ObjectiveFunctions/GradientHessianObjectiveFunction.cs b/src/Numerics/Optimization/ObjectiveFunctions/GradientHessianObjectiveFunction.cs new file mode 100644 index 00000000..24d2a914 --- /dev/null +++ b/src/Numerics/Optimization/ObjectiveFunctions/GradientHessianObjectiveFunction.cs @@ -0,0 +1,57 @@ +using System; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization.ObjectiveFunctions +{ + internal class GradientHessianObjectiveFunction : IObjectiveFunction + { + readonly Func, Tuple, Matrix>> _function; + + public GradientHessianObjectiveFunction(Func, Tuple, Matrix>> function) + { + _function = function; + } + + public IObjectiveFunction CreateNew() + { + return new GradientHessianObjectiveFunction(_function); + } + + public IObjectiveFunction Fork() + { + // no need to deep-clone values since they are replaced on evaluation + return new GradientHessianObjectiveFunction(_function) + { + Point = Point, + Value = Value, + Gradient = Gradient, + Hessian = Hessian + }; + } + + public bool IsGradientSupported + { + get { return true; } + } + + public bool IsHessianSupported + { + get { return true; } + } + + public void EvaluateAt(Vector point) + { + Point = point; + + var result = _function(point); + Value = result.Item1; + Gradient = result.Item2; + Hessian = result.Item3; + } + + public Vector Point { get; private set; } + public double Value { get; private set; } + public Vector Gradient { get; private set; } + public Matrix Hessian { get; private set; } + } +} \ No newline at end of file diff --git a/src/Numerics/Optimization/ObjectiveFunctions/GradientObjectiveFunction.cs b/src/Numerics/Optimization/ObjectiveFunctions/GradientObjectiveFunction.cs new file mode 100644 index 00000000..8d967f01 --- /dev/null +++ b/src/Numerics/Optimization/ObjectiveFunctions/GradientObjectiveFunction.cs @@ -0,0 +1,59 @@ +using System; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization.ObjectiveFunctions +{ + internal class GradientObjectiveFunction : IObjectiveFunction + { + readonly Func, Tuple>> _function; + + public GradientObjectiveFunction(Func, Tuple>> function) + { + _function = function; + } + + public IObjectiveFunction CreateNew() + { + return new GradientObjectiveFunction(_function); + } + + public IObjectiveFunction Fork() + { + // no need to deep-clone values since they are replaced on evaluation + return new GradientObjectiveFunction(_function) + { + Point = Point, + Value = Value, + Gradient = Gradient + }; + } + + public bool IsGradientSupported + { + get { return true; } + } + + public bool IsHessianSupported + { + get { return false; } + } + + public void EvaluateAt(Vector point) + { + Point = point; + + var result = _function(point); + Value = result.Item1; + Gradient = result.Item2; + } + + public Vector Point { get; private set; } + public double Value { get; private set; } + public Vector Gradient { get; private set; } + + public Matrix Hessian + { + get { throw new NotSupportedException(); } + } + } +} \ No newline at end of file diff --git a/src/Numerics/Optimization/ObjectiveFunctions/HessianObjectiveFunction.cs b/src/Numerics/Optimization/ObjectiveFunctions/HessianObjectiveFunction.cs new file mode 100644 index 00000000..54215980 --- /dev/null +++ b/src/Numerics/Optimization/ObjectiveFunctions/HessianObjectiveFunction.cs @@ -0,0 +1,59 @@ +using System; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization.ObjectiveFunctions +{ + internal class HessianObjectiveFunction : IObjectiveFunction + { + readonly Func, Tuple>> _function; + + public HessianObjectiveFunction(Func, Tuple>> function) + { + _function = function; + } + + public IObjectiveFunction CreateNew() + { + return new HessianObjectiveFunction(_function); + } + + public IObjectiveFunction Fork() + { + // no need to deep-clone values since they are replaced on evaluation + return new HessianObjectiveFunction(_function) + { + Point = Point, + Value = Value, + Hessian = Hessian + }; + } + + public bool IsGradientSupported + { + get { return false; } + } + + public bool IsHessianSupported + { + get { return true; } + } + + public void EvaluateAt(Vector point) + { + Point = point; + + var result = _function(point); + Value = result.Item1; + Hessian = result.Item2; + } + + public Vector Point { get; private set; } + public double Value { get; private set; } + public Matrix Hessian { get; private set; } + + public Vector Gradient + { + get { throw new NotSupportedException(); } + } + } +} \ No newline at end of file diff --git a/src/Numerics/Optimization/ObjectiveFunctions/LazyObjectiveFunction.cs b/src/Numerics/Optimization/ObjectiveFunctions/LazyObjectiveFunction.cs new file mode 100644 index 00000000..2639e3a9 --- /dev/null +++ b/src/Numerics/Optimization/ObjectiveFunctions/LazyObjectiveFunction.cs @@ -0,0 +1,142 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2016 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization.ObjectiveFunctions +{ + internal class LazyObjectiveFunction : IObjectiveFunction + { + readonly Func, double> _function; + readonly Func, Vector> _gradient; + readonly Func, Matrix> _hessian; + + Vector _point; + + bool _hasFunctionValue; + double _functionValue; + + bool _hasGradientValue; + Vector _gradientValue; + + bool _hasHessianValue; + Matrix _hessianValue; + + public LazyObjectiveFunction(Func, double> function, Func, Vector> gradient = null, Func, Matrix> hessian = null) + { + _function = function; + _gradient = gradient; + _hessian = hessian; + + IsGradientSupported = gradient != null; + IsHessianSupported = hessian != null; + } + + public IObjectiveFunction CreateNew() + { + return new LazyObjectiveFunction(_function, _gradient, _hessian); + } + + public IObjectiveFunction Fork() + { + // no need to deep-clone values since they are replaced on evaluation + return new LazyObjectiveFunction(_function, _gradient, _hessian) + { + _point = _point, + _hasFunctionValue = _hasFunctionValue, + _functionValue = _functionValue, + _hasGradientValue = _hasGradientValue, + _gradientValue = _gradientValue, + _hasHessianValue = _hasHessianValue, + _hessianValue = _hessianValue + }; + } + + public bool IsGradientSupported { get; private set; } + public bool IsHessianSupported { get; private set; } + + public void EvaluateAt(Vector point) + { + _point = point; + _hasFunctionValue = false; + _hasGradientValue = false; + _hasHessianValue = false; + + // don't keep references unnecessarily + _gradientValue = null; + _hessianValue = null; + } + + public Vector Point + { + get { return _point; } + } + + public double Value + { + get + { + if (!_hasFunctionValue) + { + _functionValue = _function(_point); + _hasFunctionValue = true; + } + return _functionValue; + } + } + + public Vector Gradient + { + get + { + if (!_hasGradientValue) + { + _gradientValue = _gradient(_point); + _hasGradientValue = true; + } + return _gradientValue; + } + } + + public Matrix Hessian + { + get + { + if (!_hasHessianValue) + { + _hessianValue = _hessian(_point); + _hasHessianValue = true; + } + return _hessianValue; + } + } + } +} \ No newline at end of file diff --git a/src/Numerics/Optimization/ObjectiveFunctions/LazyObjectiveFunctionBase.cs b/src/Numerics/Optimization/ObjectiveFunctions/LazyObjectiveFunctionBase.cs new file mode 100644 index 00000000..d2ea4ad6 --- /dev/null +++ b/src/Numerics/Optimization/ObjectiveFunctions/LazyObjectiveFunctionBase.cs @@ -0,0 +1,149 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2016 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization.ObjectiveFunctions +{ + public abstract class LazyObjectiveFunctionBase : IObjectiveFunction + { + Vector _point; + + protected bool HasFunctionValue { get; set; } + protected double FunctionValue { get; set; } + + protected bool HasGradientValue { get; set; } + protected Vector GradientValue { get; set; } + + protected bool HasHessianValue { get; set; } + protected Matrix HessianValue { get; set; } + + protected LazyObjectiveFunctionBase(bool gradientSupported, bool hessianSupported) + { + IsGradientSupported = gradientSupported; + IsHessianSupported = hessianSupported; + } + + public abstract IObjectiveFunction CreateNew(); + + public virtual IObjectiveFunction Fork() + { + // we need to deep-clone values since they may be updated inplace on evaluation + LazyObjectiveFunctionBase fork = (LazyObjectiveFunctionBase)CreateNew(); + fork._point = _point?.Clone(); + fork.HasFunctionValue = HasFunctionValue; + fork.FunctionValue = FunctionValue; + fork.HasGradientValue = HasGradientValue; + fork.GradientValue = GradientValue?.Clone(); + fork.HasHessianValue = HasHessianValue; + fork.HessianValue = HessianValue?.Clone(); + return fork; + } + + public bool IsGradientSupported { get; private set; } + public bool IsHessianSupported { get; private set; } + + public void EvaluateAt(Vector point) + { + _point = point; + HasFunctionValue = false; + HasGradientValue = false; + HasHessianValue = false; + } + + protected abstract void EvaluateValue(); + + protected virtual void EvaluateGradient() + { + Gradient = null; + } + + protected virtual void EvaluateHessian() + { + Hessian = null; + } + + public Vector Point + { + get { return _point; } + } + + public double Value + { + get + { + if (!HasFunctionValue) + { + EvaluateValue(); + } + return FunctionValue; + } + protected set + { + FunctionValue = value; + HasFunctionValue = true; + } + } + + public Vector Gradient + { + get + { + if (!HasGradientValue) + { + EvaluateGradient(); + } + return GradientValue; + } + protected set + { + GradientValue = value; + HasGradientValue = true; + } + } + + public Matrix Hessian + { + get + { + if (!HasHessianValue) + { + EvaluateHessian(); + } + return HessianValue; + } + protected set + { + HessianValue = value; + HasHessianValue = true; + } + } + } +} diff --git a/src/Numerics/Optimization/ObjectiveFunctions/ObjectiveFunctionBase.cs b/src/Numerics/Optimization/ObjectiveFunctions/ObjectiveFunctionBase.cs new file mode 100644 index 00000000..f212bcc0 --- /dev/null +++ b/src/Numerics/Optimization/ObjectiveFunctions/ObjectiveFunctionBase.cs @@ -0,0 +1,42 @@ +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization.ObjectiveFunctions +{ + public abstract class ObjectiveFunctionBase : IObjectiveFunction + { + protected ObjectiveFunctionBase(bool isGradientSupported, bool isHessianSupported) + { + IsGradientSupported = isGradientSupported; + IsHessianSupported = isHessianSupported; + } + + public abstract IObjectiveFunction CreateNew(); + + public virtual IObjectiveFunction Fork() + { + // we need to deep-clone values since they may be updated inplace on evaluation + ObjectiveFunctionBase objective = (ObjectiveFunctionBase)CreateNew(); + objective.Point = Point == null ? null : Point.Clone(); + objective.Value = Value; + objective.Gradient = Gradient == null ? null : Gradient.Clone(); + objective.Hessian = Hessian == null ? null : Hessian.Clone(); + return objective; + } + + public bool IsGradientSupported { get; private set; } + public bool IsHessianSupported { get; private set; } + + public void EvaluateAt(Vector point) + { + Point = point; + Evaluate(); + } + + protected abstract void Evaluate(); + + public Vector Point { get; private set; } + public double Value { get; protected set; } + public Vector Gradient { get; protected set; } + public Matrix Hessian { get; protected set; } + } +} diff --git a/src/Numerics/Optimization/ObjectiveFunctions/ValueObjectiveFunction.cs b/src/Numerics/Optimization/ObjectiveFunctions/ValueObjectiveFunction.cs new file mode 100644 index 00000000..19e75d74 --- /dev/null +++ b/src/Numerics/Optimization/ObjectiveFunctions/ValueObjectiveFunction.cs @@ -0,0 +1,59 @@ +using System; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization.ObjectiveFunctions +{ + internal class ValueObjectiveFunction : IObjectiveFunction + { + readonly Func, double> _function; + + public ValueObjectiveFunction(Func, double> function) + { + _function = function; + } + + public IObjectiveFunction CreateNew() + { + return new ValueObjectiveFunction(_function); + } + + public IObjectiveFunction Fork() + { + // no need to deep-clone values since they are replaced on evaluation + return new ValueObjectiveFunction(_function) + { + Point = Point, + Value = Value, + }; + } + + public bool IsGradientSupported + { + get { return false; } + } + + public bool IsHessianSupported + { + get { return false; } + } + + public void EvaluateAt(Vector point) + { + Point = point; + Value = _function(point); + } + + public Vector Point { get; private set; } + public double Value { get; private set; } + + public Matrix Hessian + { + get { throw new NotSupportedException(); } + } + + public Vector Gradient + { + get { throw new NotSupportedException(); } + } + } +} \ No newline at end of file diff --git a/src/Numerics/Optimization/OptimizationResult.cs b/src/Numerics/Optimization/OptimizationResult.cs new file mode 100644 index 00000000..512e11d9 --- /dev/null +++ b/src/Numerics/Optimization/OptimizationResult.cs @@ -0,0 +1,8 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace MathNet.Numerics.Optimization +{ +} diff --git a/src/Numerics/Optimization/QuadraticGradientProjectionSearch.cs b/src/Numerics/Optimization/QuadraticGradientProjectionSearch.cs new file mode 100644 index 00000000..abf97f1d --- /dev/null +++ b/src/Numerics/Optimization/QuadraticGradientProjectionSearch.cs @@ -0,0 +1,99 @@ +using System; +using System.Collections.Generic; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization +{ + public static class QuadraticGradientProjectionSearch + { + public static GradientProjectionResult Search(Vector x0, Vector gradient, Matrix hessian, Vector lowerBound, Vector upperBound) + { + List isFixed = new List(x0.Count); + List breakpoint = new List(x0.Count); + for (int ii = 0; ii < x0.Count; ++ii) + { + breakpoint.Add(0.0); + isFixed.Add(false); + if (gradient[ii] < 0) + breakpoint[ii] = (x0[ii] - upperBound[ii]) / gradient[ii]; + else if (gradient[ii] > 0) + breakpoint[ii] = (x0[ii] - lowerBound[ii]) / gradient[ii]; + else + { + if (Math.Abs(x0[ii] - upperBound[ii]) < 100 * Double.Epsilon || Math.Abs(x0[ii] - lowerBound[ii]) < 100 * Double.Epsilon) + breakpoint[ii] = 0.0; + else + breakpoint[ii] = Double.PositiveInfinity; + } + } + + var orderedBreakpoint = new List(x0.Count); + orderedBreakpoint.AddRange(breakpoint); + orderedBreakpoint.Sort(); + + // Compute initial state variables + var d = -gradient; + for (int ii = 0; ii < d.Count; ++ii) + if (breakpoint[ii] <= 0.0) + d[ii] *= 0.0; + + + int jj = -1; + var x = x0; + var f1 = gradient * d; + var f2 = 0.5 * d * hessian * d; + var sMin = -f1 / f2; + var maxS = orderedBreakpoint[0]; + + if (sMin < maxS) + return new GradientProjectionResult(x + sMin * d, 0,isFixed); + + // while minimum of the last quadratic piece observed is beyond the interval searched + while (true) + { + // update data to the beginning of the interval we're searching + jj += 1; + x = x + d * maxS; + maxS = orderedBreakpoint[jj+1] - orderedBreakpoint[jj]; + + int fixedCount = 0; + for (int ii = 0; ii < d.Count; ++ii) + if (orderedBreakpoint[jj] >= breakpoint[ii]) + { + d[ii] *= 0.0; + isFixed[ii] = true; + fixedCount += 1; + } + + if (Double.IsPositiveInfinity(orderedBreakpoint[jj + 1])) + return new GradientProjectionResult(x, fixedCount, isFixed); + + f1 = gradient * d + (x - x0) * hessian * d; + f2 = d * hessian * d; + + sMin = -f1 / f2; + + if (sMin < maxS) + return new GradientProjectionResult(x + sMin * d, fixedCount, isFixed); + else if (jj + 1 >= orderedBreakpoint.Count - 1) + { + isFixed[isFixed.Count - 1] = true; + return new GradientProjectionResult(x + maxS * d, lowerBound.Count, isFixed); + } + } + } + + public struct GradientProjectionResult + { + public GradientProjectionResult(Vector cauchyPoint, int fixedCount, List isFixed) + { + CauchyPoint = cauchyPoint; + FixedCount = fixedCount; + IsFixed = isFixed; + } + public Vector CauchyPoint { get; } + public int FixedCount { get; } + public List IsFixed { get; } + } + } +} diff --git a/src/TestData/TestData.csproj b/src/TestData/TestData.csproj index aa3d569c..2bfc155d 100644 --- a/src/TestData/TestData.csproj +++ b/src/TestData/TestData.csproj @@ -26,7 +26,7 @@ AllRules.ruleset 1591 false - 5 + 6 DEBUG;TRACE @@ -41,7 +41,7 @@ AnyCPU 1591 false - 5 + 6 ..\..\out\test-signed\Net35\ @@ -57,7 +57,7 @@ AllRules.ruleset 1591 false - 5 + 6 diff --git a/src/UnitTests/OptimizationTests/BfgsBMinimizerTests.cs b/src/UnitTests/OptimizationTests/BfgsBMinimizerTests.cs new file mode 100644 index 00000000..19ff2773 --- /dev/null +++ b/src/UnitTests/OptimizationTests/BfgsBMinimizerTests.cs @@ -0,0 +1,253 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2016 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; +using MathNet.Numerics.LinearAlgebra.Double; +using MathNet.Numerics.Optimization; +using NUnit.Framework; +using System.Linq; +using System.Text; +using System.Collections.Generic; +using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions; +using System.Collections; +using MathNet.Numerics.Optimization.ObjectiveFunctions; +using NUnit.Framework.Interfaces; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + [TestFixture] + public class BfgsBMinimizerTests + { + [Test] + public void FindMinimum_Rosenbrock_Easy() + { + var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient); + var solver = new BfgsBMinimizer (1e-5, 1e-5, 1e-5, maximumIterations: 1000); + var lowerBound = new DenseVector(new[]{ -5.0, -5.0 }); + var upperBound = new DenseVector(new[]{ 5.0, 5.0 }); + var initialGuess = new DenseVector(new[] { 1.2, 1.2 }); + + var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - RosenbrockFunction.Minimum[0]), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_Rosenbrock_Hard() + { + var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient); + var solver = new BfgsBMinimizer (1e-5, 1e-5, 1e-5, maximumIterations: 1000); + + var lowerBound = new DenseVector(new[]{ -5.0, -5.0 }); + var upperBound = new DenseVector(new[]{ 5.0, 5.0 }); + + var initialGuess = new DenseVector (new[]{ -1.2, 1.0 }); + + var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - RosenbrockFunction.Minimum[0]), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_Rosenbrock_Overton() + { + var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient); + var solver = new BfgsBMinimizer (1e-5, 1e-5, 1e-5, maximumIterations: 1000); + + var lowerBound = new DenseVector(new[]{ -5.0, -5.0 }); + var upperBound = new DenseVector(new[]{ 5.0, 5.0 }); + var initialGuess = new DenseVector (new[]{ -0.9, -0.5 }); + + var result = solver.FindMinimum (obj, lowerBound, upperBound, initialGuess); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - RosenbrockFunction.Minimum[0]), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_Rosenbrock_Easy_OneBoundary() + { + var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient); + var solver = new BfgsBMinimizer (1e-5, 1e-5, 1e-5, maximumIterations: 1000); + var lowerBound = new DenseVector(new[]{ 1.0, -5.0 }); + var upperBound = new DenseVector(new[]{ 5.0, 5.0 }); + var initialGuess = new DenseVector(new[] { 1.2, 1.2 }); + + var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - RosenbrockFunction.Minimum[0]), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_Rosenbrock_Easy_TwoBoundaries() + { + var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient); + var solver = new BfgsBMinimizer (1e-5, 1e-5, 1e-5, maximumIterations: 1000); + var lowerBound = new DenseVector(new[]{ 1.0, 1.0 }); + var upperBound = new DenseVector(new[]{ 5.0, 5.0 }); + var initialGuess = new DenseVector(new[] { 1.2, 1.2 }); + + var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - RosenbrockFunction.Minimum[0]), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_Rosenbrock_MinimumGreateerOrEqualToLowerBoundary() + { + var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient); + var solver = new BfgsBMinimizer(1e-5, 1e-5, 1e-5, maximumIterations: 1000); + + var lowerBound = new DenseVector(new[] { 2, 2.0 }); + var upperBound = new DenseVector(new[] { 5.0, 5.0 }); + + var initialGuess = new DenseVector(new[] { 2.5, 2.5 }); + + var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess); + + Assert.GreaterOrEqual(result.MinimizingPoint[0],lowerBound[0]); + Assert.GreaterOrEqual(result.MinimizingPoint[1], lowerBound[1]); + } + + [Test] + public void FindMinimum_Rosenbrock_MinimumLesserOrEqualToUpperBoundary() + { + var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient); + var solver = new BfgsBMinimizer(1e-5, 1e-5, 1e-5, maximumIterations: 1000); + + var lowerBound = new DenseVector(new[] { -2.0, -2.0 }); + var upperBound = new DenseVector(new[] { 0.5, 0.5 }); + + var initialGuess = new DenseVector(new[] { -0.9, -0.5 }); + + var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess); + + Assert.LessOrEqual(result.MinimizingPoint[0],upperBound[0]); + Assert.LessOrEqual(result.MinimizingPoint[1],upperBound[1]); + } + + [Test] + [TestCaseSource(typeof(MghTestCaseEnumerator))] + public void Mgh_Tests(TestFunctions.TestCase test_case) + { + var obj = new MghObjectiveFunction(test_case.Function, true, true); + var solver = new BfgsBMinimizer(1e-8, 1e-8, 1e-8, 1000); + + var result = solver.FindMinimum(obj, test_case.LowerBound, test_case.UpperBound, test_case.InitialGuess); + + if (test_case.MinimizingPoint != null) + { + Assert.That((result.MinimizingPoint - test_case.MinimizingPoint).L2Norm(), Is.LessThan(1e-3)); + } + + var val1 = result.FunctionInfoAtMinimum.Value; + var val2 = test_case.MinimalValue; + var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2)); + var abs_err = Math.Abs(val1 - val2); + var rel_err = abs_err / abs_min; + var success = (abs_min <= 1 && abs_err < 1e-3) || (abs_min > 1 && rel_err < 1e-3); + Assert.That(success, "Minimal function value is not as expected."); + } + + [Test] + [TestCaseSource(typeof(FdMghTestCaseEnumerator))] + public void Mgh_FiniteDifference_Tests(TestFunctions.TestCase test_case) + { + var obj1 = new MghObjectiveFunction(test_case.Function, true, true); + var obj = new ForwardDifferenceGradientObjectiveFunction(obj1, test_case.LowerBound, test_case.UpperBound, 1e-10, 1e-10); + var solver = new BfgsBMinimizer(1e-8, 1e-8, 1e-8, 1000); + + var result = solver.FindMinimum(obj, test_case.LowerBound, test_case.UpperBound, test_case.InitialGuess); + + if (test_case.MinimizingPoint != null) + { + Assert.That((result.MinimizingPoint - test_case.MinimizingPoint).L2Norm(), Is.LessThan(1e-3)); + } + + var val1 = result.FunctionInfoAtMinimum.Value; + var val2 = test_case.MinimalValue; + var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2)); + var abs_err = Math.Abs(val1 - val2); + var rel_err = abs_err / abs_min; + var success = (abs_min <= 1 && abs_err < 1e-3) || (abs_min > 1 && rel_err < 1e-3); + Assert.That(success, "Minimal function value is not as expected."); + } + + + private class BaseMghTestCaseEnumerator : IEnumerable + { + private string _prefix = ""; + + public BaseMghTestCaseEnumerator(string prefix) + { + if (prefix.EndsWith(" ")) + _prefix = prefix; + else + _prefix = prefix + " "; + } + public IEnumerator GetEnumerator() + { + return + RosenbrockFunction2.TestCases + .Concat(BealeFunction.TestCases) + .Concat(HelicalValleyFunction.TestCases) + .Concat(MeyerFunction.TestCases) + .Concat(PowellSingularFunction.TestCases) + .Concat(WoodFunction.TestCases) + .Concat(BrownAndDennisFunction.TestCases) + .Where(x => x.IsBounded) + .Select(x => new TestCaseData(x) + .SetName(_prefix + x.FullName) + ) + .GetEnumerator(); + } + + IEnumerator IEnumerable.GetEnumerator() + { + return this.GetEnumerator(); + } + } + + private class MghTestCaseEnumerator : BaseMghTestCaseEnumerator + { + public MghTestCaseEnumerator() : base("") { } + } + + private class FdMghTestCaseEnumerator : BaseMghTestCaseEnumerator + { + public FdMghTestCaseEnumerator() : base("FD") { } + } + } +} + diff --git a/src/UnitTests/OptimizationTests/BfgsMinimizerTests.cs b/src/UnitTests/OptimizationTests/BfgsMinimizerTests.cs new file mode 100644 index 00000000..8622886f --- /dev/null +++ b/src/UnitTests/OptimizationTests/BfgsMinimizerTests.cs @@ -0,0 +1,131 @@ +using System; +using System.Linq; +using MathNet.Numerics.LinearAlgebra.Double; +using MathNet.Numerics.Optimization; +using NUnit.Framework; +using System.Text; +using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions; +using System.Collections.Generic; +using System.Collections; +using NUnit.Framework.Interfaces; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + [TestFixture] + public class BfgsMinimizerTests + { + [Test] + public void FindMinimum_Rosenbrock_Easy() + { + var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient); + var solver = new BfgsMinimizer(1e-5, 1e-5, 1e-5, 1000); + var result = solver.FindMinimum(obj, new DenseVector(new[] { 1.2, 1.2 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - RosenbrockFunction.Minimum[0]), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_Rosenbrock_Hard() + { + var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient); + var solver = new BfgsMinimizer(1e-5, 1e-5, 1e-5, 1000); + var result = solver.FindMinimum(obj, new DenseVector(new[] { -1.2, 1.0 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - RosenbrockFunction.Minimum[0]), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_Rosenbrock_Overton() + { + var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient); + var solver = new BfgsMinimizer(1e-5, 1e-5, 1e-5, 1000); + var result = solver.FindMinimum(obj, new DenseVector(new[] { -0.9, -0.5 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - RosenbrockFunction.Minimum[0]), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_BigRosenbrock_Easy() + { + var obj = ObjectiveFunction.Gradient(BigRosenbrockFunction.Value, BigRosenbrockFunction.Gradient); + var solver = new BfgsMinimizer(1e-10, 1e-5, 1e-5, 1000); + var result = solver.FindMinimum(obj, new DenseVector(new[] { 1.2*100.0, 1.2*100.0 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - BigRosenbrockFunction.Minimum[0]), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - BigRosenbrockFunction.Minimum[1]), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_BigRosenbrock_Hard() + { + var obj = ObjectiveFunction.Gradient(BigRosenbrockFunction.Value, BigRosenbrockFunction.Gradient); + var solver = new BfgsMinimizer(1e-5, 1e-5, 1e-5, 1000); + var result = solver.FindMinimum(obj, new DenseVector(new[] { -1.2*100.0, 1.0*100.0 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - BigRosenbrockFunction.Minimum[0]), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - BigRosenbrockFunction.Minimum[1]), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_BigRosenbrock_Overton() + { + var obj = ObjectiveFunction.Gradient(BigRosenbrockFunction.Value, BigRosenbrockFunction.Gradient); + var solver = new BfgsMinimizer(1e-5, 1e-5, 1e-5, 1000); + var result = solver.FindMinimum(obj, new DenseVector(new[] { -0.9*100.0, -0.5*100.0 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - BigRosenbrockFunction.Minimum[0]), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - BigRosenbrockFunction.Minimum[1]), Is.LessThan(1e-3)); + } + + private class MghTestCaseEnumerator : IEnumerable + { + public IEnumerator GetEnumerator() + { + return + RosenbrockFunction2.TestCases + .Concat(BealeFunction.TestCases) + .Concat(HelicalValleyFunction.TestCases) + .Concat(MeyerFunction.TestCases) + .Concat(PowellSingularFunction.TestCases) + .Concat(WoodFunction.TestCases) + .Concat(BrownAndDennisFunction.TestCases) + .Where(x => x.IsUnbounded) + .Select(x => new TestCaseData(x) + .SetName(x.FullName) + ) + .GetEnumerator(); + } + + IEnumerator IEnumerable.GetEnumerator() + { + return this.GetEnumerator(); + } + } + + [Test] + [TestCaseSource(typeof(MghTestCaseEnumerator))] + public void Mgh_Tests(TestFunctions.TestCase test_case) + { + var obj = new MghObjectiveFunction(test_case.Function, true, true); + var solver = new BfgsMinimizer(1e-8, 1e-8, 1e-8, 1000); + + var result = solver.FindMinimum(obj, test_case.InitialGuess); + + if (test_case.MinimizingPoint != null) + { + Assert.That((result.MinimizingPoint - test_case.MinimizingPoint).L2Norm(), Is.LessThan(1e-3)); + } + + var val1 = result.FunctionInfoAtMinimum.Value; + var val2 = test_case.MinimalValue; + var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2)); + var abs_err = Math.Abs(val1 - val2); + var rel_err = abs_err / abs_min; + var success = (abs_min <= 1 && abs_err < 1e-3) || (abs_min > 1 && rel_err < 1e-3); + Assert.That(success, "Minimal function value is not as expected."); + } + } +} diff --git a/src/UnitTests/OptimizationTests/BfgsTest.cs b/src/UnitTests/OptimizationTests/BfgsTest.cs new file mode 100644 index 00000000..b2c17916 --- /dev/null +++ b/src/UnitTests/OptimizationTests/BfgsTest.cs @@ -0,0 +1,56 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2016 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using MathNet.Numerics.LinearAlgebra.Double; +using MathNet.Numerics.Optimization; +using NUnit.Framework; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + [TestFixture, Category("RootFinding")] + internal class BfgsTest + { + private const double Precision = 1e-4; + + [Test] + public void MinimizeRosenbrock() + { + CheckRosenbrock(15.0, 8.0, expectedMin: 0.0); + CheckRosenbrock(-1.2, 1.0, expectedMin: 0.0); + CheckRosenbrock(-1.2, 100.0, expectedMin: 0.0); + } + + private static void CheckRosenbrock(double a, double b, double expectedMin) + { + var x = BfgsSolver.Solve(new DenseVector(new[] { a, b }), RosenbrockFunction.Value, RosenbrockFunction.Gradient); + Numerics.Precision.AlmostEqual(expectedMin, RosenbrockFunction.Value(x), Precision); + } + } +} diff --git a/src/UnitTests/OptimizationTests/ConjugateGradientMinimizerTests.cs b/src/UnitTests/OptimizationTests/ConjugateGradientMinimizerTests.cs new file mode 100644 index 00000000..254a1db7 --- /dev/null +++ b/src/UnitTests/OptimizationTests/ConjugateGradientMinimizerTests.cs @@ -0,0 +1,101 @@ +using System; +using MathNet.Numerics.LinearAlgebra.Double; +using MathNet.Numerics.Optimization; +using NUnit.Framework; +using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions; +using System.Collections; +using System.Collections.Generic; +using System.Linq; +using NUnit.Framework.Interfaces; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + [TestFixture] + public class ConjugateGradientMinimizerTests + { + [Test] + public void FindMinimum_Rosenbrock_Easy() + { + var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient); + var solver = new ConjugateGradientMinimizer(1e-5, 1000); + var result = solver.FindMinimum(obj, new DenseVector(new[]{1.2,1.2})); + + Assert.That(Math.Abs(result.MinimizingPoint[0]-1.0), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_Rosenbrock_Hard() + { + var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient); + var solver = new ConjugateGradientMinimizer(1e-5, 1000); + var result = solver.FindMinimum(obj, new DenseVector(new[] { -1.2, 1.0 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3)); + } + + private class MghTestCaseEnumerator : IEnumerable + { + private static readonly string[] _ignore_list = + { + "Beale fun (MGH #5) unbounded", + "Meyer fun (MGH #10) unbounded", + "Powell singular fun (MGH #13) unbounded", + "Rosenbrock fun (MGH #1) hard start", + "Rosenbrock fun (MGH #1) Overton start", + }; + + private static bool in_ignore_list(string test_name) + { + return _ignore_list.Contains(test_name); + } + + public IEnumerator GetEnumerator() + { + return + RosenbrockFunction2.TestCases + .Concat(BealeFunction.TestCases) + .Concat(HelicalValleyFunction.TestCases) + .Concat(MeyerFunction.TestCases) + .Concat(PowellSingularFunction.TestCases) + .Concat(WoodFunction.TestCases) + .Concat(BrownAndDennisFunction.TestCases) + .Where(x => x.IsUnbounded) + .Select(x => new TestCaseData(x) + .SetName(x.FullName) + .IgnoreIf(in_ignore_list(x.FullName),"Algo error, not implementation error.") + ) + .GetEnumerator(); + } + + IEnumerator IEnumerable.GetEnumerator() + { + return this.GetEnumerator(); + } + } + + [Test] + [TestCaseSource(typeof(MghTestCaseEnumerator))] + public void Mgh_Tests(TestFunctions.TestCase test_case) + { + var obj = new MghObjectiveFunction(test_case.Function, true, true); + var solver = new ConjugateGradientMinimizer(1e-8, 1000); + + var result = solver.FindMinimum(obj, test_case.InitialGuess); + + if (test_case.MinimizingPoint != null) + { + Assert.That((result.MinimizingPoint - test_case.MinimizingPoint).L2Norm(), Is.LessThan(1e-3)); + } + + var val1 = result.FunctionInfoAtMinimum.Value; + var val2 = test_case.MinimalValue; + var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2)); + var abs_err = Math.Abs(val1 - val2); + var rel_err = abs_err / abs_min; + var success = (abs_min <= 1 && abs_err < 1e-3) || (abs_min > 1 && rel_err < 1e-3); + Assert.That(success, "Minimal function value is not as expected."); + } + } +} diff --git a/src/UnitTests/OptimizationTests/GoldenSectionMinimizerTests.cs b/src/UnitTests/OptimizationTests/GoldenSectionMinimizerTests.cs new file mode 100644 index 00000000..e04088d9 --- /dev/null +++ b/src/UnitTests/OptimizationTests/GoldenSectionMinimizerTests.cs @@ -0,0 +1,32 @@ +using System; +using MathNet.Numerics.Optimization; +using NUnit.Framework; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + [TestFixture] + public class GoldenSectionMinimizerTests + { + [Test] + public void Test_Works() + { + var algorithm = new GoldenSectionMinimizer(1e-5, 1000); + var f1 = new Func(x => (x - 3)*(x - 3)); + var obj = new SimpleObjectiveFunction1D(f1); + var r1 = algorithm.FindMinimum(obj, -100, 100); + + Assert.That(Math.Abs(r1.MinimizingPoint - 3.0), Is.LessThan(1e-4)); + } + + [Test] + public void Test_ExpansionWorks() + { + var algorithm = new GoldenSectionMinimizer(1e-5, 1000); + var f1 = new Func(x => (x - 3)*(x - 3)); + var obj = new SimpleObjectiveFunction1D(f1); + var r1 = algorithm.FindMinimum(obj, -5, 5); + + Assert.That(Math.Abs(r1.MinimizingPoint - 3.0), Is.LessThan(1e-4)); + } + } +} diff --git a/src/UnitTests/OptimizationTests/NelderMeadSimplexTests.cs b/src/UnitTests/OptimizationTests/NelderMeadSimplexTests.cs new file mode 100644 index 00000000..ef16d76f --- /dev/null +++ b/src/UnitTests/OptimizationTests/NelderMeadSimplexTests.cs @@ -0,0 +1,133 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2016 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using MathNet.Numerics.LinearAlgebra.Double; +using MathNet.Numerics.Optimization; +using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions; +using NUnit.Framework; +using System; +using System.Collections; +using System.Collections.Generic; +using System.Linq; +using NUnit.Framework.Interfaces; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + [TestFixture] + public class NelderMeadSimplexTests + { + [Test] + public void NMS_FindMinimum_Rosenbrock_Easy() + { + var obj = ObjectiveFunction.Value(RosenbrockFunction.Value); + var solver = new NelderMeadSimplex(1e-5, maximumIterations: 1000); + var initialGuess = new DenseVector(new[] { 1.2, 1.2 }); + + var result = solver.FindMinimum(obj, initialGuess); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3)); + } + + + [Test] + public void NMS_FindMinimum_Rosenbrock_Hard() + { + var obj = ObjectiveFunction.Value(RosenbrockFunction.Value); + var solver = new NelderMeadSimplex(1e-5, maximumIterations: 1000); + + var initialGuess = new DenseVector(new[] { -1.2, 1.0 }); + + var result = solver.FindMinimum(obj,initialGuess); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3)); + } + + private class MghTestCaseEnumerator : IEnumerable + { + private static readonly string[] _ignore_list = + { + "Meyer fun (MGH #10) unbounded", + }; + + private static bool in_ignore_list(string test_name) + { + return _ignore_list.Contains(test_name); + } + + public IEnumerator GetEnumerator() + { + return + RosenbrockFunction2.TestCases + .Concat(BealeFunction.TestCases) + .Concat(HelicalValleyFunction.TestCases) + .Concat(MeyerFunction.TestCases) + .Concat(PowellSingularFunction.TestCases) + .Concat(WoodFunction.TestCases) + .Concat(BrownAndDennisFunction.TestCases) + .Where(x => x.IsUnbounded) + .Select(x => new TestCaseData(x) + .SetName(x.FullName) + .IgnoreIf(in_ignore_list(x.FullName), "Algo error, not implementation error") + ) + .GetEnumerator(); + } + + IEnumerator IEnumerable.GetEnumerator() + { + return this.GetEnumerator(); + } + } + + [Test] + [TestCaseSource(typeof(MghTestCaseEnumerator))] + public void Mgh_Tests(TestFunctions.TestCase test_case) + { + var obj = new MghObjectiveFunction(test_case.Function, true, true); + var solver = new NelderMeadSimplex(1e-8, 1000); + + var result = solver.FindMinimum(obj, test_case.InitialGuess); + + if (test_case.MinimizingPoint != null) + { + Assert.That((result.MinimizingPoint - test_case.MinimizingPoint).L2Norm(), Is.LessThan(1e-3)); + } + + var val1 = result.FunctionInfoAtMinimum.Value; + var val2 = test_case.MinimalValue; + var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2)); + var abs_err = Math.Abs(val1 - val2); + var rel_err = abs_err / abs_min; + var success = (abs_min <= 1 && abs_err < 1e-3) || (abs_min > 1 && rel_err < 1e-3); + Assert.That(success, "Minimal function value is not as expected."); + } + } +} diff --git a/src/UnitTests/OptimizationTests/NewtonMinimizerTests.cs b/src/UnitTests/OptimizationTests/NewtonMinimizerTests.cs new file mode 100644 index 00000000..6f66c2c7 --- /dev/null +++ b/src/UnitTests/OptimizationTests/NewtonMinimizerTests.cs @@ -0,0 +1,188 @@ +using System; +using MathNet.Numerics.LinearAlgebra.Double; +using MathNet.Numerics.Optimization; +using MathNet.Numerics.Optimization.ObjectiveFunctions; +using NUnit.Framework; +using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions; +using System.Collections.Generic; +using System.Collections; +using System.Linq; +using NUnit.Framework.Interfaces; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + public class LazyRosenbrockObjectiveFunction : LazyObjectiveFunctionBase + { + public LazyRosenbrockObjectiveFunction() : base(true, true) { } + + public override IObjectiveFunction CreateNew() + { + return new LazyRosenbrockObjectiveFunction(); + } + + protected override void EvaluateValue() + { + Value = RosenbrockFunction.Value(Point); + } + + protected override void EvaluateGradient() + { + Gradient = RosenbrockFunction.Gradient(Point); + } + + protected override void EvaluateHessian() + { + Hessian = RosenbrockFunction.Hessian(Point); + } + } + + public class RosenbrockObjectiveFunction : ObjectiveFunctionBase + { + public RosenbrockObjectiveFunction() : base(true, true) { } + + public override IObjectiveFunction CreateNew() + { + return new RosenbrockObjectiveFunction(); + } + + protected override void Evaluate() + { + // here we could directly overwrite the existing matrix cells instead. + // note: values must then be initialized manually first, if null. + Value = RosenbrockFunction.Value(Point); + Gradient = RosenbrockFunction.Gradient(Point); + Hessian = RosenbrockFunction.Hessian(Point); + } + } + + [TestFixture] + public class NewtonMinimizerTests + { + [Test] + public void FindMinimum_Rosenbrock_Easy() + { + var obj = ObjectiveFunction.GradientHessian(RosenbrockFunction.Value, RosenbrockFunction.Gradient, RosenbrockFunction.Hessian); + var solver = new NewtonMinimizer(1e-5, 1000); + var result = solver.FindMinimum(obj, new DenseVector(new[] { 1.2, 1.2 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_Rosenbrock_Hard() + { + var obj = ObjectiveFunction.GradientHessian(point => Tuple.Create(RosenbrockFunction.Value(point), RosenbrockFunction.Gradient(point), RosenbrockFunction.Hessian(point))); + var solver = new NewtonMinimizer(1e-5, 1000); + var result = solver.FindMinimum(obj, new DenseVector(new[] { -1.2, 1.0 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_Rosenbrock_Overton() + { + var obj = new LazyRosenbrockObjectiveFunction(); + var solver = new NewtonMinimizer(1e-5, 1000); + var result = solver.FindMinimum(obj, new DenseVector(new[] { -0.9, -0.5 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_Linesearch_Rosenbrock_Easy() + { + var obj = new RosenbrockObjectiveFunction(); + var solver = new NewtonMinimizer(1e-5, 1000, true); + var result = solver.FindMinimum(obj, new DenseVector(new[] { 1.2, 1.2 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_Linesearch_Rosenbrock_Hard() + { + var obj = new LazyRosenbrockObjectiveFunction(); + var solver = new NewtonMinimizer(1e-5, 1000, true); + var result = solver.FindMinimum(obj, new DenseVector(new[] { -1.2, 1.0 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_Linesearch_Rosenbrock_Overton() + { + var obj = new LazyRosenbrockObjectiveFunction(); + var solver = new NewtonMinimizer(1e-5, 1000, true); + var result = solver.FindMinimum(obj, new DenseVector(new[] { -0.9, -0.5 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3)); + } + + private class MghTestCaseEnumerator : IEnumerable + { + private static readonly string[] _ignore_list = + { + "Beale fun (MGH #5) unbounded", + "Meyer fun (MGH #10) unbounded", + "Wood fun (MGH #14) unbounded", + }; + + private static bool in_ignore_list(string test_name) + { + return _ignore_list.Contains(test_name); + } + + public IEnumerator GetEnumerator() + { + return + RosenbrockFunction2.TestCases + .Concat(BealeFunction.TestCases) + .Concat(HelicalValleyFunction.TestCases) + .Concat(MeyerFunction.TestCases) + .Concat(PowellSingularFunction.TestCases) + .Concat(WoodFunction.TestCases) + .Concat(BrownAndDennisFunction.TestCases) + .Where(x => x.IsUnbounded) + .Select(x => new TestCaseData(x) + .SetName(x.FullName) + .IgnoreIf(in_ignore_list(x.FullName),"Algo error, not implementation error") + ) + .GetEnumerator(); + } + + IEnumerator IEnumerable.GetEnumerator() + { + return this.GetEnumerator(); + } + } + + [Test] + [TestCaseSource(typeof(MghTestCaseEnumerator))] + public void Mgh_Tests(TestFunctions.TestCase test_case) + { + var obj = new MghObjectiveFunction(test_case.Function, true, true); + var solver = new NewtonMinimizer(1e-8, 1000, useLineSearch: false); + + var result = solver.FindMinimum(obj, test_case.InitialGuess); + + if (test_case.MinimizingPoint != null) + { + Assert.That((result.MinimizingPoint - test_case.MinimizingPoint).L2Norm(), Is.LessThan(1e-3)); + } + + var val1 = result.FunctionInfoAtMinimum.Value; + var val2 = test_case.MinimalValue; + var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2)); + var abs_err = Math.Abs(val1 - val2); + var rel_err = abs_err / abs_min; + var success = (abs_min <= 1 && abs_err < 1e-3) || (abs_min > 1 && rel_err < 1e-3); + Assert.That(success, "Minimal function value is not as expected."); + } + } +} diff --git a/src/UnitTests/OptimizationTests/RosenbrockFunction.cs b/src/UnitTests/OptimizationTests/RosenbrockFunction.cs new file mode 100644 index 00000000..56a43f8d --- /dev/null +++ b/src/UnitTests/OptimizationTests/RosenbrockFunction.cs @@ -0,0 +1,66 @@ +using System; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.LinearAlgebra.Double; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + public static class RosenbrockFunction + { + public static double Value(Vector input) + { + return Math.Pow((1 - input[0]), 2) + 100 * Math.Pow((input[1] - input[0] * input[0]), 2); + } + + public static Vector Gradient(Vector input) + { + Vector output = new DenseVector(2); + output[0] = -2 * (1 - input[0]) + 200 * (input[1] - input[0] * input[0]) * (-2 * input[0]); + output[1] = 2 * 100 * (input[1] - input[0] * input[0]); + return output; + } + + public static Matrix Hessian(Vector input) + { + Matrix output = new DenseMatrix(2, 2); + output[0, 0] = 2 - 400 * input[1] + 1200 * input[0] * input[0]; + output[1, 1] = 200; + output[0, 1] = -400 * input[0]; + output[1, 0] = output[0, 1]; + return output; + } + + public static Vector Minimum + { + get + { + return new DenseVector(new double[] { 1, 1 }); + } + } + } + + public static class BigRosenbrockFunction + { + public static double Value(Vector input) + { + return 1000.0 + 100.0 * RosenbrockFunction.Value(input / 100.0); + } + + public static Vector Gradient(Vector input) + { + return 100.0 * RosenbrockFunction.Gradient(input / 100.0); + } + + public static Matrix Hessian(Vector input) + { + return 100.0 * RosenbrockFunction.Hessian(input / 100.0); + } + + public static Vector Minimum + { + get + { + return new DenseVector(new double[] { 100, 100 }); + } + } + } +} diff --git a/src/UnitTests/OptimizationTests/RosenbrockFunctionTests.cs b/src/UnitTests/OptimizationTests/RosenbrockFunctionTests.cs new file mode 100644 index 00000000..b1b0e5ae --- /dev/null +++ b/src/UnitTests/OptimizationTests/RosenbrockFunctionTests.cs @@ -0,0 +1,55 @@ +using System; +using MathNet.Numerics.LinearAlgebra.Double; +using NUnit.Framework; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + [TestFixture] + class RosenbrockFunctionTests + { + [Test] + public void TestGradient() + { + var input = new DenseVector(new[]{ -0.9, -0.5 } ); + + var v1 = RosenbrockFunction.Value(input); + var g = RosenbrockFunction.Gradient(input); + + var eps = 1e-5; + var eps0 = (new DenseVector(new[] { 1.0, 0.0 })) * eps; + var eps1 = (new DenseVector(new[] { 0.0, 1.0 })) * eps; + + var g0 = (RosenbrockFunction.Value(input + eps0) - RosenbrockFunction.Value(input - eps0)) / (2 * eps); + var g1 = (RosenbrockFunction.Value(input + eps1) - RosenbrockFunction.Value(input - eps1)) / (2 * eps); + + Assert.That(Math.Abs(g0 - g[0]) < 1e-3); + Assert.That(Math.Abs(g1 - g[1]) < 1e-3); + } + + [Test] + public void TestHessian() + { + var input = new DenseVector(new[] { -0.9, -0.5 }); + + var v1 = RosenbrockFunction.Value(input); + var h = RosenbrockFunction.Hessian(input); + + var eps = 1e-5; + + var eps0 = (new DenseVector(new[] { 1.0, 0.0 })) * eps; + var eps1 = (new DenseVector(new[] { 0.0, 1.0 })) * eps; + + var epsuu = (new DenseVector(new[] { 1.0, 1.0 })) * eps; + var epsud = (new DenseVector(new[] { 1.0, -1.0 })) * eps; + + var h00 = (RosenbrockFunction.Value(input + eps0) - 2*RosenbrockFunction.Value(input) + RosenbrockFunction.Value(input - eps0)) / (eps*eps); + var h11 = (RosenbrockFunction.Value(input + eps1) - 2 * RosenbrockFunction.Value(input) + RosenbrockFunction.Value(input - eps1)) / (eps * eps); + var h01 = (RosenbrockFunction.Value(input + epsuu) - RosenbrockFunction.Value(input + epsud) - RosenbrockFunction.Value(input - epsud) + RosenbrockFunction.Value(input - epsuu)) / (4*eps * eps); + + Assert.That(Math.Abs(h00 - h[0,0]) < 1e-3); + Assert.That(Math.Abs(h11 - h[1,1]) < 1e-3); + Assert.That(Math.Abs(h01 - h[0, 1]) < 1e-3); + Assert.That(Math.Abs(h01 - h[1, 0]) < 1e-3); + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestCaseDataExtensions.cs b/src/UnitTests/OptimizationTests/TestCaseDataExtensions.cs new file mode 100644 index 00000000..91d2d4ef --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestCaseDataExtensions.cs @@ -0,0 +1,20 @@ +using NUnit.Framework; +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + internal static class TestCaseDataExtensions + { + public static TestCaseData IgnoreIf(this TestCaseData input, bool do_ignore, string reason) + { + if (do_ignore) + return input.Ignore(reason); + else + return input; + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctionAdapters.cs b/src/UnitTests/OptimizationTests/TestFunctionAdapters.cs new file mode 100644 index 00000000..8b5a43ae --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctionAdapters.cs @@ -0,0 +1,54 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.Optimization; +using MathNet.Numerics.Optimization.ObjectiveFunctions; +using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions; +using MathNet.Numerics.LinearAlgebra.Double; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + public class MghObjectiveFunction : LazyObjectiveFunctionBase + { + private ITestFunction TestFunction; + + public MghObjectiveFunction(ITestFunction testFunction, bool use_gradient, bool use_hessian) + : base(use_gradient, use_hessian) + { + this.TestFunction = testFunction; + } + + public override IObjectiveFunction CreateNew() + { + return new MghObjectiveFunction(this.TestFunction, this.IsGradientSupported, this.IsHessianSupported); + } + + protected override void EvaluateValue() + { + this.Value = this.TestFunction.SsqValue(this.Point); + } + + protected override void EvaluateGradient() + { + if (this.IsGradientSupported) + { + if (this.GradientValue == null) + this.Gradient = new DenseVector(this.TestFunction.ParameterDimension); + this.TestFunction.SsqGradientByRef(this.Point, GradientValue); + } + } + + protected override void EvaluateHessian() + { + if (this.IsHessianSupported) + { + if (this.HessianValue == null) + this.Hessian = new DenseMatrix(this.TestFunction.ParameterDimension, this.TestFunction.ParameterDimension); + this.TestFunction.SsqHessianByRef(this.Point, HessianValue); + } + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctionTests.cs b/src/UnitTests/OptimizationTests/TestFunctionTests.cs new file mode 100644 index 00000000..e3108615 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctionTests.cs @@ -0,0 +1,309 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions; +using NUnit.Framework; +using MathNet.Numerics.LinearAlgebra; +using System.Collections; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + [TestFixture] + public class TestFunctionTests + { + private static IEnumerable MghCases + { + get + { + return Enumerable.Empty() + .Concat(RosenbrockFunction2.TestCases) + .Concat(BealeFunction.TestCases) + .Concat(HelicalValleyFunction.TestCases) + .Concat(MeyerFunction.TestCases) + .Concat(PowellSingularFunction.TestCases) + .Concat(WoodFunction.TestCases) + .Concat(BrownAndDennisFunction.TestCases); + } + } + + private class MghCaseEnumerator : IEnumerable + { + public string CategoryName { get; protected set; } + + public MghCaseEnumerator(string category_name) + { + this.CategoryName = category_name; + } + + public virtual IEnumerator GetEnumerator() + { + return MghCases + .Select(x => + new TestCaseData(x) + .SetName($"{x.FullName} {this.CategoryName}") + ).GetEnumerator(); + } + + IEnumerator IEnumerable.GetEnumerator() + { + return this.GetEnumerator(); + } + } + + [Test] + public void Smoke_Construction() + { + var c = new TestCase() + { + InitialGuess = new double[] { 1, 2, 3 }, + MinimizingPoint = new double[] { 1, 1, 1 }, + MinimalValue = 0 + }; + } + + private class ValueAtMinimumSource : MghCaseEnumerator + { + public ValueAtMinimumSource() : base("ValueAtMinimum") { } + + public override IEnumerator GetEnumerator() + { + return MghCases + .Where(x => x.MinimizingPoint != null) + .Select(x => + new TestCaseData(x) + .SetName($"{x.FullName} {this.CategoryName}") + ) + .GetEnumerator(); + } + } + + [Test] + [TestCaseSource(typeof(ValueAtMinimumSource))] + public void ValueAtMinimum(TestFunctions.TestCase test_case) + { + if (test_case.MinimizingPoint != null) + { + var value_at_minimum = test_case.Function.SsqValue(test_case.MinimizingPoint); + Assert.That( + Math.Abs(value_at_minimum - test_case.MinimalValue) < 1e-3, + $"Function value at minimum not as expected." + ); + } + } + + private class GradientAtStartSource : MghCaseEnumerator + { + public GradientAtStartSource() : base("GradientAtStart") { } + } + + [Test] + [TestCaseSource(typeof(GradientAtStartSource))] + public void GradientAtStart(TestFunctions.TestCase test_case) + { + var a_grad = test_case.Function.SsqGradient(test_case.InitialGuess); + var fd_grad = Vector.Build.Dense(test_case.Function.ParameterDimension, 0.0); + + for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii) + { + var h = 1e-6; + + var bump_up = test_case.InitialGuess.Clone(); + bump_up[ii] += h; + var bump_down = test_case.InitialGuess.Clone(); + bump_down[ii] -= h; + + var up_val = test_case.Function.SsqValue(bump_up); + var down_val = test_case.Function.SsqValue(bump_down); + + fd_grad[ii] = 0.5 * (up_val - down_val) / h; + } + + for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii) + { + var val1 = a_grad[ii]; + var val2 = fd_grad[ii]; + var min_abs_val = Math.Min(Math.Abs(val1), Math.Abs(val2)); + if (min_abs_val <= 1) + Assert.That(Math.Abs(val1 - val2) < 1e-3, $"Problem with gradient value at start point."); + else + Assert.That(Math.Abs(val1 - val2) / min_abs_val < 1e-3, $"Problem with gradient value at start point."); + } + } + + private class HessianAtStartSource : MghCaseEnumerator + { + public HessianAtStartSource() : base("HessianAtStart") { } + + public override IEnumerator GetEnumerator() + { + return MghCases + .Where(x => x.MinimizingPoint != null) + .Select(x => + new TestCaseData(x) + .SetName($"{x.FullName} {this.CategoryName}") + ) + .GetEnumerator(); + } + } + [Test] + [TestCaseSource(typeof(HessianAtStartSource))] + public void HessianAtStart(TestFunctions.TestCase test_case) + { + var a_hess = test_case.Function.SsqHessian(test_case.InitialGuess); + var fd_hess = Matrix.Build.Dense(test_case.Function.ParameterDimension, test_case.Function.ParameterDimension); + for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii) + { + for (int jj = 0; jj < test_case.Function.ParameterDimension; ++jj) + { + var h1 = 1e-3 * Math.Max(1.0, Math.Abs(test_case.InitialGuess[ii])); + var h2 = 1e-3 * Math.Max(1.0, Math.Abs(test_case.InitialGuess[jj])); + + var bump_uu = test_case.InitialGuess.Clone(); + bump_uu[ii] += h1; + bump_uu[jj] += h2; + + var bump_dd = test_case.InitialGuess.Clone(); + bump_dd[ii] -= h1; + bump_dd[jj] -= h2; + + var bump_ud = test_case.InitialGuess.Clone(); + bump_ud[ii] += h1; + bump_ud[jj] -= h2; + + var bump_du = test_case.InitialGuess.Clone(); + bump_du[ii] -= h1; + bump_du[jj] += h2; + + var val_uu = test_case.Function.SsqValue(bump_uu); + var val_dd = test_case.Function.SsqValue(bump_dd); + var val_ud = test_case.Function.SsqValue(bump_ud); + var val_du = test_case.Function.SsqValue(bump_du); + + fd_hess[ii, jj] = (val_uu - val_ud + val_dd - val_du) / (4 * h1 * h2); + } + } + + for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii) + { + for (int jj = 0; jj < test_case.Function.ParameterDimension; ++jj) + { + var val1 = fd_hess[ii, jj]; + var val2 = a_hess[ii, jj]; + + var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2)); + if (abs_min <= 1) + { + Assert.That(Math.Abs(val1 - val2) < 1e-3, $"Problem with hessian at start point."); + } + else + { + Assert.That(Math.Abs(val1 - val2) / abs_min < 0.05, $"Problem with hessian at start point."); + } + } + } + } + + private class ItemGradientAtStartSource : MghCaseEnumerator + { + public ItemGradientAtStartSource() : base("ItemGradientAtStart") { } + } + + [Test] + [TestCaseSource(typeof(ItemGradientAtStartSource))] + public void ItemGradientAtStart(TestFunctions.TestCase test_case) + { + for (var item_index = 0; item_index < test_case.Function.ItemDimension; ++item_index) + { + + var a_grad = test_case.Function.ItemGradient(test_case.InitialGuess, item_index); + var h = 1e-4; + var fd_grad = Vector.Build.Dense(test_case.Function.ParameterDimension, 0.0); + + for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii) + { + var bump_up = test_case.InitialGuess.Clone(); + bump_up[ii] += h; + var bump_down = test_case.InitialGuess.Clone(); + bump_down[ii] -= h; + + var up_val = test_case.Function.ItemValue(bump_up, item_index); + var down_val = test_case.Function.ItemValue(bump_down, item_index); + + fd_grad[ii] = 0.5 * (up_val - down_val) / h; + } + + for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii) + { + Assert.That(Math.Abs(fd_grad[ii] - a_grad[ii]) < 1e-3, $"Failed for parameter {ii}"); + } + } + } + + private class ItemHessianAtStartSource : MghCaseEnumerator + { + public ItemHessianAtStartSource() : base("ItemHessianAtStart") { } + } + + [Test] + [TestCaseSource(typeof(ItemHessianAtStartSource))] + public void ItemHessianAtStart(TestFunctions.TestCase test_case) + { + for (var item_index = 0; item_index < test_case.Function.ItemDimension; ++item_index) + { + var a_hess = test_case.Function.ItemHessian(test_case.InitialGuess, item_index); + var h = 1e-4; + var fd_hess = Matrix.Build.Dense(test_case.Function.ParameterDimension, test_case.Function.ParameterDimension); + for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii) + { + for (int jj = 0; jj < test_case.Function.ParameterDimension; ++jj) + { + var bump_uu = test_case.InitialGuess.Clone(); + bump_uu[ii] += h; + bump_uu[jj] += h; + + var bump_dd = test_case.InitialGuess.Clone(); + bump_dd[ii] -= h; + bump_dd[jj] -= h; + + var bump_ud = test_case.InitialGuess.Clone(); + bump_ud[ii] += h; + bump_ud[jj] -= h; + + var bump_du = test_case.InitialGuess.Clone(); + bump_du[ii] -= h; + bump_du[jj] += h; + + var val_uu = test_case.Function.ItemValue(bump_uu, item_index); + var val_dd = test_case.Function.ItemValue(bump_dd, item_index); + var val_ud = test_case.Function.ItemValue(bump_ud, item_index); + var val_du = test_case.Function.ItemValue(bump_du, item_index); + + fd_hess[ii, jj] = (val_uu - val_ud + val_dd - val_du) / (4 * h * h); + } + } + + for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii) + { + for (int jj = 0; jj < test_case.Function.ParameterDimension; ++jj) + { + var val1 = fd_hess[ii, jj]; + var val2 = a_hess[ii, jj]; + + var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2)); + if (abs_min <= 1) + { + Assert.That(Math.Abs(val1 - val2) < 1e-3, $"Problem with hessian at start point."); + } + else + { + Assert.That(Math.Abs(val1 - val2) / abs_min < 0.05, $"Problem with hessian at start point."); + } + } + } + } + } + + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/BaseTestFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/BaseTestFunction.cs new file mode 100644 index 00000000..60da6c61 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/BaseTestFunction.cs @@ -0,0 +1,125 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public abstract class BaseTestFunction : ITestFunction + { + public abstract string Description { get; } + public abstract int ParameterDimension { get; } + public abstract int ItemDimension { get; } + + public abstract double ItemValue(Vector x, int itemIndex); + public abstract void ItemGradientByRef(Vector x, int itemIndex, Vector output); + public abstract void ItemHessianByRef(Vector x, int itemIndex, Matrix output); + + public virtual Vector ItemGradient(Vector x, int itemIndex) + { + var output = new LinearAlgebra.Double.DenseVector(this.ParameterDimension); + this.ItemGradientByRef(x, itemIndex, output); + return output; + } + + public virtual Matrix ItemHessian(Vector x, int itemIndex) + { + var output = new LinearAlgebra.Double.DenseMatrix(this.ParameterDimension, this.ParameterDimension); + this.ItemHessianByRef(x, itemIndex, output); + return output; + } + + + public virtual void JacobianbyRef(Vector x, Matrix output) + { + for (int ii = 0; ii < this.ItemDimension; ++ii) + { + var grad = this.ItemGradient(x, ii); + output.SetRow(ii, grad); + } + } + + public virtual Matrix Jacobian(Vector x) + { + var output = new LinearAlgebra.Double.DenseMatrix(this.ItemDimension, this.ParameterDimension); + this.JacobianbyRef(x, output); + return output; + } + + public virtual void SsqGradientByRef(Vector x, Vector output) + { + if (output.Count != this.ParameterDimension) + throw new ArgumentException($"Output vector must match parameter dimension of function; expected {this.ParameterDimension}, got {output.Count}."); + + for (int jj = 0; jj < this.ParameterDimension; ++jj) + output[jj] = 0.0; + var tmp_grad = new LinearAlgebra.Double.DenseVector(this.ParameterDimension); + double tmp_value = 0.0; + + for (int ii = 0; ii < this.ItemDimension; ++ii) + { + tmp_value = this.ItemValue(x, ii); + this.ItemGradientByRef(x, ii, tmp_grad); + for (int jj = 0; jj < this.ParameterDimension; ++jj) + output[jj] += 2 * tmp_value * tmp_grad[jj]; + } + } + + public virtual Vector SsqGradient(Vector x) + { + var output = new LinearAlgebra.Double.DenseVector(this.ParameterDimension); + this.SsqGradientByRef(x, output); + return output; + } + + public virtual void SsqHessianByRef(Vector x, Matrix output) + { + if (output.RowCount != this.ParameterDimension || output.ColumnCount != this.ParameterDimension) + throw new ArgumentException($"Output matrix must match parameter dimension of function; expected {this.ParameterDimension}x{this.ParameterDimension}, got {output.RowCount}x{output.ColumnCount}."); + + for (int ii = 0; ii < this.ParameterDimension; ++ii) + for (int jj = 0; jj < this.ParameterDimension; ++jj) + output[ii,jj] = 0.0; + + var tmp_grad = new LinearAlgebra.Double.DenseVector(this.ParameterDimension); + var tmp_hess = new LinearAlgebra.Double.DenseMatrix(this.ParameterDimension, this.ParameterDimension); + double tmp_value = 0.0; + + for (int ii = 0; ii < this.ItemDimension; ++ii) + { + tmp_value = this.ItemValue(x, ii); + this.ItemGradientByRef(x, ii, tmp_grad); + this.ItemHessianByRef(x, ii, tmp_hess); + for (int jj = 0; jj < this.ParameterDimension; ++jj) + { + for (int kk = 0; kk < this.ParameterDimension; ++kk) + { + var increment = 2 * (tmp_value * tmp_hess[jj, kk] + tmp_grad[jj] * tmp_grad[kk]); + output[jj, kk] += increment; + } + } + } + } + + public virtual Matrix SsqHessian(Vector x) + { + var output = new LinearAlgebra.Double.DenseMatrix(this.ParameterDimension, this.ParameterDimension); + this.SsqHessianByRef(x, output); + return output; + } + + public virtual double SsqValue(Vector x) + { + double ssq = 0.0; + for (int ii = 0; ii < this.ItemDimension; ++ii) + { + var tmp = this.ItemValue(x, ii); + ssq += tmp * tmp; + } + return ssq; + } + + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/BealeFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/BealeFunction.cs new file mode 100644 index 00000000..b55208dd --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/BealeFunction.cs @@ -0,0 +1,97 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class BealeFunction : BaseTestFunction + { + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new BealeFunction(), + InitialGuess = new double[] { 1, 1 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 3, 0.5 }, + CaseName = "unbounded" + }; + yield return new TestCase() + { + Function = new BealeFunction(), + InitialGuess = new double[] { 1, 1 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 3, 0.5 }, + LowerBound = new double[] { -1000, -1000}, + UpperBound = new double[] { 1000, 1000}, + CaseName = "loose bounds" + }; + yield return new TestCase() + { + Function = new BealeFunction(), + InitialGuess = new double[] { 1, 1 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 3, 0.5 }, + LowerBound = new double[] { 0.6, 0.5 }, + UpperBound = new double[] { 10, 100 }, + CaseName = "tight bounds" + }; + } + } + + public BealeFunction() { } + + public override string Description + { + get + { + return "Beale fun (MGH #5)"; + } + } + + public override int ItemDimension + { + get + { + return 3; + } + } + + public override int ParameterDimension + { + get + { + return 2; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + int ii = itemIndex + 1; + output[0] = -1 + Math.Pow(x[1], ii); + output[1] = ii * x[0] * Math.Pow(x[1], ii - 1); + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + int ii = itemIndex + 1; + output[0, 0] = 0; + output[0, 1] = ii * Math.Pow(x[1], ii - 1); + output[1, 0] = ii * Math.Pow(x[1], ii - 1); + output[1, 1] = (ii - 1) * ii * x[0] * Math.Pow(x[1], ii - 2); + } + + private static readonly double[] y = { 1.5, 2.25, 2.625}; + + public override double ItemValue(Vector x, int itemIndex) + { + int ii = itemIndex + 1; + return y[itemIndex] - x[0] * (1 - Math.Pow(x[1], ii)); + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/BrownAndDennisFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/BrownAndDennisFunction.cs new file mode 100644 index 00000000..9e5b9c41 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/BrownAndDennisFunction.cs @@ -0,0 +1,114 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class BrownAndDennisFunction : BaseTestFunction + { + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new BrownAndDennisFunction(20), + InitialGuess = new double[] { 25, 5, -5, -1 }, + MinimalValue = 85822.2, + MinimizingPoint = null, + CaseName = "unbounded" + }; + yield return new TestCase() + { + Function = new BrownAndDennisFunction(20), + InitialGuess = new double[] { 25, 5, -5, -1 }, + MinimalValue = 85822.2, + MinimizingPoint = null, + LowerBound = new double[] { -1000, -1000, -1000, -1000 }, + UpperBound = new double[] {1000, 1000, 1000, 1000 }, + CaseName = "loose bounds" + }; + yield return new TestCase() + { + Function = new BrownAndDennisFunction(20), + InitialGuess = new double[] { 25, 5, -5, -1 }, + MinimalValue = 0.88860479e5, + MinimizingPoint = null, + LowerBound = new double[] { -10, 0, -100, -20 }, + UpperBound = new double[] { 100, 15, 0, 0.2 }, + CaseName = "tight bounds" + }; + } + } + + private readonly int _items; + + public BrownAndDennisFunction(int items) + { + if (items < 4) + throw new ArgumentException("items must be >= 4"); + _items = items; + } + + public override string Description + { + get + { + return "Brown & Dennis fun (MGH #16)"; + } + } + + public override int ItemDimension + { + get + { + return _items; + } + } + + public override int ParameterDimension + { + get + { + return 4; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + var ii = itemIndex + 1; + var t = ii / 5.0; + output[0] = 2 * (x[0] + t * x[1] - Math.Exp(t)); + output[1] = (2*ii/25.0) * (5 * x[0] + ii * x[1] - 5 * Math.Exp(t)); + output[2] = 2 * (x[2] + x[3] * Math.Sin(t) - Math.Cos(t)); + output[3] = 2 * Math.Sin(t) * (x[2] + Math.Sin(t) * x[3] - Math.Cos(t)); + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + for (int ii = 0; ii < 4; ++ii) + for (int jj = 0; jj < 4; ++jj) + output[ii, jj] = 0; + var i = itemIndex + 1; + var t = i / 5.0; + output[0, 0] = 2; + output[0, 1] = 2 * t; + output[1, 0] = 2 * t; + output[1, 1] = 2 * t * t; + output[2, 2] = 2; + output[2, 3] = 2 * Math.Sin(t); + output[3, 2] = 2 * Math.Sin(t); + output[3, 3] = 2 * Math.Pow(Math.Sin(t), 2); + } + + public override double ItemValue(Vector x, int itemIndex) + { + var ii = itemIndex + 1; + var t = ii / 5.0; + return Math.Pow(x[0] + t * x[1] - Math.Exp(t), 2.0) + Math.Pow(x[2] + x[3] * Math.Sin(t) - Math.Cos(t), 2); + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/BrownBadlyScaledFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/BrownBadlyScaledFunction.cs new file mode 100644 index 00000000..e0056244 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/BrownBadlyScaledFunction.cs @@ -0,0 +1,131 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class BrownBadlyScaledFunction : BaseTestFunction + { + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new BrownBadlyScaledFunction(), + InitialGuess = new double[] { 1, 1 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 1e6, 2e-6 }, + CaseName = "unbounded" + }; + yield return new TestCase() + { + Function = new BrownBadlyScaledFunction(), + InitialGuess = new double[] { 1, 1 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 1e6, 2e-6 }, + LowerBound = new double[] { -1e8, -1e8 }, + UpperBound = new double[] { 1e8, 1e8 }, + CaseName = "loose bounds" + }; + yield return new TestCase() + { + Function = new BrownBadlyScaledFunction(), + InitialGuess = new double[] { 1, 1 }, + MinimalValue = 0.784e3, + MinimizingPoint = new double[] { 1e6, 2e-6 }, + LowerBound = new double[] { 0, 3e-5 }, + UpperBound = new double[] { 1e6, 100 }, + CaseName = "tight bounds" + }; + } + } + + public BrownBadlyScaledFunction() { } + + public override string Description + { + get + { + return "Brown badly scaled fun (MGH #4)"; + } + } + + public override int ItemDimension + { + get + { + return 3; + } + } + + public override int ParameterDimension + { + get + { + return 2; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + switch (itemIndex) + { + case 0: + output[0] = 1; + output[1] = 0; + break; + case 1: + output[0] = 0; + output[1] = 1; + break; + case 2: + output[0] = x[1]; + output[1] = x[0]; + break; + default: + throw new ArgumentException("itemIndex must be <= 2"); + } + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + switch (itemIndex) + { + case 0: + case 1: + output[0, 0] = 0; + output[0, 1] = 0; + output[1, 0] = 0; + output[1, 1] = 0; + break; + case 2: + output[0, 0] = 0; + output[0, 1] = 1; + output[1, 0] = 1; + output[1, 1] = 0; + break; + default: + throw new ArgumentException("itemIndex must be <= 2"); + } + } + + public override double ItemValue(Vector x, int itemIndex) + { + switch (itemIndex) + { + case 0: + return x[0] - 1e6; + case 1: + return x[1] - 2e-6; + case 2: + return x[0] * x[1] - 2; + default: + throw new ArgumentException("itemIndex must be <= 2"); + } + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/FreudensteinAndRothFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/FreudensteinAndRothFunction.cs new file mode 100644 index 00000000..16cbe4c2 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/FreudensteinAndRothFunction.cs @@ -0,0 +1,98 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; +using DenseVector = MathNet.Numerics.LinearAlgebra.Double.DenseVector; +using DenseMatrix = MathNet.Numerics.LinearAlgebra.Double.DenseMatrix; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class FreudensteinAndRothFunction : BaseTestFunction + { + + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new FreudensteinAndRothFunction(), + InitialGuess = new double[] { 0.5, -2 }, + MinimizingPoint = new double[] { 5, 4 }, + MinimalValue = 0, + CaseName = "unbounded" + }; + yield return new TestCase() + { + Function = new FreudensteinAndRothFunction(), + InitialGuess = new double[] { 0.5, -2 }, + MinimizingPoint = new double[] {5, 4}, + MinimalValue = 0, + LowerBound = new double[] { -1000, -1000 }, + UpperBound = new double[] { 1000, 1000}, + CaseName = "loose bounds" + }; + } + } + + public override string Description { get { return "Freudenstein & Roth fun (MGH #2)"; } } + + public override int ParameterDimension + { + get + { + return 2; + } + } + + public override int ItemDimension + { + get + { + return 2; + } + } + + public override double ItemValue(Vector x, int itemIndex) + { + if (itemIndex == 0) + return -13 + x[0] + ((5 - x[1]) * x[1] - 2) * x[1]; + else + return -29 + x[0] + ((x[1] + 1) * x[1] - 14) * x[1]; + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + if (itemIndex == 0) + { + output[0] = 1; + output[1] = -2 + (5 - 2 * x[1]) * x[1] + (5 - x[1]) * x[1]; + } + else + { + output[0] = 1; + output[1] = -14 + x[1] * (1 + x[1]) + x[1] * (1 + 2 * x[1]); + } + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + if (itemIndex == 0) + { + output[0, 0] = 0; + output[0, 1] = 0; + output[1, 0] = 0; + output[1, 1] = 10 - 6 * x[1]; + } + else + { + output[0, 0] = 0; + output[0, 1] = 0; + output[1, 0] = 0; + output[1, 1] = 2 + 6 * x[1]; + } + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/HelicalValleyFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/HelicalValleyFunction.cs new file mode 100644 index 00000000..987eb7fe --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/HelicalValleyFunction.cs @@ -0,0 +1,178 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class HelicalValleyFunction : BaseTestFunction + { + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new HelicalValleyFunction(), + InitialGuess = new double[] { -1, 0, 0 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 1, 0, 0 }, + CaseName = "unbounded" + }; + yield return new TestCase() + { + Function = new HelicalValleyFunction(), + InitialGuess = new double[] { -1, 0, 0 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 1, 0, 0 }, + LowerBound = new double[] { -1000, -1000, -1000 }, + UpperBound = new double[] { 1000, 1000, 1000 }, + CaseName = "loose bounds" + }; + yield return new TestCase() + { + Function = new HelicalValleyFunction(), + InitialGuess = new double[] { -1, 0, 0 }, + MinimalValue = 0.99042212, + LowerBound = new double[] { -100, -1, -1 }, + UpperBound = new double[] { 0.8, 1, 1 }, + CaseName = "tight bounds" + }; + } + } + + public override string Description + { + get + { + return "Helical valley fun (MGH #7)"; + } + } + + public override int ItemDimension + { + get + { + return 3; + } + } + + public override int ParameterDimension + { + get + { + return 3; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + switch (itemIndex) + { + case 0: + output[0] = -100 * theta10(x[0], x[1]); + output[1] = -100 * theta01(x[0], x[1]); + output[2] = 10; + break; + case 1: + output[0] = (10 * x[0]) / Math.Sqrt(x[0]*x[0] + x[1]*x[1]); + output[1] = (10 * x[1]) / Math.Sqrt(x[0]*x[0] + x[1]*x[1]); + output[2] = 0; + break; + case 2: + output[0] = 0; + output[1] = 0; + output[2] = 1; + break; + default: + throw new ArgumentException("itemIndex must be <= 2"); + } + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + switch (itemIndex) + { + case 0: + output[0, 0] = -100 * theta20(x[0], x[1]); + output[0, 1] = -100 * theta11(x[0], x[1]); + output[0, 2] = 0; + output[1, 0] = -100 * theta11(x[0], x[1]); + output[1, 1] = -100 * theta02(x[0], x[1]); + output[1, 2] = 0; + output[2, 0] = 0; + output[2, 1] = 0; + output[2, 2] = 0; + break; + case 1: + output[0, 0] = (10 * x[1]*x[1]) / Math.Pow(x[0]*x[0] + x[1]*x[1],1.5); + output[0, 1] = (-10 * x[0] * x[1]) / Math.Pow(x[0]*x[0] + x[1]*x[1],1.5); + output[0, 2] = 0; + output[1, 0] = (-10 * x[0] * x[1]) / Math.Pow(x[0] * x[0] + x[1] * x[1], 1.5); + output[1, 1] = (10 * x[0]*x[0]) / Math.Pow(x[0] * x[0] + x[1] * x[1], 1.5); + output[1, 2] = 0; + output[2, 0] = 0; + output[2, 1] = 0; + output[2, 2] = 0; + break; + case 2: + for (int ii = 0; ii < 2; ++ii) + for (int jj = 0; jj < 2; ++jj) + output[ii, jj] = 0; + break; + default: + throw new ArgumentException("itemIndex must be <= 2"); + } + } + + private static double theta(double x1, double x2) + { + if (x1 >= 0) + return 0.5 * Math.Atan(x2 / x1) / Math.PI; + else + return 0.5 * Math.Atan(x2 / x1) / Math.PI + 0.5; + } + + private static double theta10(double x1, double x2) + { + return -(x2 / (2 * Math.PI * Math.Pow(x1,2) + 2 * Math.PI * Math.Pow(x2,2))); + } + + private static double theta01(double x1, double x2) + { + return x1 / (2 * Math.PI * x1*x1 + 2 * Math.PI * x2*x2); + } + + private static double theta20(double x1,double x2) + { + return (x1 * x2) / (Math.PI * Math.Pow(x1 * x1 + x2 * x2, 2)); + } + + private static double theta11(double x1, double x2) + { + return (-x1 * x1 + x2 * x2) / (2 * Math.PI * Math.Pow(x1 * x1 + x2 * x2, 2)); + } + + private static double theta02(double x1, double x2) + { + return -((x1 * x2) / (Math.PI * Math.Pow(x1*x1 + x2*x2, 2))); + } + + public override double ItemValue(Vector x, int itemIndex) + { + switch (itemIndex) + { + case 0: + return 10 * (x[2] - 10 * theta(x[0], x[1])); + case 1: + return 10 * (Math.Sqrt(x[0] * x[0] + x[1] * x[1]) - 1); + case 2: + return x[2]; + default: + throw new ArgumentException("itemIndex must be <= 2"); + } + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/ITestFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/ITestFunction.cs new file mode 100644 index 00000000..563faec2 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/ITestFunction.cs @@ -0,0 +1,69 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; +using DenseVector = MathNet.Numerics.LinearAlgebra.Double.DenseVector; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class TestCase + { + public string CaseName; + public ITestFunction Function; + public DenseVector InitialGuess; + public DenseVector LowerBound; + public DenseVector UpperBound; + public double MinimalValue; + public DenseVector MinimizingPoint; + + public bool IsBounded + { + get + { + return this.LowerBound != null && this.UpperBound != null; + } + } + + public bool IsUnbounded + { + get + { + return this.IsUnboundedOverride ?? this.LowerBound == null || this.UpperBound == null; + } + } + + public bool? IsUnboundedOverride; + + public string FullName + { + get + { + return $"{this.Function.Description} {this.CaseName}"; + } + } + } + + public interface ITestFunction + { + string Description { get; } + int ParameterDimension { get; } + int ItemDimension { get; } + + double ItemValue(Vector x, int itemIndex); + Vector ItemGradient(Vector x, int itemIndex); + void ItemGradientByRef(Vector x, int itemIndex, Vector output); + Matrix ItemHessian(Vector x, int itemIndex); + void ItemHessianByRef(Vector x, int itemIndex, Matrix output); + + Matrix Jacobian(Vector x); + void JacobianbyRef(Vector x, Matrix output); + + double SsqValue(Vector x); + Vector SsqGradient(Vector x); + void SsqGradientByRef(Vector x, Vector output); + Matrix SsqHessian(Vector x); + void SsqHessianByRef(Vector x, Matrix output); + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/JennrichAndSampsonFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/JennrichAndSampsonFunction.cs new file mode 100644 index 00000000..58cc30eb --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/JennrichAndSampsonFunction.cs @@ -0,0 +1,102 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class JennrichAndSampsonFunction : BaseTestFunction + { + private readonly int _m; + + public JennrichAndSampsonFunction(int itemDimension) + { + if (itemDimension < 2) + throw new ArgumentException("itemDimension must be at least 2."); + _m = itemDimension; + } + + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new JennrichAndSampsonFunction(10), + InitialGuess = new double[] { 0.3, 0.4 }, + MinimalValue = 124.362, + MinimizingPoint = new double[] { 0.2578, 0.2578 }, + CaseName = "unbounded" + }; + //yield return new TestCase() + //{ + // Function = new JennrichAndSampsonFunction(10), + // LowerBound = new double[] { 0.6, 0.5 }, + // UpperBound = new double[] { 10, 50 }, + // StartPoint = new double[] { 1.0, 1.0 }, + // MinimizingInput = null, + // MinimizingValue = 0, + // CaseName = "tight bounds" + //}; + yield return new TestCase() + { + Function = new JennrichAndSampsonFunction(10), + LowerBound = new double[] { -50, -50 }, + UpperBound = new double[] { 50, 50 }, + InitialGuess = new double[] { 0.3, 0.4 }, + MinimizingPoint = null, + MinimalValue = 0, + CaseName = "loose bounds" + }; + } + } + + public override string Description + { + get + { + return $"Jennrich & Sampson fun (MGH #6) (n={this.ItemDimension})"; + } + } + + public override int ItemDimension + { + get + { + return _m; + } + } + + public override int ParameterDimension + { + get + { + return 2; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + int ii = itemIndex + 1; + output[0] = -(Math.Exp(ii * x[0]) * ii); + output[1] = -(Math.Exp(ii * x[1]) * ii); + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + int ii = itemIndex + 1; + output[0, 0] = -(Math.Exp(ii * x[0]) * ii*ii); + output[0, 1] = 0; + output[1, 0] = 0; + output[1, 1] = -(Math.Exp(ii * x[1]) * ii*ii); + } + + public override double ItemValue(Vector x, int itemIndex) + { + int ii = itemIndex + 1; + return 2 + 2 * ii - (Math.Exp(ii * x[0]) + Math.Exp(ii * x[1])); + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/MeyerFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/MeyerFunction.cs new file mode 100644 index 00000000..961da226 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/MeyerFunction.cs @@ -0,0 +1,99 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class MeyerFunction : BaseTestFunction + { + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new MeyerFunction(), + InitialGuess = new double[] { 0.02, 4000, 250 }, + MinimalValue = 87.9458, + MinimizingPoint = null, + CaseName = "unbounded" + }; + yield return new TestCase() + { + Function = new MeyerFunction(), + InitialGuess = new double[] { 0.02, 4000, 250 }, + MinimalValue = 87.9458, + MinimizingPoint = null, + LowerBound = new double[] { -1e6, -1e6, -1e6 }, + UpperBound = new double[] { 1e6, 1e6, 1e6 }, + CaseName = "loose bounds" + }; + } + } + + public MeyerFunction() { } + + public override string Description + { + get + { + return "Meyer fun (MGH #10)"; + } + } + + public override int ItemDimension + { + get + { + return 16; + } + } + + public override int ParameterDimension + { + get + { + return 3; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + int ii = itemIndex + 1; + output[0] = Math.Exp(x[1] / (45.0 + 5 * ii + x[2])); + output[1] = (Math.Exp(x[1] / (45.0 + 5 * ii + x[2])) * x[0]) / (45 + 5 * ii + x[2]); + output[2] = -(Math.Exp(x[1] / (45.0 + 5 * ii + x[2])) * x[0] * x[1]) / Math.Pow(45 + 5 * ii + x[2], 2); + + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + var ii = itemIndex + 1; + + var t0 = (45.0 + 5 * ii + x[2]); + var t1 = Math.Exp(x[1] / t0); + + output[0, 0] = 0; + output[0, 1] = t1 / t0; + output[0, 2] = -t1 * x[1] / Math.Pow(t0, 2); + output[1, 0] = t1 / t0; + output[1, 1] = t1 * x[0] / Math.Pow(t0, 2); + output[1, 2] = -t1 * x[0] * (t0 + x[1]) / Math.Pow(t0, 3); + output[2, 0] = -t1 * x[1] / Math.Pow(t0, 2); + output[2, 1] = -t1 * x[0] * (t0 + x[1]) / Math.Pow(t0, 3); + output[2, 2] = t1 * x[0] * x[1] * (2*t0 + x[1]) / Math.Pow(t0, 4); + } + + private static readonly double[] y = { 34780, 28610, 23650, 19630, 16370, 13720, 11540, 9744, 8261, 7030, 6005, 5147, 4427, 3820, 3307, 2872 }; + + public override double ItemValue(Vector x, int itemIndex) + { + var ii = itemIndex + 1; + var t = 45.0 + 5 * ii; + return x[0] * Math.Exp(x[1] / (t + x[2])) - y[itemIndex]; + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/PowellBadlyScaledFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/PowellBadlyScaledFunction.cs new file mode 100644 index 00000000..8885a39d --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/PowellBadlyScaledFunction.cs @@ -0,0 +1,111 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.LinearAlgebra.Double; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class PowellBadlyScaledFunction : BaseTestFunction + { + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new PowellBadlyScaledFunction(), + InitialGuess = new double[] { 0, 1 }, + MinimizingPoint = new double[] { 1.098e-5, 9.106 }, + MinimalValue = 0, + CaseName = "unbounded" + }; + yield return new TestCase() + { + Function = new PowellBadlyScaledFunction(), + InitialGuess = new double[] { 0, 1 }, + MinimizingPoint = new double[] { 1.098e-5, 9.106 }, + MinimalValue = 0, + LowerBound = new double[] { -1000, -1000 }, + UpperBound = new double[] { 1000, 1000 }, + CaseName = "loose bounds" + }; + yield return new TestCase() + { + Function = new PowellBadlyScaledFunction(), + LowerBound = new double[] { 0, 1 }, + UpperBound = new double[] { 1, 9 }, + InitialGuess = new double[] { 0, 1 }, + MinimalValue = 0.15125900e-9, + CaseName = "tight bounds" + }; + } + } + + public override string Description + { + get + { + return "Powell badly scaled fun (MGH #3)"; + } + } + + public override int ItemDimension + { + get + { + return 2; + } + } + + public override int ParameterDimension + { + get + { + return 2; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + if (itemIndex == 0) + { + output[0] = 10000 * x[1]; + output[1] = 10000 * x[0]; + } + else if (itemIndex == 1) + { + output[0] = -Math.Exp(-x[0]); + output[1] = -Math.Exp(-x[1]); + } + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + if (itemIndex == 0) + { + output[0, 0] = 0; + output[0, 1] = 10000; + output[1, 0] = 10000; + output[1, 1] = 0; + } + else + { + output[0, 0] = Math.Exp(-x[0]); + output[0, 1] = 0; + output[1, 0] = 0; + output[1,1] = Math.Exp(-x[1]); + } + } + + public override double ItemValue(Vector x, int itemIndex) + { + if (itemIndex == 0) + return 10000.0 * x[0] * x[1] - 1; + else + return Math.Exp(-x[0]) + Math.Exp(-x[1]) - 1.0001; + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/PowellSingularFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/PowellSingularFunction.cs new file mode 100644 index 00000000..130b057a --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/PowellSingularFunction.cs @@ -0,0 +1,141 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class PowellSingularFunction : BaseTestFunction + { + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new PowellSingularFunction(), + InitialGuess = new double[] { 3, -1, 0, 1 }, + MinimalValue = 0, + MinimizingPoint = new double[] {0,0,0,0}, + CaseName = "unbounded" + }; + yield return new TestCase() + { + Function = new PowellSingularFunction(), + InitialGuess = new double[] { 3, -1, 0, 1 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 0, 0, 0, 0 }, + LowerBound = new double[] {-1000, -1000, -1000, -1000}, + UpperBound = new double[] { 1000, 1000, 1000, 1000 }, + CaseName = "loose bounds" + }; + } + } + + + public PowellSingularFunction() { } + + public override string Description + { + get + { + return "Powell singular fun (MGH #13)"; + } + } + + public override int ItemDimension + { + get + { + return 4; + } + } + + public override int ParameterDimension + { + get + { + return 4; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + switch (itemIndex) + { + case 0: + output[0] = 1; + output[1] = 10; + output[2] = 0; + output[3] = 0; + break; + case 1: + output[0] = 0; + output[1] = 0; + output[2] = Math.Sqrt(5); + output[3] = -Math.Sqrt(5); + break; + case 2: + output[0] = 0; + output[1] = 2*(x[1]-2*x[2]); + output[2] = -4*x[1] + 8*x[2]; + output[3] = 0; + break; + case 3: + output[0] = 2*Math.Sqrt(10)*(x[0] - x[3]); + output[1] = 0; + output[2] = 0; + output[3] = -2*Math.Sqrt(10)*(x[0] - x[3]); + break; + default: + throw new ArgumentException("itemIndex must be <= 3"); + } + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + for (int ii = 0; ii < 4; ++ii) + for (int jj = 0; jj < 4; ++jj) + output[ii, jj] = 0; + switch(itemIndex) + { + case 0: + case 1: + break; + case 2: + output[1, 1] = 2; + output[1, 2] = -4; + output[2, 1] = -4; + output[2, 2] = 8; + break; + case 3: + output[0, 0] = 2 * Math.Sqrt(10); + output[0, 3] = -2 * Math.Sqrt(10); + output[3, 0] = -2 * Math.Sqrt(10); + output[3, 3] = 2 * Math.Sqrt(10); + break; + default: + throw new ArgumentException("itemIndex must be <= 3"); + } + } + + public override double ItemValue(Vector x, int itemIndex) + { + switch (itemIndex) + { + case 0: + return x[0] + 10 * x[1]; + case 1: + return Math.Sqrt(5) * (x[2] - x[3]); + case 2: + return Math.Pow(x[1] - 2 * x[2], 2); + case 3: + return Math.Sqrt(10.0) * Math.Pow(x[0] - x[3], 2); + default: + throw new ArgumentException("itemIndex must be <= 3"); + } + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/RosenbrockFunction2.cs b/src/UnitTests/OptimizationTests/TestFunctions/RosenbrockFunction2.cs new file mode 100644 index 00000000..c693b507 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/RosenbrockFunction2.cs @@ -0,0 +1,156 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class RosenbrockFunction2 : BaseTestFunction + { + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new RosenbrockFunction2(), + InitialGuess = new double[] { -1.2, 1 }, + MinimizingPoint = new double[] { 1, 1 }, + MinimalValue = 0, + LowerBound = new double[] { -1000, -1000 }, + UpperBound = new double[] { 1000, 1000 }, + CaseName = "hard start", + IsUnboundedOverride = true + }; + yield return new TestCase() + { + Function = new RosenbrockFunction2(), + InitialGuess = new double[] { 1.2, 1.2 }, + MinimizingPoint = new double[] { 1, 1 }, + MinimalValue = 0, + LowerBound = new double[] { -5, -5 }, + UpperBound = new double[] { 5, 5 }, + CaseName = "easy start" + }; + yield return new TestCase() + { + Function = new RosenbrockFunction2(), + InitialGuess = new double[] { -0.9, -0.5 }, + MinimizingPoint = new double[] { 1, 1 }, + MinimalValue = 0, + LowerBound = new double[] { -5, -5 }, + UpperBound = new double[] { 5, 5 }, + CaseName = "Overton start", + IsUnboundedOverride = true + }; + yield return new TestCase() + { + Function = new RosenbrockFunction2(), + InitialGuess = new double[] { 1.2, 1.2 }, + MinimizingPoint = new double[] { 1, 1 }, + MinimalValue = 0, + LowerBound = new double[] { 1, -5 }, + UpperBound = new double[] { 5, 5 }, + CaseName = "easy one active bound" + }; + yield return new TestCase() + { + Function = new RosenbrockFunction2(), + InitialGuess = new double[] { 1.2, 1.2 }, + MinimizingPoint = new double[] { 1, 1 }, + MinimalValue = 0, + LowerBound = new double[] { 1, 1 }, + UpperBound = new double[] { 5, 5 }, + CaseName = "easy two active bounds" + }; + yield return new TestCase() + { + Function = new RosenbrockFunction2(), + InitialGuess = new double[] { 2.5, 2.5 }, + MinimizingPoint = new double[] { 2, 4 }, + MinimalValue = 1, + LowerBound = new double[] { 2, 2 }, + UpperBound = new double[] { 5, 5 }, + CaseName = "min on lower bound, not local" + }; + yield return new TestCase() + { + Function = new RosenbrockFunction2(), + InitialGuess = new double[] { -0.9, -0.5 }, + MinimizingPoint = new double[] { 0.5, 0.25 }, + MinimalValue = 0.25, + LowerBound = new double[] { -2, -2 }, + UpperBound = new double[] { 0.5, 0.5 }, + CaseName = "min on upper bound, not local" + }; + + + } + } + public RosenbrockFunction2() { } + + public override string Description + { + get + { + return "Rosenbrock fun (MGH #1)"; + } + } + + public override int ItemDimension + { + get + { + return 2; + } + } + + public override int ParameterDimension + { + get + { + return 2; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + if (itemIndex == 0) + { + output[0] = -20 * x[0]; + output[1] = 10; + } else + { + output[0] = -1; + output[1] = 0; + } + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + if (itemIndex == 0) + { + output[0, 0] = -20; + output[0, 1] = 0; + output[1, 0] = 0; + output[1, 1] = 0; + } else + { + output[0, 0] = 0; + output[0, 1] = 0; + output[1, 0] = 0; + output[1, 1] = 0; + } + } + + public override double ItemValue(Vector x, int itemIndex) + { + if (itemIndex == 0) + return 10 * (x[1] - x[0] * x[0]); + else + return 1 - x[0]; + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestFunctions/WoodFunction.cs b/src/UnitTests/OptimizationTests/TestFunctions/WoodFunction.cs new file mode 100644 index 00000000..4e8f014e --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestFunctions/WoodFunction.cs @@ -0,0 +1,164 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions +{ + public class WoodFunction : BaseTestFunction + { + public static IEnumerable TestCases + { + get + { + yield return new TestCase() + { + Function = new WoodFunction(), + InitialGuess = new double[] { -3, -1, -3, -1 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 1, 1, 1, 1 }, + CaseName = "unbounded" + }; + yield return new TestCase() + { + Function = new WoodFunction(), + InitialGuess = new double[] { -3, -1, -3, -1 }, + MinimalValue = 0, + MinimizingPoint = new double[] { 1, 1, 1, 1 }, + LowerBound = new double[] { -1000, -1000, -1000, -1000 }, + UpperBound = new double[] { 1000, 1000, 1000, 1000 }, + CaseName = "loose bounds" + }; + yield return new TestCase() + { + Function = new WoodFunction(), + InitialGuess = new double[] { -3, -1, -3, -1 }, + MinimalValue = 1.5567008, + MinimizingPoint = null, + LowerBound = new double[] { -100, -100, -100, -100 }, + UpperBound = new double[] { 0, 10, 100, 100 }, + CaseName = "tight bounds" + }; + } + } + + public WoodFunction() { } + + public override string Description + { + get + { + return "Wood fun (MGH #14)"; + } + } + + public override int ItemDimension + { + get + { + return 6; + } + } + + public override int ParameterDimension + { + get + { + return 4; + } + } + + public override void ItemGradientByRef(Vector x, int itemIndex, Vector output) + { + switch (itemIndex) + { + case 0: + output[0] = -20 * x[0]; + output[1] = 10; + output[2] = 0; + output[3] = 0; + break; + case 1: + output[0] = -1; + output[1] = 0; + output[2] = 0; + output[3] = 0; + break; + case 2: + output[0] = 0; + output[1] = 0; + output[2] = -6 * Math.Sqrt(10) * x[2]; + output[3] = 3 * Math.Sqrt(10); + break; + case 3: + output[0] = 0; + output[1] = 0; + output[2] = -1; + output[3] = 0; + break; + case 4: + output[0] = 0; + output[1] = Math.Sqrt(10); + output[2] = 0; + output[3] = Math.Sqrt(10); + break; + case 5: + output[0] = 0; + output[1] = 1.0 / Math.Sqrt(10); + output[2] = 0; + output[3] = -1.0 / Math.Sqrt(10); + break; + default: + throw new ArgumentException("itemIndex must be <= 5"); + } + } + + public override void ItemHessianByRef(Vector x, int itemIndex, Matrix output) + { + for (int ii = 0; ii < 4; ++ii) + for (int jj = 0; jj < 4; ++jj) + output[ii, jj] = 0; + switch (itemIndex) + { + case 0: + output[0, 0] = -20; + break; + case 1: + break; + case 2: + output[2, 2] = -6 * Math.Sqrt(10); + break; + case 3: + case 4: + case 5: + break; + default: + throw new ArgumentException("itemIndex must be <= 5"); + + } + } + + public override double ItemValue(Vector x, int itemIndex) + { + switch (itemIndex) + { + case 0: + return 10 * (x[1] - x[0] * x[0]); + case 1: + return 1 - x[0]; + case 2: + return Math.Sqrt(90) * (x[3] - x[2] * x[2]); + case 3: + return 1 - x[2]; + case 4: + return Math.Sqrt(10) * (x[1] + x[3] - 2); + case 5: + return (x[1] - x[3]) / Math.Sqrt(10); + default: + throw new ArgumentException("itemIndex must be <= 5"); + } + } + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 25595c6b..a4c28930 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -30,7 +30,7 @@ AllRules.ruleset 1591 false - 5 + 6 DEBUG;TRACE @@ -46,7 +46,7 @@ AnyCPU 1591 false - 5 + 6 ..\..\out\test-signed\Net40\ @@ -62,7 +62,7 @@ AllRules.ruleset 1591 false - 5 + 6 @@ -366,8 +366,32 @@ + + + + + + + + + + + + + + + + + + + + + + + + @@ -416,10 +440,16 @@ + + + + + +