Brent Minimizer

src/Numerics.Tests/OptimizationTests/BrentMinimizerTests.cs

@@ -0,0 +1,60 @@
+// <copyright file="BrentMinimizerTests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// Copyright (c) 2009-2017 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.
+// </copyright>
+
+using System;
+using MathNet.Numerics.Optimization;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+    [TestFixture]
+    public class BrentMinimizerTests
+    {
+        [Test]
+        public void Test_Works()
+        {
+            var f1 = new Func<double, double>(x => (x - 3) * (x - 3));
+            var obj = ObjectiveFunction.ScalarValue(f1);
+            var r1 = BrentMinimizer.Minimum(obj, -100, 100);
+
+            Assert.That(Math.Abs(r1.MinimizingPoint - 3.0), Is.LessThan(1e-4));
+        }
+
+        [Test]
+        public void Test_ExpansionWorks()
+        {
+            var algorithm = new BrentMinimizer(1e-5, 1000);
+            var f1 = new Func<double, double>(x => (x - 3) * (x - 3));
+            var obj = ObjectiveFunction.ScalarValue(f1);
+            var r1 = algorithm.FindMinimum(obj, -2, 2);
+
+            Assert.That(Math.Abs(r1.MinimizingPoint - 3.0), Is.LessThan(1e-4));
+        }
+    }
+}

src/Numerics.Tests/OptimizationTests/GoldenSectionMinimizerTests.cs

@@ -53,7 +53,7 @@ public void Test_ExpansionWorks()
             var algorithm = new GoldenSectionMinimizer(1e-5, 1000);
             var f1 = new Func<double, double>(x => (x - 3)*(x - 3));
             var obj = ObjectiveFunction.ScalarValue(f1);
-            var r1 = algorithm.FindMinimum(obj, -5, 5);
+            var r1 = algorithm.FindMinimum(obj, -1, 1);
 
             Assert.That(Math.Abs(r1.MinimizingPoint - 3.0), Is.LessThan(1e-4));
         }

src/Numerics/Optimization/BrentMinimizer.cs

@@ -0,0 +1,193 @@
+// <copyright file="BrentMinimizer.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// Copyright (c) 2009-2017 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.
+// </copyright>
+
+using MathNet.Numerics.Optimization.ObjectiveFunctions;
+using System;
+
+namespace MathNet.Numerics.Optimization
+{
+    public class BrentMinimizer
+    {
+        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 BrentMinimizer(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 ScalarMinimizationResult FindMinimum(IScalarObjectiveFunction objective, double lowerBound, double upperBound)
+        {
+            return Minimum(objective, lowerBound, upperBound, XTolerance, MaximumIterations, MaximumExpansionSteps, LowerExpansionFactor, UpperExpansionFactor);
+        }
+
+        public static ScalarMinimizationResult Minimum(IScalarObjectiveFunction objective, double lowerBound, double upperBound, double xTolerance = 1e-5,
+            int maxIterations = 1000, int maxExpansionSteps = 10, double lowerExpansionFactor = 2.0, double upperExpansionFactor = 2.0)
+        {
+            int maxfun = maxIterations;
+
+            if (lowerBound > upperBound)
+                throw new OptimizationException("Lower bound must be lower than upper bound.");
+
+            double sqrt_eps = Math.Sqrt(2.2e-16);
+
+            // This is not the golden_mean, but golden angle. Not sure why.
+            //  https://en.wikipedia.org/wiki/Golden_angle
+            double golden_angle = 0.5 * (3.0 - Math.Sqrt(5.0));
+
+            double a = lowerBound;
+            double b = upperBound;
+            double fulc = a + golden_angle * (b - a);
+
+            double nfc = fulc, xf = fulc;
+            double rat = 0.0, e = 0.0;
+            double x = xf;
+            var evaluation = objective.Evaluate(x);
+            double fx = evaluation.Value;
+            int num = 1;
+
+            double fu = double.PositiveInfinity;
+            double ffulc = fx, fnfc = fx;
+            double xm = 0.5 * (a + b);
+            double tol1 = sqrt_eps * Math.Abs(xf) + xTolerance / 3.0;
+            double tol2 = 2.0 * tol1;
+
+            while (Math.Abs(xf - xm) > (tol2 - 0.5 * (b - a)))
+            {
+                bool golden = true;
+
+                // Check for parabolic fit
+                if (Math.Abs(e) > tol1)
+                {
+                    golden = false;
+                    double r = (xf - nfc) * (fx - ffulc);
+                    double q = (xf - fulc) * (fx - fnfc);
+                    double p = (xf - fulc) * q - (xf - nfc) * r;
+                    q = 2.0 * (q - r);
+                    if (q > 0.0)
+                        p = -p;
+                    q = Math.Abs(q);
+                    r = e;
+                    e = rat;
+
+                    // Check for acceptability of parabola
+                    if ((Math.Abs(p) < Math.Abs(0.5 * q * r)) && (p > q * (a - xf)) && (p < q * (b - xf)))
+                    {
+                        rat = (p + 0.0) / q;
+                        x = xf + rat;
+
+                        if (((x - a) < tol2) || ((b - x) < tol2))
+                        {
+                            int si_2 = Math.Sign(xm - xf) + ((xm - xf) == 0 ? 1 : 0);
+                            rat = tol1 * si_2;
+                        }
+                    }
+                    else     // do a golden-section step
+                        golden = true;
+                }
+
+                if (golden)  // do a golden-section step
+                {
+                    if (xf >= xm)
+                        e = a - xf;
+                    else
+                        e = b - xf;
+                    rat = golden_angle * e;
+                }
+
+                int si = Math.Sign(rat) + (rat == 0 ? 1 : 0);
+                x = xf + si * Math.Max(Math.Abs(rat), tol1);
+
+                evaluation = objective.Evaluate(x);
+                fu = evaluation.Value;
+                num += 1;
+
+                if (fu <= fx)
+                {
+                    if (x >= xf)
+                        a = xf;
+                    else
+                        b = xf;
+
+                    fulc = nfc; ffulc = fnfc;
+                    nfc = xf; fnfc = fx;
+                    xf = x; fx = fu;
+                }
+                else
+                {
+                    if (x < xf)
+                        a = x;
+                    else
+                        b = x;
+
+                    if ((fu <= fnfc) || (nfc == xf))
+                    {
+                        fulc = nfc; ffulc = fnfc;
+                        nfc = x; fnfc = fu;
+                    }
+                    else if ((fu <= ffulc) || (fulc == xf) || (fulc == nfc))
+                    {
+                        fulc = x; ffulc = fu;
+                    }
+                }
+
+                xm = 0.5 * (a + b);
+                tol1 = sqrt_eps * Math.Abs(xf) + xTolerance / 3.0;
+                tol2 = 2.0 * tol1;
+
+                if (num >= maxfun)
+                    break;
+            }
+
+            var exitCondition = ExitCondition.BoundTolerance;
+
+            if (num >= maxfun)
+                exitCondition = ExitCondition.ExceedIterations;
+            else if (double.IsNaN(xf) || double.IsNaN(fx) || double.IsNaN(fu))
+                exitCondition = ExitCondition.InvalidValues;
+
+            return new ScalarMinimizationResult(new ScalarValueObjectiveFunctionEvaluation(xf, fx), num, exitCondition);
+        }
+
+        static void ValueChecker(double value, double point)
+        {
+            if (Double.IsNaN(value) || Double.IsInfinity(value))
+            {
+                throw new Exception("Objective function returned non-finite value.");
+            }
+        }
+    }
+}