Very preliminary Brent and Powell minimizations.

src/Numerics/Optimization/BrentMinimizer.cs

@@ -5,7 +5,301 @@ using System.Text;
 
 namespace MathNet.Numerics.Optimization
 {
+    /// <summary>
+    /// Options for Brent Minimization.
+    /// </summary>
+    public class BrentOptions
+    {
+        public int MaximumIterations = 1000;
+
+        public double FunctionTolerance = 1e-4;
+    }
+    
+    /// <summary>
+    /// Minimizes f(p) where p is a model parameter scalar, i.e. a line-search.
+    /// </summary>
     public class BrentMinimizer
     {
+        const double verySmallNumber = 1e-21, goldenRatio = 1.618034, minimumTolerance = 1.0e-11;
+        const double growLimit = 110.0, conjugateGradient = 0.3819660;
+        const int maxIterations = 500;
+
+        int bracketFunctionCalls;
+        public int FunctionCalls, Iterations;
+        double minPoint, minFunction;
+        Func<double, double> function;
+        double pointA, pointB, pointC;
+        double functionA, functionB, functionC;
+
+        public double Tolerance { get; set; }
+
+        public BrentMinimizer(Func<double, double> function)
+        {
+            this.function = function;
+            Tolerance = 1e-4;
+        }
+
+        public int Search(out double minPoint, out double minFunction)
+        {
+            bracketFunctionCalls = 0;
+            UpdateBracketInterval(0, 1);
+            FunctionCalls = 0;
+            Iterations = 0;
+            int result = BrentMinimize();
+            FunctionCalls += bracketFunctionCalls;
+            minPoint = this.minPoint;
+            minFunction = this.minFunction;
+            return result;
+        }
+
+        private int UpdateBracketInterval(double pointAStart, double pointBStart)
+        {
+            int iterations = 0;
+            pointA = pointAStart;
+            pointB = pointBStart;
+            int maxIterations;
+            maxIterations = 1000;
+            functionA = function(pointA);
+            functionB = function(pointB);
+            double temp;
+            if (functionA < functionB) // Swap points over 
+            {
+                temp = functionA; functionA = functionB; functionB = temp;
+                temp = pointA; pointA = pointB; pointB = temp;
+            }
+            pointC = pointB + goldenRatio * (pointB - pointA);
+            functionC = function(pointC);
+            bracketFunctionCalls = 3; iterations = 0;
+            double temp1, temp2, value, denom;
+            while (functionC < functionB)
+            {
+                double pointW, functionW, wlim;
+                temp1 = (pointB - pointA) * (functionB - functionC);
+                temp2 = (pointB - pointC) * (functionB - functionA);
+                value = temp2 - temp1;
+                if (Math.Abs(value) < verySmallNumber) denom = 2.0 * verySmallNumber;
+                else denom = 2.0 * value;
+                pointW = pointB - ((pointB - pointC) * temp2 - (pointB - pointA) * temp1) / denom;
+                wlim = pointB + growLimit * (pointC - pointB);
+                if (iterations > maxIterations) return 1;
+                iterations++;
+                if ((pointW - pointC) * (pointB - pointW) > 0.0)
+                {
+                    functionW = function(pointW);
+                    bracketFunctionCalls++;
+                    if (functionW < functionC)
+                    {
+                        pointA = pointB; pointB = pointW;
+                        functionA = functionB; functionB = functionW;
+                        return 0;
+                    }
+                    else if (functionW > functionB)
+                    {
+                        pointC = pointW; functionC = functionW;
+                        return 0;
+                    }
+                    pointW = pointC + goldenRatio * (pointC - pointB);
+                    functionW = function(pointW);
+                    bracketFunctionCalls++;
+                }
+                else if ((pointW - wlim) * (wlim - pointC) >= 0.0)
+                {
+                    pointW = wlim;
+                    functionW = function(pointW);
+                    bracketFunctionCalls++;
+                }
+                else if ((pointW - wlim) * (pointC - pointW) > 0.0)
+                {
+                    functionW = function(pointW);
+                    bracketFunctionCalls++;
+                    if (functionW < functionC)
+                    {
+                        pointB = pointC; pointC = pointW;
+                        pointW = pointC + goldenRatio * (pointC - pointB);
+                        functionB = functionC; functionC = functionW;
+                        functionW = function(pointW);
+                        bracketFunctionCalls++;
+                    }
+                }
+                else
+                {
+                    pointW = pointC + goldenRatio * (pointC - pointB);
+                    functionW = function(pointW);
+                    bracketFunctionCalls++;
+                }
+                pointA = pointB; pointB = pointC; pointC = pointW;
+                functionA = functionB; functionB = functionC; functionC = functionW;
+            }
+            return 0;
+        }
+
+        // Find the minimum of the function using the Brent method with the current 
+        // bracketing interval. 
+        private int BrentMinimize()
+        {
+            int result = 0;
+            double x, w, v, fx, fw, fv;
+            double a, b, deltax, rat;
+            double cg = conjugateGradient;
+            x = w = v = pointB; // x is the point with lowest function value encountered
+            fw = fv = fx = function(x);
+            if (pointA < pointC)
+            {
+                a = pointA; b = pointC;
+            }
+            else
+            {
+                a = pointC; b = pointA;
+            }
+            deltax = 0.0;
+            FunctionCalls = 1;
+            Iterations = 0;
+            rat = 0;
+            while (Iterations < maxIterations)
+            {
+                double tol1, tol2, xmin, fval, xmid;
+                double temp1, temp2, p;
+                double u, fu, dx_temp;
+                tol1 = Tolerance * Math.Abs(x) + minimumTolerance;
+                tol2 = 2.0 * tol1;
+                xmid = 0.5 * (a + b);
+                if (Math.Abs(x - xmid) < (tol2 - 0.5 * (b - a)))        // check for convergence 
+                {
+                    xmin = x; fval = fx;
+                    result = 1;
+                    break;
+                }
+                if (Math.Abs(deltax) <= tol1)
+                {
+                    if (x >= xmid) deltax = a - x;        // do a golden section step 
+                    else deltax = b - x;
+                    rat = cg * deltax;
+                }
+                else // do a parabolic step 
+                {
+                    temp1 = (x - w) * (fx - fv);
+                    temp2 = (x - v) * (fx - fw);
+                    p = (x - v) * temp2 - (x - w) * temp1;
+                    temp2 = 2.0 * (temp2 - temp1);
+                    if (temp2 > 0.0) p = -p;
+                    temp2 = Math.Abs(temp2);
+                    dx_temp = deltax;
+                    deltax = rat;
+                    // check parabolic fit 
+                    if ((p > temp2 * (a - x)) && (p < temp2 * (b - x))
+                                        && (Math.Abs(p) < Math.Abs(0.5 * temp2 * dx_temp)))
+                    {
+                        rat = p * 1.0 / temp2; // if parabolic step is useful. 
+                        u = x + rat;
+                        if (((u - a) < tol2) || ((b - u) < tol2))
+                        {
+                            if ((xmid - x) >= 0) rat = tol1;
+                            else rat = -tol1;
+                        }
+                    }
+                    else
+                    {
+                        if (x >= xmid) deltax = a - x; // if it's not do a golden section step 
+                        else deltax = b - x;
+                        rat = cg * deltax;
+                    }
+                }
+
+                if (Math.Abs(rat) < tol1) // update by at least tol1 
+                {
+                    if (rat >= 0) u = x + tol1;
+                    else u = x - tol1;
+                }
+                else
+                {
+                    u = x + rat;
+                }
+                fu = function(u);
+                FunctionCalls++;
+                if (fu > fx)        // if it's bigger than current 
+                {
+                    if (u < x) a = u;
+                    else b = u;
+                    if ((fu <= fw) || (w == x))
+                    {
+                        v = w; w = u; fv = fw; fw = fu;
+                    }
+                    else if ((fu <= fv) || (v == x) || (v == w))
+                    {
+                        v = u; fv = fu;
+                    }
+                }
+                else
+                {
+                    if (u >= x) a = x;
+                    else b = x;
+                    v = w; w = x; x = u;
+                    fv = fw; fw = fx; fx = fu;
+                }
+                Iterations++;
+            }
+            this.minPoint = x;
+            this.minFunction = fx;
+            return result;
+        }
+    }
+
+    /// <summary>
+    /// Minimizes f(p u) where p is a model parameter scalar and u is a direction vector.
+    /// </summary>
+    public class MultiDimensionalBrent
+    {
+        private Func<double[], double> function;
+        private int functionCalls;
+        double[] point;
+        BrentMinimizer lineSearch;
+
+        public MultiDimensionalBrent(Func<double[], double> function)
+        {
+            this.function = function;
+            lineSearch = new BrentMinimizer(this.PointAlongLine);
+            functionCalls = 0;
+        }
+
+        public double Tolerance
+        {
+            get { return lineSearch.Tolerance; }
+            set { lineSearch.Tolerance = value; }
+        }
+
+        public int FunctionCalls
+        {
+            get { return lineSearch.FunctionCalls; }
+        }
+
+        public double[] StartingPoint { get; set; }
+        public double[] Direction { get; set; }
+
+
+        public void SetDimension(int N)
+        {
+            point = new double[N];
+        }
+
+        public int Search(out double[] minPoint, out double minFunction)
+        {
+            double path;
+            int result = lineSearch.Search(out path, out minFunction);
+            for (int i = 0; i < StartingPoint.Length; ++i)
+            {
+                point[i] = StartingPoint[i] + path * Direction[i];
+            }
+            minPoint = point;
+            return result;
+        }
+
+        public double PointAlongLine(double path)
+        {
+            for (int i = 0; i < StartingPoint.Length; ++i)
+            {
+                point[i] = StartingPoint[i] + path * Direction[i];
+            }
+            return function(point);
+        }
     }
 }

