From 9610bf032f37a129504b9fc468129af0ef3e013d Mon Sep 17 00:00:00 2001 From: joemoorhouse Date: Sun, 10 Nov 2013 21:36:24 +0000 Subject: [PATCH] Very preliminary Brent and Powell minimizations. --- src/Numerics/Optimization/BrentMinimizer.cs | 294 ++++++++++++++++++ src/Numerics/Optimization/PowellMinimizer.cs | 165 +++++++++- .../FunctionMinimizationTests.cs | 61 ++++ .../NonLinearLeastSquaresTest.cs | 1 - src/UnitTests/UnitTests.csproj | 1 + 5 files changed, 520 insertions(+), 2 deletions(-) create mode 100644 src/UnitTests/OptimizationTests/FunctionMinimizationTests.cs diff --git a/src/Numerics/Optimization/BrentMinimizer.cs b/src/Numerics/Optimization/BrentMinimizer.cs index 6b6de901..3d029903 100644 --- a/src/Numerics/Optimization/BrentMinimizer.cs +++ b/src/Numerics/Optimization/BrentMinimizer.cs @@ -5,7 +5,301 @@ using System.Text; namespace MathNet.Numerics.Optimization { + /// + /// Options for Brent Minimization. + /// + public class BrentOptions + { + public int MaximumIterations = 1000; + + public double FunctionTolerance = 1e-4; + } + + /// + /// Minimizes f(p) where p is a model parameter scalar, i.e. a line-search. + /// 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 function; + double pointA, pointB, pointC; + double functionA, functionB, functionC; + + public double Tolerance { get; set; } + + public BrentMinimizer(Func 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; + } + } + + /// + /// Minimizes f(p u) where p is a model parameter scalar and u is a direction vector. + /// + public class MultiDimensionalBrent + { + private Func function; + private int functionCalls; + double[] point; + BrentMinimizer lineSearch; + + public MultiDimensionalBrent(Func 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); + } } } diff --git a/src/Numerics/Optimization/PowellMinimizer.cs b/src/Numerics/Optimization/PowellMinimizer.cs index a1377516..1927e44a 100644 --- a/src/Numerics/Optimization/PowellMinimizer.cs +++ b/src/Numerics/Optimization/PowellMinimizer.cs @@ -5,8 +5,171 @@ using System.Text; namespace MathNet.Numerics.Optimization { - public class PowellSolver + /// + /// Minimizes f(p) where p is a vector of model parameters using the Powell method. + /// + 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 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 f, + double[] pStart) + { + // Need to minimize sum of squares of residuals; create this function: + Func 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 + } } } diff --git a/src/UnitTests/OptimizationTests/FunctionMinimizationTests.cs b/src/UnitTests/OptimizationTests/FunctionMinimizationTests.cs new file mode 100644 index 00000000..52975f64 --- /dev/null +++ b/src/UnitTests/OptimizationTests/FunctionMinimizationTests.cs @@ -0,0 +1,61 @@ +// +// 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. +// + +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 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); + } + + } +} diff --git a/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs b/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs index 325275ef..98bc6276 100644 --- a/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs +++ b/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); } - } } diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 0776c7da..8c2292b6 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -353,6 +353,7 @@ +