Multidimensional Optimization

Finding the minimum or maximum of a function in many variables is one of the most common problems in numerical computing. It is not surprising that a wide range of techniques exist, each with its advantages and disadvantages. Numerics.NET implements a large number of these techniques, to offer solutions for the broadest possible range of applications.

The Downhill Simplex Method of Nelder and Mead

The Nelder-Mead algorithm, often also called the downhill simplex method, is a simple algorithm that produces reasonable results when no derivatives are available. A simplex is a generalization of a triangle (2 dimensional) or tetrahedron (3 dimensional) to arbitrary dimensions. A simplex in n dimensions has n+1 vertices. This is the smallest number of points needed to define a region in n-dimensional space. The algorithm works by letting this simplex 'travel' downhill on the surface defined by the objective function through reflection. If no lower point can be found, the simplex is contracted towards the lowest point.

The algorithm defines its own versions of the ValueTest and SolutionTest properties. The value test checks the objective-value spread across the simplex. The solution test checks coordinate widths across every vertex, using an absolute tolerance by default. The default convergence collection requires both enabled tests to pass, including any required successive passes. Equal values at two vertices do not imply that the simplex is small. The gradient test does not apply.

Convergence of the Nelder-Mead method is linear at best. It can fail in certain pathological situations. If derivatives are available, and the derivatives are continuous, then one of the other methods should be used.

Conjugate Gradient Optimizers

Two algorithms are available that search successively in a direction from a well-defined set of directions. For each direction, a one-dimensional search for an extremum is performed along the direction.

The conjugate gradient method, implemented by the ConjugateGradientOptimizer class, uses a series of search directions that are optimal for a quadratic objective function. The method requires that the gradient of the objective function be available. There are four variants of the formula to generate a new search direction, giving rise to the Fletcher-Reeves and Polak-Ribière methods. The variant is selected by passing the appropriate ConjugateGradientMethod to the constructor, or setting the Method property. It is found that, in general, the Polak-Ribière method tends to converge somewhat faster. The positive variant clamps a negative direction coefficient to zero and is the default. The Hager-Zhang variant is also available.

Powell's method, implemented by the PowellOptimizer class, is a variation of the conjugate gradient method that does not require knowledge of the gradient of the objective function. Instead, the directions are generated from the geometrical properties of the successive approximations. A drawback of Powell's method is that it usually requires more function evaluations, as the line searches need to be exact for the method to converge.

Newton Steps Without a Dense Hessian

Use NewtonConjugateGradientOptimizer when first derivatives and Hessian-vector products are available. A Hessian-vector product tells the method how the gradient changes in one direction. It can be much cheaper to compute than all second derivatives. The optimizer uses these products to find Newton steps with conjugate gradients. This is a different algorithm from the first-derivative conjugate-gradient optimizer described above.

Supply the objective and its gradient through the existing optimizer callbacks, plus Hessian. The function receives the current point and reusable linear operator, initially null. Return an explicit matrix or an operator supplying Hessian-vector multiplication. Keep it valid until the next callback; overwrite changed matrix entries, including entries becoming zero, when reusing a matrix. Gradients, products and Hessian entries use the original objective sign, including for maximization.

The following example solves a 1,000-variable quadratic problem with bounds. Its Hessian is the identity, so a product only needs to copy the direction. The solution has every component equal to one.