src/Numerics/Optimization/PowellMinimizer.cs

@@ -5,8 +5,171 @@ using System.Text;
 
 namespace MathNet.Numerics.Optimization
 {
-    public class PowellSolver
+    /// <summary>
+    /// Minimizes f(p) where p is a vector of model parameters using the Powell method.
+    /// </summary>
+    public class PowellMinimizer
     {
+        public double PointTolerance { get; set; }
+        public double FunctionTolerance { get; set; }
+        int? maxIterations, maxFunctionCalls;
+        double[] minimumPoint;
+        double functionAtMinimum;
+        public int FunctionCalls { get; set; }
+        public int Iterations { get; set; }
+        public int? MaxIterations { get; set; }
+        public int? MaxFunctionCalls { get; set; }
+        public double[] MinimumPoint { get { return minimumPoint; } }
+        public double FunctionValueAtMinimum { get { return functionAtMinimum; } }
 
+        Func<double[], double> function;
+        MultiDimensionalBrent powellLineSearch;
+
+        public PowellMinimizer()
+        {
+            //this.function = function;
+            //powellLineSearch = new MultiDimensionalBrent(function);
+            PointTolerance = 1e-4;
+            FunctionTolerance = 1e-4;
+            MaxIterations = null;
+            MaxFunctionCalls = null;
+            FunctionCalls = 0;
+            Iterations = 0;
+            minimumPoint = null;
+            functionAtMinimum = 0;
+        }
+
+        public double[] CurveFit(double[] x, double[] y, Func<double, double[], double> f,
+            double[] pStart)
+        {
+            // Need to minimize sum of squares of residuals; create this function:
+            Func<double[], double> function = (p) =>
+            {
+                double sum = 0;
+                for (int i = 0; i < x.Length; ++i)
+                {
+                    double temp = y[i] - f(x[i], p);
+                    sum += temp * temp;
+                }
+                return sum;
+            };
+            this.function = function;
+            powellLineSearch = new MultiDimensionalBrent(function);
+            this.Minimize(pStart);
+            return minimumPoint;
+        }
+
+        public int Minimize(double[] p)
+        {
+            // Set line search valuer to use main valuer: 
+            // (this valuer takes starting point, direciton and length and 
+            // returns scalar): 
+            // 
+            int N = p.Length; // number of dimensions 
+            powellLineSearch.SetDimension(N);
+            double fval;
+            FunctionCalls = 0;
+            Iterations = 0;
+            if (maxIterations == null) maxIterations = N * 1000;
+            if (maxFunctionCalls == null) maxFunctionCalls = N * 1000;
+
+            // An array of N directions:
+            double[][] direc = new double[N][];
+            for (int i = 0; i < N; ++i)
+            {
+                direc[i] = new double[N];
+                direc[i][i] = 1.0;
+            }
+            double[] x = p;
+            double[] x1 = (double[])x.Clone(); 
+            powellLineSearch.Tolerance = PointTolerance * 100; // Set tolerance 
+
+            fval = function(x);
+            FunctionCalls++;
+            double[] x2 = new double[N];
+            double fx;
+            double[] direc1 = new double[N];
+            double[] xnew;
+            while (true)
+            {; 
+                fx = fval;
+                int bigind = 0;
+                double delta = 0.0;
+                double fx2;
+                for (int i = 0; i < N; ++i)
+                {
+                    direc1 = direc[i];
+                    fx2 = fval;
+                    powellLineSearch.StartingPoint = x;
+                    powellLineSearch.Direction = direc1;
+
+                    powellLineSearch.Search(out xnew, out fval);
+                    // Do a linesearch with specified starting point and direction. 
+                    FunctionCalls += powellLineSearch.FunctionCalls;
+
+                    for (int j = 0; j < N; ++j) x[j] = xnew[j];
+
+                    if ((fx2 - fval) > delta)
+                    {
+                        delta = fx2 - fval;
+                        bigind = i;
+                    }
+                }
+                Iterations++;
+                if (2.0 * (fx - fval) <= FunctionTolerance * ((Math.Abs(fx) + Math.Abs(fval)) + 1e-20)) break;
+                if (FunctionCalls >= maxFunctionCalls) break;
+                if (Iterations >= maxIterations) break;
+
+                // Construct the extrapolated point  
+                direc1 = new double[N];
+                for (int i = 0; i < N; ++i)
+                {
+                    direc1[i] = x[i] - x1[i];
+                    x2[i] = 2.0 * x[i] - x1[i];
+                    x1[i] = x[i];
+                }
+
+                fx2 = function(x2);
+                FunctionCalls++;
+                if (fx > fx2)
+                {
+                    double t = 2.0 * (fx + fx2 - 2.0 * fval);
+                    double temp = (fx - fval - delta);
+                    t *= temp * temp;
+                    temp = fx - fx2;
+                    t -= delta * temp * temp;
+                    if (t < 0.0)
+                    {
+                        powellLineSearch.StartingPoint = x;
+                        powellLineSearch.Direction = direc1;
+                        powellLineSearch.Search(out xnew, out fval);
+
+                        FunctionCalls += powellLineSearch.FunctionCalls;
+                        direc1 = new double[N];
+                        for (int i = 0; i < N; ++i)
+                        {
+                            direc1[i] = xnew[i] - x[i];
+                            x[i] = xnew[i];
+                        }
+
+                        direc[bigind] = direc[N - 1];
+                        direc[N - 1] = direc1;
+                    }
+                }
+            }
+            minimumPoint = (double[])x.Clone();
+            functionAtMinimum = fx;
+            
+            // Find out what happened: 
+            if (FunctionCalls >= maxFunctionCalls)
+            {
+                return 1; // Max function calls exceeded 
+            }
+            if (Iterations >= maxIterations)
+            {
+                return 2; // Max iterations exceeded 
+            }
+            return 0; // all good 
+        }
     }
 }

