Optimization: clarify objective function forking and mutation, fixes unit tests

src/Numerics/Optimization/IObjectiveFunction.cs

@@ -2,9 +2,8 @@
 
 namespace MathNet.Numerics.Optimization
 {
-    public interface IObjectiveFunction
+    public interface IObjectiveFunctionEvaluation
     {
-        void EvaluateAt(Vector<double> point);
         IObjectiveFunction Fork();
 
         Vector<double> Point { get; }
@@ -16,4 +15,9 @@ namespace MathNet.Numerics.Optimization
         bool IsHessianSupported { get; }
         Matrix<double> Hessian { get; }
     }
+
+    public interface IObjectiveFunction : IObjectiveFunctionEvaluation
+    {
+        void EvaluateAt(Vector<double> point);
+    }
 }

src/Numerics/Optimization/Implementation/WeakWolfeLineSearch.cs

@@ -19,8 +19,9 @@ namespace MathNet.Numerics.Optimization.Implementation
         }
 
         // Implemented following http://www.math.washington.edu/~burke/crs/408/lectures/L9-weak-Wolfe.pdf
-        public LineSearchOutput FindConformingStep(IObjectiveFunction objective, IObjectiveFunction startingPoint, Vector<double> searchDirection, double initialStep)
+        public LineSearchOutput FindConformingStep(IObjectiveFunctionEvaluation startingPoint, Vector<double> searchDirection, double initialStep)
         {
+            var objective = startingPoint.Fork();
             if (!(objective is CheckedObjectiveFunction))
             {
                 objective = new CheckedObjectiveFunction(objective, ValidateValue, ValidateGradient, null);
@@ -30,21 +31,21 @@ namespace MathNet.Numerics.Optimization.Implementation
             double upperBound = Double.PositiveInfinity;
             double step = initialStep;
 
+            Vector<double> initialPoint = startingPoint.Point;
             double initialValue = startingPoint.Value;
             Vector<double> initialGradient = startingPoint.Gradient;
 
             double initialDd = searchDirection * initialGradient;
 
             int ii;
-            IObjectiveFunction candidateEval = objective.Fork();
             MinimizationOutput.ExitCondition reasonForExit = MinimizationOutput.ExitCondition.None;
             for (ii = 0; ii < _maximumIterations; ++ii)
             {
-                candidateEval.EvaluateAt(startingPoint.Point + searchDirection * step);
+                objective.EvaluateAt(initialPoint + searchDirection * step);
 
-                double stepDd = searchDirection * candidateEval.Gradient;
+                double stepDd = searchDirection * objective.Gradient;
 
-                if (candidateEval.Value > initialValue + _c1 * step * initialDd)
+                if (objective.Value > initialValue + _c1 * step * initialDd)
                 {
                     upperBound = step;
                     step = 0.5 * (lowerBound + upperBound);
@@ -63,9 +64,9 @@ namespace MathNet.Numerics.Optimization.Implementation
                 if (!Double.IsInfinity(upperBound))
                 {
                     double maxRelChange = 0.0;
-                    for (int jj = 0; jj < candidateEval.Point.Count; ++jj)
+                    for (int jj = 0; jj < objective.Point.Count; ++jj)
                     {
-                        double tmp = Math.Abs(searchDirection[jj] * (upperBound - lowerBound)) / Math.Max(Math.Abs(candidateEval.Point[jj]), 1.0);
+                        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)
@@ -86,7 +87,7 @@ namespace MathNet.Numerics.Optimization.Implementation
                 throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", _maximumIterations));
             }
 
-            return new LineSearchOutput(candidateEval, ii, step, reasonForExit);
+            return new LineSearchOutput(objective, ii, step, reasonForExit);
         }
 
         bool Conforms(IObjectiveFunction startingPoint, Vector<double> searchDirection, double step, IObjectiveFunction endingPoint)

src/Numerics/Optimization/NewtonMinimizer.cs

@@ -21,38 +21,41 @@ namespace MathNet.Numerics.Optimization
         public MinimizationOutput FindMinimum(IObjectiveFunction objective, Vector<double> 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.");
+            }
 
             if (!(objective is CheckedObjectiveFunction))
+            {
                 objective = new CheckedObjectiveFunction(objective, ValidateObjective, ValidateGradient, ValidateHessian);
-
+            }
-            IObjectiveFunction initialEval = objective.Fork();
-            initialEval.EvaluateAt(initialGuess);
 
             // Check that we're not already done
-            if (ExitCriteriaSatisfied(initialGuess, initialEval.Gradient))
+            objective.EvaluateAt(initialGuess);
-                return new MinimizationOutput(initialEval, 0, MinimizationOutput.ExitCondition.AbsoluteGradient);
+            if (ExitCriteriaSatisfied(objective.Gradient))
+            {
+                return new MinimizationOutput(objective, 0, MinimizationOutput.ExitCondition.AbsoluteGradient);
+            }
 
             // Set up line search algorithm
             var lineSearcher = new WeakWolfeLineSearch(1e-4, 0.9, 1e-4, maxIterations: 1000);
 
-            // Declare state variables
-            IObjectiveFunction candidatePoint = initialEval;
-
             // Subsequent steps
             int iterations = 0;
             int totalLineSearchSteps = 0;
             int iterationsWithNontrivialLineSearch = 0;
             bool tmpLineSearch = false;
-            while (!ExitCriteriaSatisfied(candidatePoint.Point, candidatePoint.Gradient) && iterations < MaximumIterations)
+            while (!ExitCriteriaSatisfied(objective.Gradient) && iterations < MaximumIterations)
             {
-                var searchDirection = candidatePoint.Hessian.LU().Solve(-candidatePoint.Gradient);
+                var searchDirection = objective.Hessian.LU().Solve(-objective.Gradient);
-                if (searchDirection * candidatePoint.Gradient >= 0)
+                if (searchDirection * objective.Gradient >= 0)
                 {
-                    searchDirection = -candidatePoint.Gradient;
+                    searchDirection = -objective.Gradient;
                     tmpLineSearch = true;
                 }
 
@@ -61,7 +64,7 @@ namespace MathNet.Numerics.Optimization
                     LineSearchOutput result;
                     try
                     {
-                        result = lineSearcher.FindConformingStep(objective, candidatePoint, searchDirection, 1.0);
+                        result = lineSearcher.FindConformingStep(objective.Fork(), searchDirection, 1.0);
                     }
                     catch (Exception e)
                     {
@@ -69,11 +72,11 @@ namespace MathNet.Numerics.Optimization
                     }
                     iterationsWithNontrivialLineSearch += result.Iterations > 0 ? 1 : 0;
                     totalLineSearchSteps += result.Iterations;
-                    candidatePoint = result.FunctionInfoAtMinimum;
+                    objective = result.FunctionInfoAtMinimum;
                 }
                 else
                 {
-                    candidatePoint.EvaluateAt(candidatePoint.Point + searchDirection);
+                    objective.EvaluateAt(objective.Point + searchDirection);
                 }
 
                 tmpLineSearch = false;
@@ -82,12 +85,14 @@ namespace MathNet.Numerics.Optimization
             }
 
             if (iterations == MaximumIterations)
+            {
                 throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations));
+            }
 
-            return new MinimizationWithLineSearchOutput(candidatePoint, iterations, MinimizationOutput.ExitCondition.AbsoluteGradient, totalLineSearchSteps, iterationsWithNontrivialLineSearch);
+            return new MinimizationWithLineSearchOutput(objective, iterations, MinimizationOutput.ExitCondition.AbsoluteGradient, totalLineSearchSteps, iterationsWithNontrivialLineSearch);
         }
 
-        private bool ExitCriteriaSatisfied(Vector<double> candidatePoint, Vector<double> gradient)
+        private bool ExitCriteriaSatisfied(Vector<double> gradient)
         {
             return gradient.Norm(2.0) < GradientTolerance;
         }