const int n = 1000;
NewtonConjugateGradientOptimizer optimizer = new NewtonConjugateGradientOptimizer
{
    InitialGuess = Vector.Zeros<double>(n),
    LowerBounds = Vector.FromFunction(n, i => 0.0),
    UpperBounds = Vector.FromFunction(n, i => 2.0),
    ObjectiveFunction = x =>
    {
        double value = 0;
        for (int i = 0; i < n; i++) value += .5 * (x[i] - 1) * (x[i] - 1);
        return value;
    },
    FastGradientFunction = (x, g) =>
    {
        for (int i = 0; i < n; i++) g[i] = x[i] - 1;
        return g;
    },
    // The Hessian is the identity: its product copies the direction.
    Hessian = (x, previous) => previous ?? LinearOperator<double>.Create(n, n,
        (operation, direction, result) =>
        {
            result = result ?? Vector.Zeros<double>(n);
            direction.CopyTo(result); return result;
        }),
    HessianDiagonal = (x, result) => result.SetValue(1)
};
Vector<double> solution = optimizer.FindMinimum();
if (optimizer.Status != AlgorithmStatus.Converged)
    throw new InvalidOperationException("Newton-CG did not converge.");
// Every solution component is 1.

Const n As Integer = 1000
Dim optimizer As New NewtonConjugateGradientOptimizer()
optimizer.InitialGuess = Vector.Zeros(Of Double)(n)
optimizer.LowerBounds = Vector.FromFunction(n, Function(i) 0.0)
optimizer.UpperBounds = Vector.FromFunction(n, Function(i) 2.0)
optimizer.ObjectiveFunction = Function(x)
                                  Dim value As Double = 0
                                  For i As Integer = 0 To n - 1
                                      value += 0.5 * (x(i) - 1) * (x(i) - 1)
                                  Next
                                  Return value
                              End Function
optimizer.FastGradientFunction = Function(x, g)
                                      For i As Integer = 0 To n - 1
                                          g(i) = x(i) - 1
                                      Next
                                      Return g
                                  End Function
' The Hessian is the identity: its product copies the direction.
optimizer.Hessian = Function(x, previous)
                        Return If(previous, LinearOperator(Of Double).Create(n, n,
                            Function(operation, direction, result)
                                If result Is Nothing Then result = Vector.Zeros(Of Double)(n)
                                direction.CopyTo(result)
                                Return result
                            End Function))
                    End Function
optimizer.HessianDiagonal = Sub(x, result) result.SetValue(1)
Dim solution As Vector(Of Double) = optimizer.FindMinimum()
If optimizer.Status <> AlgorithmStatus.Converged Then Throw New InvalidOperationException("Newton-CG did not converge.")
' Every solution component is 1.

No code example is currently available or this language may not be supported.

let n = 1000
let optimizer = new NewtonConjugateGradientOptimizer()
optimizer.InitialGuess <- Vector.Zeros<float>(n)
optimizer.LowerBounds <- Vector.FromFunction(n, Func<int,float>(fun i -> 0.0))
optimizer.UpperBounds <- Vector.FromFunction(n, Func<int,float>(fun i -> 2.0))
optimizer.ObjectiveFunction <- Func<_,_>(fun x ->
    let mutable value = 0.0
    for i in 0..n-1 do value <- value + 0.5*(x[i]-1.0)*(x[i]-1.0)
    value)
optimizer.FastGradientFunction <- Func<_,_,_>(fun x g ->
    for i in 0..n-1 do g[i] <- x[i]-1.0
    g)
// The Hessian is the identity: its product copies the direction.
optimizer.Hessian <- Func<_,_,_>(fun x previous ->
    if not (isNull previous) then previous
    else LinearOperator<float>.Create(n, n, Func<_,_,_,_>(fun operation direction result ->
        let output = if isNull result then Vector.Zeros<float>(n) :> Vector<float> else result
        direction.CopyTo(output) |> ignore
        output)))
optimizer.HessianDiagonal <- Action<_,_>(fun x result -> result.SetValue(1.0) |> ignore)
let solution = optimizer.FindMinimum()
if optimizer.Status <> AlgorithmStatus.Converged then invalidOp "Newton-CG did not converge."
// Every solution component is 1.

Optional LowerBounds and UpperBounds are copied at the start of each solve. The initial point is moved inside these bounds before evaluating the objective. For general constraints, use NonlinearProgram instead.

