Find a Maximum or Minimum

Finding a maximum or minimum is an extremely common practical scenario. Here's some code that finds the the maximum in the first period of the sine function:

using System;
using Meta.Numerics;
using Meta.Numerics.Analysis;

Interval range = Interval.FromEndpoints(0.0, 6.28);
Extremum min = FunctionMath.FindMaximum(x => Math.Sin(x), range);
Console.WriteLine($"Found maximum of {min.Value} at {min.Location}.");
Console.WriteLine($"Required {min.EvaluationCount} evaluations.");

To find the minimum instead of the maximum use FunctionMath.FindMinimum.

Why is the location not accurate to full precision?

The location should returned in the example above should be very close to pi/2 = 1.5707963267948966, but will likely still off in the last 8 or so digits. Isn't Meta.Numerics supposed to achieve nearly full (~16-digit) precision?

It is, but for minimum-finding we run up against a fundamental limitation of floating-point. A typical minimum is approximately quadratic, so near the minimum δy ~ (δx)^2. To detect that we have moved off the minimum, we need δy ~ ε, where ε is the floating point precision (about 1.0E-16). But that implies δx ~ √ε (about 1.0E-8). So, for typical minima, no algorithm using standard floating point arithmetic is going to do better than about half of full precision.

What about multi-dimensional minima?

In real-world practical scenarios, you often want to find the minimum of a function that depends on several (or several dozen) variables. Fortunately, Meta.Numerics can do that, too. Here is some code that finds the minimum of the well-known Rosenbrock test function:

using Meta.Numerics.Matrices;

MultiExtremum rosenbrock = MultiFunctionMath.FindLocalMinimum(
    x => MoreMath.Sqr(2.0 - x[0]) + 100.0 * MoreMath.Sqr(x[1] - x[0] * x[0]),
    new ColumnVector(0.0, 0.0)
);
ColumnVector xm = rosenbrock.Location;
Console.WriteLine($"Found minimum of {rosenbrock.Value} at ({xm[0]},{xm[1]}).");
Console.WriteLine($"Required {rosenbrock.EvaluationCount} evaluations.");

Notice that you are required to supply a starting point for the multi-dimensional search.

What if I don't know a good starting point?

The textbook answer is that you should always find a good starting point, but real life doesn't give you textbook problems.

Meta.Numerics' answer is that you need to at least tell us the region over which you are willing to search. Here is some code that finds the global minimum of the Levi function, a test function designed to have many local minima in search region:

Func<IReadOnlyList<double>, double> leviFunction = z => {
    double x = z[0];
    double y = z[1];
    return(
        MoreMath.Sqr(MoreMath.Sin(3.0 * Math.PI * x)) +
        MoreMath.Sqr(x - 1.0) * (1.0 + MoreMath.Sqr(MoreMath.Sin(3.0 * Math.PI * y))) +
        MoreMath.Sqr(y - 1.0) * (1.0 + MoreMath.Sqr(MoreMath.Sin(2.0 * Math.PI * y)))
    );
};

Interval[] leviRegion = new Interval[] {
    Interval.FromEndpoints(-10.0, +10.0),
    Interval.FromEndpoints(-10.0, +10.0)
};

MultiExtremum levi = MultiFunctionMath.FindGlobalMinimum(leviFunction, leviRegion);
ColumnVector zm = levi.Location;
Console.WriteLine($"Found minimum of {levi.Value} at ({zm[0]},{zm[1]}).");
Console.WriteLine($"Required {levi.EvaluationCount} evaluations.");

The FindGlobalMinimum function does a very different sort of searching than FindLocalMinimum. The FindLocalMinimum function models the behavior of your function near the current best-guess location and then uses the minimum of that model as the next best-guess. The FindGlobalMinimum maintains a collection of pretty-good minima within the search region, then perturbs and hybridizes those minima to try to improve its collection.

FindGlobalMinimum requires many more function evaluations to converge on a minimum (and by default targets a much lower precision so as not to take inordinately long), but it is much less likely to get caught in a non-global minimum than FindLocalMinimum.

How do I set termination criteria and monitor progress?

Like other Meta.Numerics.Analysis APIs, our optimization APIs can accept a settings object that define an evaluation budget, a target precision, and a listener that can monitor progress. Note that, for multi-dimensional optimization, the target precision refers to the y-value: the optimization terminates when the minimum y-value stops changing by more than the tolerance that the precision settings define.

Here is a cool example that puts together all our machinery to solve the famous Thompson Problem: how will n identical charges free to move on the unit sphere arrange themselves so as to minimize to total electrostatic potential?

// Select a number of charges
// The dimension of the minimization problem is 2n
int n = 6;

// Define a function that takes 2n polar coordinates in the form
// phi_1, theta_1, phi_2, theta_2, ..., phi_n, theta_n, and computes
// the sum of the potential energy 1/d for all pairs.
Func<IReadOnlyList<double>, double> thompson = (IReadOnlyList<double> v) => {
    double e = 0.0;
    // Iterate over all distinct pairs of points
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < i; j++) {
            // Compute the chord length between charges i and j
            double dx = Math.Cos(v[2 * j]) * Math.Cos(v[2 * j + 1]) - Math.Cos(v[2 * i]) * Math.Cos(v[2 * i + 1]);
            double dy = Math.Cos(v[2 * j]) * Math.Sin(v[2 * j + 1]) - Math.Cos(v[2 * i]) * Math.Sin(v[2 * i + 1]);
            double dz = Math.Sin(v[2 * j]) - Math.Sin(v[2 * i]);
            double d = Math.Sqrt(dx * dx + dy * dy + dz * dz);
            e += 1.0 / d;
        }
    }
    return (e);
};

// Limit all polar coordinates to their standard ranges.
Interval[] box = new Interval[2 * n];
for (int i = 0; i < n; i++) {
    box[2 * i] = Interval.FromEndpoints(-Math.PI, Math.PI);
    box[2 * i + 1] = Interval.FromEndpoints(-Math.PI / 2.0, Math.PI / 2.0);
}

// Use settings to monitor progress toward a rough solution.
MultiExtremumSettings roughSettings = new MultiExtremumSettings() {
    RelativePrecision = 0.05, AbsolutePrecision = 0.0,
    Listener = r => {
        Console.WriteLine($"Minimum {r.Value} after {r.EvaluationCount} evaluations.");
    }
};
MultiExtremum rough = MultiFunctionMath.FindGlobalMinimum(thompson, box, roughSettings);

// Use settings to monitor progress toward a refined solution.
MultiExtremumSettings refinedSettings = new MultiExtremumSettings() {
    RelativePrecision = 1.0E-5, AbsolutePrecision = 0.0,
    Listener = r => {
        Console.WriteLine($"Minimum {r.Value} after {r.EvaluationCount} evaluations.");
    }
};
MultiExtremum refined = MultiFunctionMath.FindLocalMinimum(thompson, rough.Location, refinedSettings);

Console.WriteLine($"Minimum potential energy {refined.Value}.");
Console.WriteLine($"Required {rough.EvaluationCount} + {refined.EvaluationCount} evaluations.");

Notice that we have used the very useful trick of first calling FindGlobalMinimum with a very low precision target to get close to the global minimum, then calling FindLocalMinimum with a high precision target and the minimum found by FindGlobalMinimum as a starting point to vastly improve the accuracy of the minimum.