@ -1,152 +1,173 @@
using System ; // <copyright file="BfgsMinimizer.cs" company="Math.NET">
using System.Collections.Generic ; // Math.NET Numerics, part of the Math.NET Project
using System.Linq ; // http://numerics.mathdotnet.com
using System.Text ; // 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.LinearAlgebra ; using MathNet.Numerics.LinearAlgebra ;
namespace MathNet.Numerics.Optimization namespace MathNet.Numerics.Optimization
{ {
public class BfgsMinimizer public class BfgsMinimizer
{ {
public double GradientTolerance { get ; set ; } public double GradientTolerance { get ; set ; }
public double ParameterTolerance { get ; set ; } public double ParameterTolerance { get ; set ; }
public int MaximumIterations { get ; set ; } public int MaximumIterations { get ; set ; }
public BfgsMinimizer ( double gradient_tolerance , double parameter_tolerance , int maximum_iterations = 1 0 0 0 ) public BfgsMinimizer ( double gradientTolerance , double parameterTolerance , int maximumIterations = 1 0 0 0 )
{ {
this . GradientTolerance = gradient_tolerance ; GradientTolerance = gradientT olerance ;
this . ParameterTolerance = parameter_tolerance ; ParameterTolerance = parameterT olerance ;
this . MaximumIterations = maximum_iterations ; MaximumIterations = maximumI terations ;
} }
public MinimizationOutput FindMinimum ( IObjectiveFunction objective , Vector < double > initial_guess ) public MinimizationOutput FindMinimum ( IObjectiveFunction objective , Vector < double > initialG uess )
{ {
if ( ! objective . GradientSupported ) if ( ! objective . GradientSupported )
throw new IncompatibleObjectiveException ( "Gradient not supported in objective function, but required for BFGS minimization." ) ; throw new IncompatibleObjectiveException ( "Gradient not supported in objective function, but required for BFGS minimization." ) ;
if ( ! ( objective is ObjectiveChecker ) ) if ( ! ( objective is ObjectiveChecker ) )
objective = new ObjectiveChecker ( objective , this . ValidateObjective , this . ValidateGradient , null ) ; objective = new ObjectiveChecker ( objective , ValidateObjective , ValidateGradient , null ) ;
IEvaluation initialEval = objective . Evaluate ( initialGuess ) ;
IEvaluation initial_eval = objective . Evaluate ( initial_guess ) ;
// Check that we're not already done
// Check that we're not already done
ExitCondition current_exit_c ondition = this . ExitCriteriaSatisfied ( initial_e val , null ) ; ExitCondition currentExitCondition = ExitCriteriaSatisfied ( initialEval , null ) ;
if ( current_exit_c ondition ! = ExitCondition . None ) if ( currentExitC ondition ! = ExitCondition . None )
return new MinimizationOutput ( initial_e val , 0 , current_exit_c ondition ) ; return new MinimizationOutput ( initialE val , 0 , currentExitC ondition ) ;
// Set up line search algorithm
// Set up line search algorithm
var line_s earcher = new WeakWolfeLineSearch ( 1e-4 , 0.9 , this . ParameterTolerance , max_iterations : 1 0 0 0 ) ; var lineS earcher = new WeakWolfeLineSearch ( 1e-4 , 0.9 , ParameterTolerance , maxIterations : 1 0 0 0 ) ;
// Declare state variables
// Declare state variables
IEvaluation candidate_point , previous_point ; Vector < double > gradient ;
double step_size ;
Vector < double > gradient , step , search_direction ;
Matrix < double > inverse_pseudo_hessian ;
// First step
// First step
inverse_pseudo_h essian= Matrix < double > . Build . DiagonalIdentity ( initial_g uess . Count ) ; Matrix < double > inversePseudoHessian = Matrix < double > . Build . DiagonalIdentity ( initialGuess . Count ) ;
search_direction = - initial_e val. Gradient ; Vector < double > searchDirection = - initialE val. Gradient ;
step_s ize= 1 0 0 * this . GradientTolerance / ( search_direction * search_d irection ) ; double stepS ize= 1 0 0 * GradientTolerance / ( searchDirection * searchD irection) ;
LineSearchOutput result ; LineSearchOutput result ;
try try
{ {
result = line_s earcher . FindConformingStep ( objective , initial_eval , search_direction , step_s ize ) ; result = lineS earcher . FindConformingStep ( objective , initialEval , searchDirection , stepS ize ) ;
} }
catch ( Exception e ) catch ( Exception e )
{ {
throw new InnerOptimizationException ( "Line search failed." , e ) ; throw new InnerOptimizationException ( "Line search failed." , e ) ;
} }
previous_point = initial_e val ; IEvaluation previousPoint = initialE val ;
candidate_p oint= result . FunctionInfoAtMinimum ; IEvaluation candidateP oint= result . FunctionInfoAtMinimum ;
gradient = candidate_p oint . Gradient ; gradient = candidateP oint . Gradient ;
step = candidate_p oint . Point - initial_g uess ; Vector < double > step = candidateP oint . Point - initialG uess ;
step_s ize = result . FinalStep ; stepS ize = result . FinalStep ;
// Subsequent steps
// Subsequent steps
int iterations ; int iterations ;
int total_line_search_steps = result . Iterations ; int totalLineSearchSteps = result . Iterations ;
int iterations_with_nontrivial_line_search = result . Iterations > 0 ? 0 : 1 ; int iterationsWithNontrivialLineSearch = result . Iterations > 0 ? 0 : 1 ;
int steepest_descent_resets = 0 ; for ( iterations = 1 ; iterations < MaximumIterations ; + + iterations )
for ( iterations = 1 ; iterations < this . MaximumIterations ; + + iterations )
{ {
var y = candidate_p oint . Gradient - previous_p oint . Gradient ; var y = candidateP oint . Gradient - previousP oint . Gradient ;
double sy = step * y ; double sy = step * y ;
inverse_pseudo_hessian = inverse_pseudo_h essian + ( ( sy + y * inverse_pseudo_hessian * y ) / Math . Pow ( sy , 2.0 ) ) * step . OuterProduct ( step ) - ( ( inverse_pseudo_hessian * y . ToColumnMatrix ( ) ) * step . ToRowMatrix ( ) + step . ToColumnMatrix ( ) * ( y . ToRowMatrix ( ) * inverse_pseudo_h essian ) ) * ( 1.0 / sy ) ; inversePseudoHessian = inversePseudoH essian + ( ( sy + y * inversePseudoHessian * y ) / Math . Pow ( sy , 2.0 ) ) * step . OuterProduct ( step ) - ( ( inversePseudoHessian * y . ToColumnMatrix ( ) ) * step . ToRowMatrix ( ) + step . ToColumnMatrix ( ) * ( y . ToRowMatrix ( ) * inversePseudoH essian) ) * ( 1.0 / sy ) ;
search_direction = - inverse_pseudo_hessian * candidate_p oint . Gradient ; searchDirection = - inversePseudoHessian * candidateP oint . Gradient ;
if ( search_direction * candidate_p oint . Gradient > = - this . GradientTolerance * this . GradientTolerance ) if ( searchDirection * candidateP oint . Gradient > = - GradientTolerance * GradientTolerance )
{ {
search_direction = - candidate_point . Gradient ; searchDirection = - candidatePoint . Gradient ;
inverse_pseudo_hessian = Matrix < double > . Build . DiagonalIdentity ( initial_guess . Count ) ; inversePseudoHessian = Matrix < double > . Build . DiagonalIdentity ( initialGuess . Count ) ;
steepest_descent_resets + = 1 ;
} }
try try
{ {
result = line_s earcher . FindConformingStep ( objective , candidate_point , search_d irection , 1.0 ) ; result = lineS earcher . FindConformingStep ( objective , candidatePoint , searchD irection , 1.0 ) ;
} }
catch ( Exception e ) catch ( Exception e )
{ {
throw new InnerOptimizationException ( "Line search failed." , e ) ; throw new InnerOptimizationException ( "Line search failed." , e ) ;
} }
iterations_with_nontrivial_line_s earch + = result . Iterations > 0 ? 1 : 0 ; iterationsWithNontrivialLineS earch + = result . Iterations > 0 ? 1 : 0 ;
total_line_search_s teps + = result . Iterations ; totalLineSearchS teps + = result . Iterations ;
step_s ize = result . FinalStep ; stepS ize = result . FinalStep ;
step = result . FunctionInfoAtMinimum . Point - candidate_p oint . Point ; step = result . FunctionInfoAtMinimum . Point - candidateP oint . Point ;
previous_point = candidate_p oint ; previousPoint = candidateP oint ;
candidate_p oint = result . FunctionInfoAtMinimum ; candidateP oint = result . FunctionInfoAtMinimum ;
current_exit_c ondition = this . ExitCriteriaSatisfied ( candidate_point , previous_p oint ) ; currentExitC ondition = ExitCriteriaSatisfied ( candidatePoint , previousP oint ) ;
if ( current_exit_c ondition ! = ExitCondition . None ) if ( currentExitC ondition ! = ExitCondition . None )
break ; break ;
} }
if ( iterations = = this . MaximumIterations ) if ( iterations = = MaximumIterations )
throw new MaximumIterationsException ( String . Format ( "Maximum iterations ({0}) reached." , this . MaximumIterations ) ) ; throw new MaximumIterationsException ( String . Format ( "Maximum iterations ({0}) reached." , MaximumIterations ) ) ;
return new MinimizationWithLineSearchOutput ( candidate_p oint , iterations , current_exit_condition , total_line_search_steps , iterations_with_nontrivial_line_s earch ) ; return new MinimizationWithLineSearchOutput ( candidateP oint , iterations , currentExitCondition , totalLineSearchSteps , iterationsWithNontrivialLineS earch ) ;
} }
private ExitCondition ExitCriteriaSatisfied ( IEvaluation candidate_point , IEvaluation last_p oint ) ExitCondition ExitCriteriaSatisfied ( IEvaluation candidatePoint , IEvaluation lastP oint )
{ {
Vector < double > rel_g rad = new MathNet . Numerics . LinearAlgebra . Double . DenseVector ( candidate_p oint . Point . Count ) ; Vector < double > relG rad = new LinearAlgebra . Double . DenseVector ( candidateP oint . Point . Count ) ;
double relative_g radient = 0.0 ; double relativeG radient = 0.0 ;
double normalizer = Math . Max ( Math . Abs ( candidate_p oint . Value ) , 1.0 ) ; double normalizer = Math . Max ( Math . Abs ( candidateP oint . Value ) , 1.0 ) ;
for ( int ii = 0 ; ii < rel_g rad . Count ; + + ii ) for ( int ii = 0 ; ii < relG rad . Count ; + + ii )
{ {
double tmp = candidate_p oint . Gradient [ ii ] * Math . Max ( Math . Abs ( candidate_p oint . Point [ ii ] ) , 1.0 ) / normalizer ; double tmp = candidateP oint . Gradient [ ii ] * Math . Max ( Math . Abs ( candidateP oint . Point [ ii ] ) , 1.0 ) / normalizer ;
relative_g radient = Math . Max ( relative_g radient , Math . Abs ( tmp ) ) ; relativeG radient = Math . Max ( relativeG radient , Math . Abs ( tmp ) ) ;
} }
if ( relative_gradient < this . GradientTolerance ) if ( relativeGradient < GradientTolerance )
{ {
return ExitCondition . RelativeGradient ; return ExitCondition . RelativeGradient ;
} }
if ( last_p oint ! = null ) if ( lastP oint ! = null )
{ {
double most_p rogress = 0.0 ; double mostP rogress = 0.0 ;
for ( int ii = 0 ; ii < candidate_p oint . Point . Count ; + + ii ) for ( int ii = 0 ; ii < candidateP oint . Point . Count ; + + ii )
{ {
var tmp = Math . Abs ( candidate_p oint . Point [ ii ] - last_p oint . Point [ ii ] ) / Math . Max ( Math . Abs ( last_p oint . Point [ ii ] ) , 1.0 ) ; var tmp = Math . Abs ( candidateP oint . Point [ ii ] - lastP oint . Point [ ii ] ) / Math . Max ( Math . Abs ( lastP oint . Point [ ii ] ) , 1.0 ) ;
most_p rogress = Math . Max ( most_p rogress , tmp ) ; mostP rogress = Math . Max ( mostP rogress , tmp ) ;
} }
if ( most_p rogress < this . ParameterTolerance ) if ( mostP rogress < ParameterTolerance )
{ {
return ExitCondition . LackOfProgress ; return ExitCondition . LackOfProgress ;
} }
} }
return ExitCondition . None ; return ExitCondition . None ;
} }
private void ValidateGradient ( IEvaluation eval ) void ValidateGradient ( IEvaluation eval )
{ {
foreach ( var x in eval . Gradient ) foreach ( var x in eval . Gradient )
{ {
@ -155,7 +176,7 @@ namespace MathNet.Numerics.Optimization
} }
} }
private void ValidateObjective ( IEvaluation eval ) void ValidateObjective ( IEvaluation eval )
{ {
if ( Double . IsNaN ( eval . Value ) | | Double . IsInfinity ( eval . Value ) ) if ( Double . IsNaN ( eval . Value ) | | Double . IsInfinity ( eval . Value ) )
throw new EvaluationException ( "Non-finite objective function returned." , eval ) ; throw new EvaluationException ( "Non-finite objective function returned." , eval ) ;
Christoph Ruegg
13 years ago