A preconditioner helps conjugate gradients find a Newton step with fewer products. Automatic selection uses incomplete Cholesky for a sparse Hessian matrix, an available diagonal otherwise, and limited secant history when neither is available. Use PreconditioningMethod to choose a method explicitly. Sparse IC/MIC need an explicit sparse Hessian matrix. Sparse matrices supply both entries and products. Keep the stored pattern fixed within a solve, including entries that are currently zero; numerical values can change at every point.

Inspect the returned status. Successful termination requires bound feasibility, projected stationarity and complementarity within OptimalityTolerance, whose default is 1e-6. The ordinary value, step and gradient tests are disabled by default. If a custom test stops at a point that fails the mandatory checks, the status is ConvergedToFalseSolution. These checks do not guarantee parameter accuracy or global optimality for nonconvex functions.

Iteration and objective-evaluation limits, cancellation and repeated solves use the optimizer's usual controls. History and workspaces reset per solve. Callback exceptions propagate; incomplete or nonfinite derivative outputs fail the solve. Optional UseLineSearch enables backtracking when there are no finite bounds. Finite bounds retain the trust-region method.

Quasi-Newton Methods

The full Newton method requires the calculation of second derivatives and the solution of a system of equations in every iteration. This can be very expensive and may not lead to great results, especially when the current point is far from the actual extremum.

Quasi-Newton methods use a numerical approximation to the inverse of the Hessian matrix that is maintained through each iteration. Different variations of the method use different ways to keep the approximation up to date. There are two main variations: the Davidon-Fletcher-Powell method (commonly abbreviated to DFP) and the Broyden-Fletcher-Goldfard-Shanno method (BFGS). The method is selected by passing the appropriate QuasiNewtonMethod to the constructor, or setting the Method property. The BFGS method is usually somewhat superior to the DFP method. It is the default.

The BFGS and DFP methods use a dense approximation to the inverse Hessian. As the number of variables grows larger, both the memory requirements and the computation time per iteration grow considerably. To address this, a limited memory variant of the BFGS method was developed. This is implemented by the LimitedMemoryBfgsOptimizer class. It uses only a small, fixed number of search directions to build the Hessian approximation. Set the history length with the constructor. If not provided, a default value of 3 is used. More history can help on curved valleys or ill-conditioned problems. Explicit histories of 10 or 20 are useful choices to compare, at the cost of more storage and work per iteration. Check callback counts and the original gradient when comparing them.

The function for which an extremum is to be found is called the objective function. It is accessed through the ObjectiveFunction property, which is a delegate of type Func<T, TResult>. The delegate takes a Vector<T> argument and returns a real number.

The following code demonstrates how to write methods to implement objective functions. The objective function implemented here is the Rosenbrock function,

f(x₁,x₂) = (1 - x₁)² + 105(x₂-x₁²)²

which is a famous test function for optimization.

double Rosenbrock(Vector<double> x)
{
    double a = 1 - x[0];
    double b = (x[1] - x[0] * x[0]);
    return a * a + 105 * b * b;
}

Function Rosenbrock(x As Vector(Of Double)) As Double
    Dim a As Double = 1 - x(0)
    Dim b As Double = (x(1) - x(0) * x(0))
    Return a * a + 105 * b * b
End Function

No code example is currently available or this language may not be supported.

let Rosenbrock (x : Vector<float>) =
    let a = 1.0 - x[0]
    let b = (x[1] - x[0] * x[0])
    a * a + 105.0 * b * b

Some algorithms require the gradient (vector of partial derivatives) of the objective function. There are two ways to supply this function. The first is through the GradientFunction property. This property is of type Func<T, TResult>. It is a method that takes a Vector<T> argument and also returns a Vector<T>. The drawback of using this property is that a new Vector<T> instance is created on every call.

The following code demonstrates how to write a method that calculates the gradient of an objective function. Once again, we use the Rosenbrock function as our example:

