Optimization: Add interval expansion to GoldenSectionMinimizer

src/Numerics/Optimization/GoldenSectionMinimizer.cs

@@ -6,16 +6,25 @@ namespace MathNet.Numerics.Optimization
     {
         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)
+        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)
         {
-            double middlePointX = lowerBound + (upperBound - lowerBound) / (1 + Constants.GoldenRatio);
+            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);
@@ -24,8 +33,28 @@ namespace MathNet.Numerics.Optimization
             ValueChecker(middle.Value, middlePointX);
             ValueChecker(upper.Value, upperBound);
 
-            if (upperBound <= lowerBound)
-                throw new OptimizationException("Lower bound must be lower than upper bound.");
+            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);
+                lower = objective.Evaluate(lowerBound);
+                middle = objective.Evaluate(middlePointX);
+                upper = objective.Evaluate(upperBound);
+
+                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.");
@@ -71,7 +100,7 @@ namespace MathNet.Numerics.Optimization
             return new MinimizationResult1D(middle, iterations, MinimizationResult.ExitCondition.BoundTolerance);
         }
 
-        private void ValueChecker(double value, double point)
+        void ValueChecker(double value, double point)
         {
             if (Double.IsNaN(value) || Double.IsInfinity(value))
                 throw new Exception("Objective function returned non-finite value.");

src/UnitTests/OptimizationTests/TestGoldenSectionMinimizer.cs

@@ -22,11 +22,11 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
         public void Test_ExpansionWorks()
         {
             var algorithm = new GoldenSectionMinimizer(1e-5, 1000);
-            var f1 = new Func<double, double>(x => (x - 3) * (x - 3));
+            var f1 = new Func<double, double>(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));
         }
     }
-}
+}