src/UnitTests/OptimizationTests/FunctionMinimizationTests.cs

@@ -0,0 +1,61 @@
+// <copyright file="FunctionMinimizationTests.cs" company="Math.NET">
+// 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-2013 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 FunctionMinimizationTests
+    {
+        [Test]
+        public void PowellCurveFit()
+        {
+            // y = b1*(1-exp[-b2*x])  +  e
+            var xin = new double[] { 1, 2, 3, 5, 7, 10 };
+            var yin = new double[] { 109, 149, 149, 191, 213, 224 };
+
+            Func<double, double[], double> function = (x, p) => p[0] * (1 - Math.Exp(-p[1] * x));
+
+            var minimizer = new PowellMinimizer();
+
+            var popt = minimizer.CurveFit(xin, yin, function, new double[] { 1, 1 }); // 100, 0.75
+
+            double[] expected = new double[] { 2.1380940889E+02, 5.4723748542E-01 };
+
+            double residual = 0;
+            for (int i = 0; i < yin.Length; ++i) residual += (yin[i] - function(xin[i], popt)) * (yin[i] - function(xin[i], popt));
+            //Assert.AreEqual(3, Brent.FindRoot(f2, 2.1, 3.4, 0.001, 50), 0.001);
+        }
+
+    }
+}

src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs

@@ -61,6 +61,5 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
             for (int i = 0; i < yin.Length; ++i) residual += (yin[i] - function(xin[i], popt)) * (yin[i] - function(xin[i], popt)); 
             //Assert.AreEqual(3, Brent.FindRoot(f2, 2.1, 3.4, 0.001, 50), 0.001);
         }
-
     }
 }

src/UnitTests/UnitTests.csproj

@@ -353,6 +353,7 @@
     <Compile Include="NumberTheoryTests\GcdRelatedTest.cs" />
     <Compile Include="NumberTheoryTests\GcdRelatedTestBigInteger.cs" />
     <Compile Include="NumberTheoryTests\IntegerTheoryTest.cs" />
+    <Compile Include="OptimizationTests\FunctionMinimizationTests.cs" />
     <Compile Include="OptimizationTests\NonLinearLeastSquaresTest.cs" />
     <Compile Include="RootFindingTests\BisectionTest.cs" />
     <Compile Include="PermutationTest.cs" />