Vector<double> RosenbrockGradient(Vector<double> x)
{
    Vector<double> f = Vector.Zeros<double>(2);
    double a = 1 - x[0];
    double b = (x[1] - x[0] * x[0]);
    f[0] = -2 * a - 420 * x[0] * b;
    f[1] = 210 * b;
    return f;
}

Function RosenbrockGradient(x As Vector(Of Double)) As Vector(Of Double)
    Dim f = Vector.Zeros(Of Double)(2)
    Dim a As Double = 1 - x(0)
    Dim b As Double = (x(1) - x(0) * x(0))
    f(0) = -2 * a - 420 * x(0) * b
    f(1) = 210 * b
    Return f
End Function

No code example is currently available or this language may not be supported.

let RosenbrockGradient (x : Vector<float>) =
    let f = Vector.Zeros<float>(2)
    let a = 1.0 - x[0]
    let b = (x[1] - x[0] * x[0])
    f[0] <- -2.0 * a - 420.0 * x[0] * b
    f[1] <- 210.0 * b
    f

Alternatively, the gradient can be specified using the FastGradientFunction property. This property is of type Func<T1, T2, TResult>. It is a method that takes two Vector<T> arguments. The first is the input vector. The second is the Vector<T> that is to hold the result. The return value is a reference to this result vector.

The following code demonstrates how to write a method that calculates the gradient of an objective function as a Func<T1, T2, TResult> delegate.

Vector<double> RosenbrockFastGradient(Vector<double> x, Vector<double> f)
{
    if (f == null)
        f = Vector.Zeros<double>(2);
    double a = 1 - x[0];
    double b = (x[1] - x[0] * x[0]);
    f[0] = -2 * a - 420 * x[0] * b;
    f[1] = 210 * b;
    return f;
}

Function RosenbrockFastGradient(x As Vector(Of Double), f As Vector(Of Double)) As Vector(Of Double)
    If (f Is Nothing) Then f = Vector.Zeros(Of Double)(2)
    Dim a As Double = 1 - x(0)
    Dim b As Double = (x(1) - x(0) * x(0))
    f(0) = -2 * a - 420 * x(0) * b
    f(1) = 210 * b
    Return f
End Function

No code example is currently available or this language may not be supported.

let RosenbrockFastGradient (x : Vector<float>) (f : Vector<float>) =
    let a = 1.0 - x[0]
    let b = (x[1] - x[0] * x[0])
    f[0] <- -2.0 * a - 420.0 * x[0] * b
    f[1] <- 210.0 * b
    f

If the gradient function is not available, or if it is very expensive to calculate, it may be preferable to either approximate it numerically, or use a method that does not require the gradient function.

The next example puts this all together. It creates a BFGS optimizer and sets it up to minimize the Rosenbrock function:

QuasiNewtonOptimizer bfgs =
    new QuasiNewtonOptimizer(QuasiNewtonMethod.Bfgs);
bfgs.ObjectiveFunction = Rosenbrock;
bfgs.FastGradientFunction = RosenbrockFastGradient;

Dim bfgs As QuasiNewtonOptimizer =
    New QuasiNewtonOptimizer(QuasiNewtonMethod.Bfgs)
bfgs.ObjectiveFunction = AddressOf Rosenbrock
bfgs.FastGradientFunction = AddressOf RosenbrockFastGradient

No code example is currently available or this language may not be supported.

let bfgs = QuasiNewtonOptimizer(QuasiNewtonMethod.Bfgs)
bfgs.ObjectiveFunction <- Func<_,_>(Rosenbrock)
bfgs.FastGradientFunction <- Func<_,_,_>(RosenbrockFastGradient)

When using .NET 4.0, it is possible to have the gradient computed using Automatic Differentiation. To do so, specify the objective function by setting the SymbolicObjectiveFunction property to a lambda expression that represents the objective function. The lambda expression should take a Vector<T> as its input.

