bdodson
11 years ago
committed by
Erik Ovegard
Erik Ovegard
5 changed files with 301 additions and 0 deletions
@ -0,0 +1,108 @@ |
|||
// <copyright file="BfgsSolver.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 System.Collections.Generic; |
|||
using MathNet.Numerics.LinearAlgebra; |
|||
using MathNet.Numerics.LinearAlgebra.Double; |
|||
using System.Linq; |
|||
using System.Text; |
|||
|
|||
namespace MathNet.Numerics.RootFinding |
|||
{ |
|||
/// <summary>
|
|||
/// Broyden-Fletcher-Goldfarb-Shanno solver for finding function minima
|
|||
/// See http://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm
|
|||
/// Inspired by implementation: https://github.com/PatWie/CppNumericalSolvers/blob/master/src/BfgsSolver.cpp
|
|||
/// </summary>
|
|||
public static class BfgsSolver |
|||
{ |
|||
private const double gradientTolerance = 1e-5; |
|||
private const int maxIterations = 100000; |
|||
|
|||
/// <summary>
|
|||
/// Finds a minimum of a function by the BFGS quasi-Newton method
|
|||
/// This uses the function and it's gradient (partial derivatives in each direction) and approximates the Hessian
|
|||
/// </summary>
|
|||
/// <param name="initialGuess">An initial guess</param>
|
|||
/// <param name="functionValue">Evaluates the function at a point</param>
|
|||
/// <param name="functionGradient">Evaluates the gradient of the function at a point</param>
|
|||
/// <returns>The minimum found</returns>
|
|||
public static Vector<double> Solve(Vector initialGuess, Func<Vector<double>, double> functionValue, Func<Vector<double>, Vector<double>> functionGradient) |
|||
{ |
|||
int dim = initialGuess.Count; |
|||
int iter = 0; |
|||
// H represents the approximation of the inverse hessian matrix
|
|||
// it is updated via the Sherman–Morrison formula (http://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula)
|
|||
Matrix<double> H = DenseMatrix.CreateIdentity(dim); |
|||
|
|||
Vector<double> x = initialGuess; |
|||
Vector<double> x_old = x; |
|||
Vector<double> grad; |
|||
|
|||
do |
|||
{ |
|||
// search along the direction of the gradient
|
|||
grad = functionGradient(x); |
|||
Vector<double> p = -1 * H * grad; |
|||
double rate = WolfeRule.LineSearch(x, p, functionValue, functionGradient); |
|||
x = x + rate * p; |
|||
Vector<double> grad_old = grad; |
|||
|
|||
// update the gradient
|
|||
grad = functionGradient(x); |
|||
|
|||
Vector<double> s = x - x_old; |
|||
Vector<double> y = grad - grad_old; |
|||
|
|||
double rho = 1.0 / (y * s); |
|||
if (iter == 0) |
|||
{ |
|||
// set up an initial hessian
|
|||
H = (y * s) / (y * y) * DenseMatrix.CreateIdentity(dim); |
|||
} |
|||
|
|||
var sM = s.ToColumnMatrix(); |
|||
var yM = y.ToColumnMatrix(); |
|||
|
|||
// Update the estimate of the hessian
|
|||
H = H |
|||
- rho * (sM * (yM.TransposeThisAndMultiply(H)) + (H * yM).TransposeAndMultiply(sM)) |
|||
+ rho * rho * (y.DotProduct(H * y) + 1.0 / rho) * (sM.TransposeAndMultiply(sM)); |
|||
x_old = x; |
|||
iter++; |
|||
} |
|||
while ((grad.InfinityNorm() > gradientTolerance) && (iter < maxIterations)); |
|||
|
|||
return x; |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,113 @@ |
|||
// <copyright file="WolfeRule.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 System.Collections.Generic; |
|||
using System.Linq; |
|||
using MathNet.Numerics.LinearAlgebra; |
|||
|
|||
namespace MathNet.Numerics.RootFinding |
|||
{ |
|||
/// <summary>
|
|||
/// Performs an inexact line search. This is used as a part of quasi-Newton optimization methods to figure
|
|||
/// out how far to move along a certain gradient.
|
|||
/// See http://en.wikipedia.org/wiki/Wolfe_conditions
|
|||
/// Inspired by implementation: https://github.com/PatWie/CppNumericalSolvers/blob/master/src/linesearch/WolfeRule.h
|
|||
/// </summary>
|
|||
internal static class WolfeRule |
|||
{ |
|||
/// <summary>
|
|||
/// Searches along a line to satisfy the Wolfe conditions (inexact search for minimum)
|
|||
/// </summary>
|
|||
/// <param name="x0">Starting point of search</param>
|
|||
/// <param name="z">Search direction</param>
|
|||
/// <param name="functionValue">Evaluates the function being minimized</param>
|
|||
/// <param name="functionGradient">Evaluates the gradient of the function</param>
|
|||
/// <param name="alphaInit">Initial value for the coefficient of z (distance to travel in z direction)</param>
|
|||
/// <returns></returns>
|
|||
public static double LineSearch( |
|||
Vector<double> x0, |
|||
Vector<double> z, |
|||
Func<Vector<double>, double> functionValue, |
|||
Func<Vector<double>, Vector<double>> functionGradient, |
|||
float alphaInit = 1) |
|||
{ |
|||
Vector<double> x = x0; |
|||
|
|||
// evaluate phi(0)
|
|||
double phi0 = functionValue(x0); |
|||
|
|||
// evaluate phi'(0)
|
|||
Vector<double> grad = functionGradient(x); |
|||
double phi0_dash = z * grad; |
|||
|
|||
double alpha = alphaInit; |
|||
|
|||
bool decrease_direction = true; |
|||
|
|||
// 200 guesses max
|
|||
for (int iter = 0; iter < 200; ++iter) { |
|||
|
|||
// new guess for phi(alpha)
|
|||
Vector<double> x_candidate = x + alpha * z; |
|||
double phi = functionValue(x_candidate); |
|||
|
|||
// decrease condition invalid --> shrink interval
|
|||
if (phi > phi0 + 0.0001 * alpha * phi0_dash) |
|||
{ |
|||
alpha *= 0.5; |
|||
decrease_direction = false; |
|||
} |
|||
else |
|||
{ |
|||
// valid decrease --> test strong wolfe condition
|
|||
Vector<double> grad2 = functionGradient(x_candidate); |
|||
double phi_dash = z * grad2; |
|||
|
|||
// curvature condition invalid ?
|
|||
if ((phi_dash < 0.9 * phi0_dash) || !decrease_direction) { |
|||
// increase interval
|
|||
alpha *= 4.0; |
|||
} |
|||
else { |
|||
// both condition are valid --> we are happy
|
|||
x = x_candidate; |
|||
grad = grad2; |
|||
phi0 = phi; |
|||
return alpha; |
|||
} |
|||
} |
|||
} |
|||
|
|||
|
|||
return alpha; |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,77 @@ |
|||
// <copyright file="BfgsTest.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 System.Collections.Generic; |
|||
using System.Linq; |
|||
using System.Text; |
|||
using System.Threading.Tasks; |
|||
using MathNet.Numerics.LinearAlgebra; |
|||
using NUnit.Framework; |
|||
using MathNet.Numerics.LinearAlgebra.Double; |
|||
using MathNet.Numerics.RootFinding; |
|||
|
|||
namespace MathNet.Numerics.UnitTests.RootFindingTests |
|||
{ |
|||
[TestFixture, Category("RootFinding")] |
|||
internal class BfgsTest |
|||
{ |
|||
private const double precision = 1e-4; |
|||
|
|||
[Test] |
|||
public void MinimizeRosenbrock() |
|||
{ |
|||
CheckRosenbrock(15.0, 8.0, expectedMin: 0.0); |
|||
CheckRosenbrock(-1.2, 1.0, expectedMin: 0.0); |
|||
CheckRosenbrock(-1.2, 100.0, expectedMin: 0.0); |
|||
} |
|||
|
|||
private static void CheckRosenbrock(double a, double b, double expectedMin) |
|||
{ |
|||
var x = BfgsSolver.Solve(new DenseVector(new[] { a, b }), Rosenbrock, RosenbrockGradient); |
|||
Precision.AlmostEqual(expectedMin, Rosenbrock(x), precision); |
|||
} |
|||
|
|||
private static double Rosenbrock(Vector<double> x) |
|||
{ |
|||
double t1 = (1 - x[0]); |
|||
double t2 = (x[1] - x[0] * x[0]); |
|||
return t1 * t1 + 100 * t2 * t2; |
|||
} |
|||
|
|||
private static Vector<double> RosenbrockGradient(Vector<double> x) |
|||
{ |
|||
var grad = new DenseVector(2); |
|||
grad[0] = -2 * (1 - x[0]) + 200 * (x[1] - x[0] * x[0]) * (-2 * x[0]); |
|||
grad[1] = 200 * (x[1] - x[0] * x[0]); |
|||
return grad; |
|||
} |
|||
} |
|||
} |
|||
Loading…
Reference in new issue