optimize

High level D API for Levenberg-Marquardt Algorithm.

Computes the argmin over x of sum_i(f(x_i)^2) using the Levenberg-Marquardt algorithm, and an estimate of the Jacobian of f at x.

The function f should take an input vector of length n, and fill an output vector of length m.

The function g is the Jacobian of f, and should fill a row-major m x n matrix.

LeastSquaresResult!T optimize(LeastSquaresSettings!T settings, size_t m, Slice!(T*) x, Slice!(const(T)*) l, Slice!(const(T)*) u)
LeastSquaresResult!T
optimize
(
alias f
alias g = null
alias tm = null
T
)
(
scope ref LeastSquaresSettings!T settings
,
size_t m
,
Slice!(T*) x
,
Slice!(const(T)*) l
,
Slice!(const(T)*) u
)
if (
(
is(T == float) ||
is(T == double)
)
)
LeastSquaresResult!T optimize(LeastSquaresSettings!T settings, size_t m, Slice!(T*) x, Slice!(const(T)*) l, Slice!(const(T)*) u, TaskPool taskPool)

Parameters

f: n -> m function
g: m × n Jacobian (optional)
tm: thread manager for finite difference jacobian approximation in case of g is null (optional)
settings LeastSquaresSettings!T: Levenberg-Marquardt data structure

Throws

LeastSquaresException

Examples

With Jacobian

import mir.ndslice.allocation: slice;
import mir.ndslice.slice: sliced;
import mir.blas: nrm2;

LeastSquaresSettings!double settings;
auto x = [100.0, 100].sliced;
auto l = x.shape.slice(-double.infinity);
auto u = x.shape.slice(+double.infinity);
optimize!(
    (x, y)
    {
        y[0] = x[0];
        y[1] = 2 - x[1];
    },
    (x, J)
    {
        J[0, 0] = 1;
        J[0, 1] = 0;
        J[1, 0] = 0;
        J[1, 1] = -1;
    },
)(settings, 2, x, l, u);

assert(nrm2((x - [0, 2].sliced).slice) < 1e-8);

Using Jacobian finite difference approximation computed using in multiple threads.

import mir.ndslice.allocation: slice;
import mir.ndslice.slice: sliced;
import mir.blas: nrm2;
import std.parallelism: taskPool;

LeastSquaresSettings!double settings;
auto x = [-1.2, 1].sliced;
auto l = x.shape.slice(-double.infinity);
auto u = x.shape.slice(+double.infinity);
settings.optimize!(
    (x, y) // Rosenbrock function
    {
        y[0] = 10 * (x[1] - x[0]^^2);
        y[1] = 1 - x[0];
    },
)(2, x, l, u, taskPool);

// import std.stdio;
// writeln(settings);
// writeln(x);

assert(nrm2((x - [1, 1].sliced).slice) < 1e-6);

Rosenbrock

import mir.algorithm.iteration: all;
import mir.ndslice.allocation: slice;
import mir.ndslice.slice: Slice, sliced;
import mir.blas: nrm2;

LeastSquaresSettings!double settings;
auto x = [-1.2, 1].sliced;
auto l = x.shape.slice(-double.infinity);
auto u = x.shape.slice(+double.infinity);

alias rosenbrockRes = (x, y)
{
    y[0] = 10 * (x[1] - x[0]^^2);
    y[1] = 1 - x[0];
};

alias rosenbrockJac = (x, J)
{
    J[0, 0] = -20 * x[0];
    J[0, 1] = 10;
    J[1, 0] = -1;
    J[1, 1] = 0;
};

static class FFF
{
    static auto opCall(Slice!(const(double)*) x, Slice!(double*, 2) J)
    {
        rosenbrockJac(x, J);
    }
}

settings.optimize!(rosenbrockRes, FFF)(2, x, l, u);

// import std.stdio;

// writeln(settings.iterations, " ", settings.fCalls, " ", settings.gCalls, " x = ", x);

assert(nrm2((x - [1, 1].sliced).slice) < 1e-8);

/////

settings = settings.init;
x[] = [150.0, 150.0];
l[] = [10.0, 10.0];
u[] = [200.0, 200.0];

settings.optimize!(rosenbrockRes, rosenbrockJac)(2, x, l, u);

// writeln(settings.iterations, " ", settings.fCalls, " ", settings.gCalls, " ", x);
assert(nrm2((x - [10, 100].sliced).slice) < 1e-5);
assert(x.all!"a >= 10");

import mir.blas: nrm2;
import mir.math.common;
import mir.ndslice.allocation: slice;
import mir.ndslice.topology: linspace, map;
import mir.ndslice.slice: sliced;
import mir.random;
import mir.random.algorithm;
import mir.random.variable;
import std.parallelism: taskPool;

alias model = (x, p) => p[0] * map!exp(-x * p[1]);

auto p = [1.0, 2.0];

auto xdata = [20].linspace([0.0, 10.0]);
auto rng = Random(12345);
auto ydata = slice(model(xdata, p) + 0.01 * rng.randomSlice(normalVar, xdata.shape));

auto x = [0.5, 0.5].sliced;
auto l = x.shape.slice(-double.infinity);
auto u = x.shape.slice(+double.infinity);

LeastSquaresSettings!double settings;
settings.optimize!((p, y) => y[] = model(xdata, p) - ydata)(ydata.length, x, l, u);

assert((x - [1.0, 2.0].sliced).slice.nrm2 < 0.05);

import mir.algorithm.iteration: all;
import mir.ndslice.allocation: slice;
import mir.ndslice.topology: map, repeat, iota;
import mir.ndslice.slice: sliced;
import mir.random;
import mir.random.variable;
import mir.random.algorithm;
import mir.math.common;

alias model = (x, p) => p[0] * map!exp(-x / p[1]) + p[2];

auto xdata = iota([100], 1);
auto rng = Random(12345);
auto ydata = slice(model(xdata, [10.0, 10.0, 10.0]) + 0.1 * rng.randomSlice(normalVar, xdata.shape));

LeastSquaresSettings!double settings;

auto x = [15.0, 15.0, 15.0].sliced;
auto l = [5.0, 11.0, 5.0].sliced;
auto u = x.shape.slice(+double.infinity);

settings.optimize!((p, y) => y[] = model(xdata, p) - ydata)
    (ydata.length, x, l, u);

assert(all!"a >= b"(x, l));

// import std.stdio;

// writeln(x);
// writeln(settings.iterations, " ", settings.fCalls, " ", settings.gCalls);

settings = settings.init;
x[] = [5.0, 5.0, 5.0];
l[] = -double.infinity;
u[] = [15.0, 9.0, 15.0];
settings.optimize!((p, y) => y[] = model(xdata, p) - ydata)
    (ydata.length, x, l , u);

assert(x.all!"a <= b"(u));

// writeln(x);
// writeln(settings.iterations, " ", settings.fCalls, " ", settings.gCalls);

import mir.blas: nrm2;
import mir.math.common: sqrt;
import mir.ndslice.allocation: slice;
import mir.ndslice.slice: sliced;

LeastSquaresSettings!double settings;
auto x = [0.001, 0.0001].sliced;
auto l = [-0.5, -0.5].sliced;
auto u = [0.5, 0.5].sliced;
settings.optimize!(
    (x, y)
    {
        y[0] = sqrt(1 - (x[0] ^^ 2 + x[1] ^^ 2));
    },
)(1, x, l, u);

assert(nrm2((x - u).slice) < 1e-8);

optimize

Parameters

Throws

Examples

See Also

Meta

Source