Optimization algorithms in multiple dimensions commonly have two or three criteria that may each signal that an extremum has been found.

The ValueTest property returns a SimpleConvergenceTest object that represents the convergence test based on the value of the objective function. The test is successful if the change in the value of the objective function is less than the tolerance. The test returns a value of Divergent if the value of the objective function is infinite, and BadFunction when the value is NaN.

There is some danger in using this test, since some algorithms may not return a new best estimate for the extremum on each iteration. For this reason, the test is not active by default. To make it active, its Enabled property must be set to true.

The SolutionTest property returns a VectorConvergenceTest<T> object that represents the convergence test based on the estimated extremum. The test is successful if the change in the approximate extremum is less than the tolerance. The test returns a value of Divergent if one of the elements of the solution is infinite, and BadFunction when one of the values is NaN.

The VectorConvergenceTest has many options to specify the details of the convergence test. By default, the error is calculated as the maximum of the errors in each of the elements of the change of the extremum from the previous iteration compared to the corresponding elements of the extremum.

Methods that use gradients have a third test available. The GradientTest property returns a VectorConvergenceTest<T> based on the gradient of the objective function. The test succeeds if the norm of the gradient vector is less than the tolerance. This test is inactive for algorithms that do not use gradients.

For dense BFGS/DFP, unbounded L-BFGS and first-derivative conjugate gradients, the untouched defaults require the gradient test to pass. A small relative position change alone does not establish convergence, especially when coordinates are large. If rounding prevents any movement while the gradient remains too large, the status is RoundOffError. Explicit position-test settings, a disabled gradient test or a custom convergence collection preserve the selected stopping policy. A small gradient is a local stopping condition; it does not guarantee global optimality or accurate parameters.

Specific algorithms may supply their own convergence tests, or may modify the exact meaning of the tests discussed here. For example, the Nelder-Mead method uses slight variations of the value and solution tests that are more convenient for that particular algorithm.

The example below illustrates some of the options mentioned here. It sets the tolerance of the gradient test to 100 times the machine precision, which will succeed if the sum of the absolute values is less than the tolerance. It also explicitly specifies that the value test is not to be used.

bfgs.ValueTest.Enabled = false;
bfgs.GradientTest.Norm = VectorConvergenceNorm.SumOfAbsoluteValues;
bfgs.GradientTest.Tolerance = 100 * MachineConstants.Epsilon;

bfgs.ValueTest.Enabled = False
bfgs.GradientTest.Norm = VectorConvergenceNorm.SumOfAbsoluteValues
bfgs.GradientTest.Tolerance = 100 * MachineConstants.Epsilon

No code example is currently available or this language may not be supported.

bfgs.ValueTest.Enabled <- false
bfgs.GradientTest.Norm <- VectorConvergenceNorm.SumOfAbsoluteValues
bfgs.GradientTest.Tolerance <- 100.0 * MachineConstants.Epsilon

Full CMA-ES can spend substantial time factoring its covariance matrix. CovarianceRefreshInterval sets a positive maximum number of generations between factor refreshes. The default, 1, refreshes every generation. Larger values reuse sampling and whitening factors while updating covariance and paths on every generation. Intervals such as 2 or 5 can reduce decomposition work and allocation, but can require more objective evaluations or change final accuracy. Compare results on your objective before choosing one.

Active covariance updates refresh every generation to keep their covariance metric current. Diagonal covariance also refreshes every generation because that takes linear work. These policies can refresh sooner than a supplied maximum. Each new solve starts with fresh factors; covariance is checked for finite values even on generations that reuse them.

Nelder-Mead's existing factor properties also allow dimension-dependent settings. For a problem with n variables, n at least 2, keep ReflectionFactor at -1 and consider ExpansionFactor = 1 + 2/n, ContractionFactor = 3/4 - 1/(2n), and ShrinkageFactor = 1 - 1/n. These adaptive factors can reduce work in higher dimensions; they are not uniformly better on smaller problems. The classic factors remain the defaults. Geometric and fitness stopping still do not certify stationarity.

