optimize

High level D API for Levenberg-Marquardt Algorithm.

Computes the argmin over x of sum_i(f(x_i)^2) using the Modified 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.

  1. void optimize(scope ref LeastSquaresLM!T lm)
    void
    optimize
    (
    alias g = null
    alias tm = null
    T
    )
    (
    scope ref LeastSquaresLM!T lm
    )
    if (
    (
    is(T == float) ||
    is(T == double)
    )
    &&
    __traits(compiles, optimizeImpl!(f, g, tm, T))
    )
  2. void optimize(scope ref LeastSquaresLM!T lm, 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)

lm
Type: LeastSquaresLM!T

Levenberg-Marquardt data structure

Throws

Examples

With Jacobian

1 import mir.ndslice.allocation: slice;
2 import mir.ndslice.slice: sliced;
3 import mir.blas: nrm2;
4 
5 auto lm = LeastSquaresLM!double(2, 2);
6 lm.x[] = [100, 100];
7 lm.optimize!(
8     (x, y)
9     {
10         y[0] = x[0];
11         y[1] = 2 - x[1];
12     },
13     (x, J)
14     {
15         J[0, 0] = 1;
16         J[0, 1] = 0;
17         J[1, 0] = 0;
18         J[1, 1] = -1;
19     },
20 );
21 
22 assert(nrm2((lm.x - [0, 2].sliced).slice) < 1e-8);

Using Jacobian finite difference approximation computed using in multiple threads.

1 import mir.ndslice.allocation: slice;
2 import mir.ndslice.slice: sliced;
3 import mir.blas: nrm2;
4 import std.parallelism: taskPool;
5 
6 auto lm = LeastSquaresLM!double(2, 2);
7 lm.x[] = [-1.2, 1];
8 lm.optimize!(
9     (x, y) // Rosenbrock function
10     {
11         y[0] = 10 * (x[1] - x[0]^^2);
12         y[1] = 1 - x[0];
13     },
14 )(taskPool);
15 
16 assert(nrm2((lm.x - [1, 1].sliced).slice) < 1e-6);

Rosenbrock

1 import mir.algorithm.iteration: all;
2 import mir.ndslice.allocation: slice;
3 import mir.ndslice.slice: sliced;
4 import mir.blas: nrm2;
5 
6 auto lm = LeastSquaresLM!double(2, 2, Yes.lowerBounds, Yes.upperBounds);
7 lm.x[] = [-1.2, 1];
8 
9 alias rosenbrockRes = (x, y)
10 {
11     y[0] = 10 * (x[1] - x[0]^^2);
12     y[1] = 1 - x[0];
13 };
14 
15 alias rosenbrockJac = (x, J)
16 {
17     J[0, 0] = -20 * x[0];
18     J[0, 1] = 10;
19     J[1, 0] = -1;
20     J[1, 1] = 0;
21 };
22 
23 static class FFF
24 {
25     static auto opCall(Slice!(const(double)*) x, Slice!(double*, 2) J)
26     {
27         rosenbrockJac(x, J);
28     }
29 }
30 
31 lm.optimize!(rosenbrockRes, FFF);
32 
33 // import std.stdio;
34 
35 // writeln(lm.iterCt, " ", lm.fCalls, " ", lm.gCalls);
36 
37 assert(nrm2((lm.x - [1, 1].sliced).slice) < 1e-8);
38 
39 /////
40 
41 lm.reset;
42 lm.lower[] = [10.0, 10.0];
43 lm.upper[] = [200.0, 200.0];
44 lm.x[] = [150.0, 150.0];
45 
46 lm.optimize!(rosenbrockRes, rosenbrockJac);
47 
48 // writeln(lm.iterCt, " ", lm.fCalls, " ", lm.gCalls);
49 assert(nrm2((lm.x - [10, 10].sliced).slice) < 1e-8);
50 assert(lm.x.all!"a >= 10");
1 import mir.blas: nrm2;
2 import mir.math.common;
3 import mir.ndslice.allocation: slice;
4 import mir.ndslice.topology: linspace, map;
5 import mir.ndslice.slice: sliced;
6 import mir.random;
7 import mir.random.algorithm;
8 import mir.random.variable;
9 import std.parallelism: taskPool;
10 
11 alias model = (x, p) => p[0] * map!exp(-x * p[1]);
12 
13 auto p = [1.0, 2.0];
14 
15 auto xdata = [20].linspace([0.0, 10.0]);
16 auto rng = Random(12345);
17 auto ydata = slice(model(xdata, p) + 0.01 * rng.randomSlice(normalVar, xdata.shape));
18 
19 auto lm = LeastSquaresLM!double(xdata.length, 2);
20 lm.x[] = [0.5, 0.5];
21 
22 lm.optimize!((p, y) => y[] = model(xdata, p) - ydata)();
23 
24 assert((lm.x - [1.0, 2.0].sliced).slice.nrm2 < 0.05);
1 import mir.algorithm.iteration: all;
2 import mir.ndslice.allocation: slice;
3 import mir.ndslice.topology: map, repeat, iota;
4 import mir.ndslice.slice: sliced;
5 import mir.random;
6 import mir.random.variable;
7 import mir.random.algorithm;
8 import mir.math.common;
9 
10 alias model = (x, p) => p[0] * map!exp(-x / p[1]) + p[2];
11 
12 auto xdata = iota([100], 1);
13 auto rng = Random(12345);
14 auto ydata = slice(model(xdata, [10.0, 10.0, 10.0]) + 0.1 * rng.randomSlice(normalVar, xdata.shape));
15 
16 auto lm = LeastSquaresLM!double(xdata.length, 3, [5.0, 11.0, 5.0], double.infinity.repeat(3).slice.field);
17 
18 lm.x[] = [15.0, 15.0, 15.0];
19 lm.optimize!((p, y) => 
20     y[] = model(xdata, p) - ydata);
21 
22 assert(all!"a >= b"(lm.x, lm.lower));
23 
24 // import std.stdio;
25 
26 // writeln(lm.x);
27 // writeln(lm.iterCt, " ", lm.fCalls, " ", lm.gCalls);
28 
29 lm.reset;
30 lm.x[] = [5.0, 5.0, 5.0];
31 lm.upper[] = [15.0, 9.0, 15.0];
32 lm.optimize!((p, y) => y[] = model(xdata, p) - ydata);
33 
34 assert(all!"a <= b"(lm.x, lm.upper));
35 
36 // writeln(lm.x);
37 // writeln(lm.iterCt, " ", lm.fCalls, " ", lm.gCalls);
1 import mir.blas: nrm2;
2 import mir.math.common: sqrt;
3 import mir.ndslice.allocation: slice;
4 import mir.ndslice.slice: sliced;
5 
6 auto lm = LeastSquaresLM!double(1, 2, [-0.5, -0.5], [0.5, 0.5]);
7 lm.x[] = [0.001, 0.0001];
8 lm.optimize!(
9     (x, y)
10     {
11         y[0] = sqrt(1 - (x[0] ^^ 2 + x[1] ^^ 2));
12     },
13 );
14 
15 assert(nrm2((lm.x - lm.upper).slice) < 1e-8);

See Also

Meta