The differential-evolution, particle-swarm and CMA-ES optimizers use a population of points. Their default convergence collection requires all enabled tests to pass: small changes in the best value for successive iterations, a small range of population values, and a small geometric spread. The value-change and fitness-range tests use objective units and absolute tolerances. The geometric test uses the Euclidean norm of coordinate widths across every population member. Adding a constant to an objective no longer changes stopping through division by the objective mean.

Floating-point values still limit what a search can distinguish. If a large constant hides small improvements through rounding, even the selection of candidates can change. Adjust the objective's scale when those differences matter. Tolerances depend on the units of the objective and variables. Tighter geometry can improve accuracy, but these tests do not certify stationarity or global optimality. The diagonal CMA-ES variant does not learn correlations and can need more work on rotated valleys.

Use the population optimizer's ValueTest, PopulationConvergenceTest and DiversityTest to configure stopping. The default best-value test requires 20 successive passes. Differential evolution and particle swarm also expose StagnationLimit for this count. A changed quantifier, a disabled test or a supplied test evaluator retains the caller's selected policy. Disabling every test leaves the iteration and evaluation limits in control.

Simulated annealing keeps its current walk separate from the best point found. Accepted uphill moves become the starting point for the next proposal. Extremum, BestSolution and BestValue report the best evaluated point, including on budget exits; values use the original objective sign. Local step and value tests are disabled by default. StagnationLimit counts proposals without a new best value, and cooling occurs after each complete temperature level. Temperature or stagnation termination is a search policy, not a numerical certificate. Callers can explicitly enable other tests.

Many algorithms for Multidimensional optimization work by successively solving a one-dimensional optimization problem along a search direction, and altering the search direction in each iteration. These algorithms all inherit from DirectionalOptimizer. Although general one-dimensional optimization algorithms may be used, they usually don't work well. They are either wasteful because they are too precise, or the results fail to meet certain criteria that are required to ensure convergence of the algorithm.

Specialized line search algorithms are therefore in order. Numerics.NET provides three of them:

All optimization algorithms that use line searches automatically select the appropriate line search algorithm, and set the search parameters to appropriate values. The unit line search is useful for experimental purposes, or when the objective function is known to be well behaved so that a line search is not necessary. You can use the LineSearch property to access the line search object or assign a different search before solving. Both limited-memory BFGS methods use a supplied search. With variable bounds, trial points must stay in the feasible step range: the optimizer prevents evaluations outside that range and rejects an infeasible returned step. Set the property to null to select the method's default on the next solve. The search cannot be replaced during a calculation.

Dense BFGS uses More-Thuente with curvature factor 0.8. DFP and first-derivative conjugate gradients use factor 0.1 to control curvature more tightly. Unbounded L-BFGS retains its existing More-Thuente driver and factor 0.9. A supplied search retains its settings. If it accepts a step with weak or nonfinite curvature, quasi-Newton methods discard that update and recover a descent direction when needed. Dense approximations begin at identity; unbounded L-BFGS accepts only finite positive history scales.

RestartIterations and RestartIterations let you discard accumulated curvature at regular intervals. For example, a value of 20 resets the approximation before the next direction after iterations 20, 40, and so on. The newest curvature pair is discarded along with the older information. The approximation restarts from an identity matrix. This can require more iterations, so use a reset interval when you have a reason to do so.

Zero disables scheduled resets. The existing defaults remain: int.MaxValue for the dense optimizer and zero for the bounded optimizer. Set the interval before solving; it cannot change during a calculation. Each new solve starts with a fresh approximation regardless of the interval.

You must specify an initial guess for the solution by setting the InitialGuess property to a suitable Vector<T> value. Any vector type may be used.

The FindExtremum method finds the extremum. It returns a DenseVector<T> that contains the best approximation to the extremum. This value is also available through the Extremum

bfgs.InitialGuess = Vector.Create(-1.2, 1.0);
bfgs.ExtremumType = ExtremumType.Minimum;
bfgs.FindExtremum();

bfgs.InitialGuess = Vector.Create(-1.2, 1.0)
bfgs.ExtremumType = ExtremumType.Minimum
bfgs.FindExtremum()

No code example is currently available or this language may not be supported.

bfgs.InitialGuess <- Vector.Create(-1.2, 1.0)
bfgs.ExtremumType <- ExtremumType.Minimum
bfgs.FindExtremum() |> ignore

After the FindExtremum method returns, a number of properties give more details about the algorithm. The Status property returns an AlgorithmStatus value that indicates how the algorithm terminated. The possible values are listed below:

Each of the ConvergenceTest objects discussed earlier exposes an Error property that returns the estimated error used in the test. The IterationsNeeded and EvaluationsNeeded properties return the number of iterations and the number of evaluations of the objective function, respectively. Algorithms that use gradients also have a GradientEvaluationsNeeded.

The final example prints a summary of the results of running the BFGS algorithm on the Rosenbrock function:

Console.WriteLine("BFGS Method:");
Console.WriteLine("  Solution: {0}", bfgs.Extremum);
Console.WriteLine("  Estimated error: {0}", bfgs.SolutionTest.Error);
Console.WriteLine("  Gradient residual: {0}", bfgs.GradientTest.Error);
Console.WriteLine("  # iterations: {0}", bfgs.IterationsNeeded);
Console.WriteLine("  # function evaluations: {0}",
    bfgs.EvaluationsNeeded);
Console.WriteLine("  # gradient evaluations: {0}",
    bfgs.GradientEvaluationsNeeded);

Console.WriteLine("BFGS Method:")
Console.WriteLine("  Solution: {0}", bfgs.Extremum)
Console.WriteLine("  Estimated error: {0}", bfgs.SolutionTest.Error)
Console.WriteLine("  Gradient residual: {0}", bfgs.GradientTest.Error)
Console.WriteLine("  # iterations: {0}", bfgs.IterationsNeeded)
Console.WriteLine("  # function evaluations: {0}",
    bfgs.EvaluationsNeeded)
Console.WriteLine("  # gradient evaluations: {0}",
    bfgs.GradientEvaluationsNeeded)

No code example is currently available or this language may not be supported.

printfn "BFGS Method:"
printfn $"  Solution: {bfgs.Extremum}"
printfn $"  Estimated error: {bfgs.SolutionTest.Error}"
printfn $"  Gradient residual: {bfgs.GradientTest.Error}"
printfn $"  # iterations: {bfgs.IterationsNeeded}"
printfn $"  # function evaluations: {bfgs.EvaluationsNeeded}"
printfn $"  # gradient evaluations: {bfgs.GradientEvaluationsNeeded}"

References

Fletcher, R. Practical Methods of Optimization, John Wiley & Sons, New York, 1987.

Fletcher, R. and Reeves, C.M., Function Minimization by Conjugate Gradients, Comp, J. 7, 149-154, 1964.

Moré and D. Thuente, Line Search Algorithms with Guaranteed Sufficient Decrease, ACM Transactions on Mathematical Software , Vol. 20, No. 3, pp. 286--307, 1994.

Nelder, J. A. and Mead, R. "A Simplex Method for Function Minimization." Comput. J. 7, 308-313, 1965.

Polak, E., and Ribière, G., Note sur la Convergence de Methodes de Directions Conjugees, Revue Francaise d'Informatique et de Recherche Operationnelle, 3, 35-43, 1969.

Powell, M.J.D., An Efficient Method for Finding the Minimum of a Function of Several Variables Without Calculating Derivatives, Comp. J. 7, 155-162, 1964.