mir.optim.least_squares

enum LeastSquaresStatus: int;

maxIterations

Maximum number of iterations reached

furtherImprovement

The algorithm cann't improve the solution

xConverged

Stationary values

gConverged

Stationary gradient

fConverged

Good (small) residual

badBounds

badGuess

badMinStepQuality

badGoodStepQuality

badStepQuality

badLambdaParams

numericError

template LeastSquaresFunction(T)

template LeastSquaresJacobian(T)

Delegates for low level D API.

template LeastSquaresFunctionBetterC(T)

Delegates for low level C API.

template LeastSquaresJacobianBetterC(T)

struct LeastSquaresSettings(T) if (is(T == double) || is(T == float));

Least-Squares iteration settings.

uint maxIterations;

Maximum number of iterations

uint maxAge;

Maximum jacobian model age (0 for default selection)

T jacobianEpsilon;

epsilon for finite difference Jacobian approximation

T absTolerance;

Absolute tolerance for step size, L2 norm

T relTolerance;

Relative tolerance for step size, L2 norm

T gradTolerance;

Absolute tolerance for gradient, L-inf norm

T maxGoodResidual;

The algorithm stops iteration when the residual value is less or equal to maxGoodResidual.

T maxStep;

maximum norm of iteration step

T maxLambda;

minimum trust region radius

T minLambda;

maximum trust region radius

T minStepQuality;

for steps below this quality, the trust region is shrinked

T goodStepQuality;

for steps above this quality, the trust region is expanded

T lambdaIncrease;

lambda is multiplied by this factor after step below min quality

T lambdaDecrease;

lambda is multiplied by this factor after good quality steps

BoxQPSettings!T qpSettings;

Bound constrained convex quadratic problem settings

struct LeastSquaresResult(T) if (is(T == double) || is(T == float));

Least-Squares results.

LeastSquaresStatus status;

Computation status

uint iterations;

Successful step count

uint fCalls;

Number of the function calls

uint gCalls;

Number of the Jacobian calls

T residual;

Final residual

T lambda;

LMA variable for (inverse of) initial trust region radius

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

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.

Throws:

LeastSquaresException

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)
LeastSquaresSettings!T settings	Levenberg-Marquardt data structure
TaskPool taskPool	task Pool with .parallel method for finite difference jacobian approximation in case of g is null (optional)

See Also:

optimizeLeastSquares

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);

Examples:

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);

Examples:

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");

Examples:

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);

Examples:

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);

Examples:

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);

LeastSquaresResult!T optimizeLeastSquares(alias f, alias g = null, alias tm = null, T)(ref scope LeastSquaresSettings!T settings, size_t m, Slice!(T*) x, Slice!(const(T)*) l, Slice!(const(T)*) u);

High level nothtow D API for Levenberg-Marquardt Algorithm.

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

Returns:

optimization status.

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)
LeastSquaresSettings!T settings	Levenberg-Marquardt data structure
size_t m	length (dimension) of y = f(x)
Slice!(T*) x	initial (in) and final (out) X value
Slice!(const(T)*) l	lower X bound
Slice!(const(T)*) u	upper X bound

See Also:

optimize

pure nothrow @nogc @safe string leastSquaresStatusString(LeastSquaresStatus st);

Status string for low (extern) and middle (nothrow) levels D API.

Parameters:

LeastSquaresStatus st

optimization status

Returns:

description for

alias LeastSquaresTask = void delegate(uint totalThreads, uint threadId, uint i) pure nothrow @nogc @safe;

alias LeastSquaresTaskBetterC = extern (C) void function(ref scope const(void delegate(uint, uint, uint) pure nothrow @nogc @safe), uint totalThreads, uint threadId, uint i) pure nothrow @nogc @safe;

alias LeastSquaresThreadManager = void delegate(uint count, scope void delegate(uint totalThreads, uint threadId, uint i) pure nothrow @nogc @safe task) pure nothrow @nogc @safe;

Thread manager delegate type for low level extern(D) API.

pure nothrow @nogc @trusted LeastSquaresResult!double optimizeLeastSquaresD(ref scope LeastSquaresSettings!double settings, size_t m, Slice!(double*) x, Slice!(const(double)*) l, Slice!(const(double)*) u, Slice!(double*) work, Slice!(lapackint*) iwork, scope LeastSquaresFunction!double f, scope LeastSquaresJacobian!double g = null, scope LeastSquaresThreadManager tm = null);

pure nothrow @nogc @trusted LeastSquaresResult!float optimizeLeastSquaresS(ref scope LeastSquaresSettings!float settings, size_t m, Slice!(float*) x, Slice!(const(float)*) l, Slice!(const(float)*) u, Slice!(float*) work, Slice!(lapackint*) iwork, scope LeastSquaresFunction!float f, scope LeastSquaresJacobian!float g = null, scope LeastSquaresThreadManager tm = null);

template optimizeLeastSquares(T : double)

template optimizeLeastSquares(T : float)

Low level extern(D) instatiation.

Parameters:

LeastSquaresSettings!double settings	Levenberg-Marquardt data structure
size_t m	length (dimension) of y = f(x)
Slice!(double*) x	initial (in) and final (out) X value
Slice!(const(double)*) l	lower X bound
Slice!(const(double)*) u	upper X bound
LeastSquaresFunction!double f	n -> m function
LeastSquaresJacobian!double g	m × n Jacobian (optional)
LeastSquaresThreadManager tm	thread manager for finite difference jacobian approximation in case of g is null (optional)
Slice!(double*) work	floating point workspace length of at least mir_least_squares_work_length
Slice!(lapackint*) iwork	floating point workspace length of at least mir_least_squares_iwork_length

pure nothrow @nogc @safe size_t mir_least_squares_work_length(size_t m, size_t n);

pure nothrow @nogc @safe size_t mir_least_squares_iwork_length(size_t m, size_t n);

pure nothrow @nogc @trusted immutable(char)* mir_least_squares_status_string(LeastSquaresStatus st);

Status string for extern(C) API.

Parameters:

LeastSquaresStatus st

optimization status

Returns:

description for

alias LeastSquaresThreadManagerBetterC = extern (C) void function(scope void* context, uint count, ref scope const(void delegate(uint, uint, uint) pure nothrow @nogc @safe) taskContext, scope extern (C) void function(ref scope const(void delegate(uint, uint, uint) pure nothrow @nogc @safe), uint totalThreads, uint threadId, uint i) pure nothrow @nogc @safe task) pure nothrow @nogc @system;

Thread manager function type for low level extern(C) API.

pure nothrow @nogc @system LeastSquaresResult!double mir_optimize_least_squares_d(ref scope LeastSquaresSettings!double settings, size_t m, size_t n, double* x, const(double)* l, const(double)* u, Slice!(double*) work, Slice!(lapackint*) iwork, scope void* fContext, scope LeastSquaresFunctionBetterC!double f, scope void* gContext = null, scope LeastSquaresJacobianBetterC!double g = null, scope void* tmContext = null, scope LeastSquaresThreadManagerBetterC tm = null);

pure nothrow @nogc @system LeastSquaresResult!float mir_optimize_least_squares_s(ref scope LeastSquaresSettings!float settings, size_t m, size_t n, float* x, const(float)* l, const(float)* u, Slice!(float*) work, Slice!(lapackint*) iwork, scope void* fContext, scope LeastSquaresFunctionBetterC!float f, scope void* gContext = null, scope LeastSquaresJacobianBetterC!float g = null, scope void* tmContext = null, scope LeastSquaresThreadManagerBetterC tm = null);

template mir_optimize_least_squares(T : double)

template mir_optimize_least_squares(T : float)

Low level extern(C) wrapper instatiation.

Parameters:

LeastSquaresSettings!double settings	Levenberg-Marquardt data structure
void* fContext	context for the function
LeastSquaresFunctionBetterC!double f	n -> m function
void* gContext	context for the Jacobian (optional)
LeastSquaresJacobianBetterC!double g	m × n Jacobian (optional)
LeastSquaresThreadManagerBetterC tm	thread manager for finite difference jacobian approximation in case of g is null (optional)
size_t m	length (dimension) of y = f(x)
size_t n	length (dimension) of X
double* x	initial (in) and final (out) X value
const(double)* l	lower X bound
const(double)* u	upper X bound
LeastSquaresFunctionBetterC!double f	n -> m function
void* fContext	f context
LeastSquaresJacobianBetterC!double g	m × n Jacobian (optional)
void* gContext	g context
LeastSquaresThreadManagerBetterC tm	thread manager for finite difference jacobian approximation in case of g is null (optional)
void* tmContext	tm context
Slice!(double*) work	floating point workspace length of at least mir_least_squares_work_length
Slice!(lapackint*) iwork	floating point workspace length of at least mir_least_squares_iwork_length

pure nothrow @nogc @safe void mir_least_squares_init_d(ref LeastSquaresSettings!double settings);

pure nothrow @nogc @safe void mir_least_squares_init_s(ref LeastSquaresSettings!float settings);

template mir_least_squares_init(T : double)

template mir_least_squares_init(T : float)

Initialize LM data structure with default params for iteration.

Parameters:

LeastSquaresSettings!double settings

Levenberg-Marquart data structure

pure nothrow @nogc @safe void mir_least_squares_reset_d(ref LeastSquaresSettings!double settings);

pure nothrow @nogc @safe void mir_least_squares_reset_s(ref LeastSquaresSettings!float settings);

template mir_least_squares_reset(T : double)

template mir_least_squares_reset(T : float)

Resets all counters and flags, fills x, y, upper, lower, vecors with default values.

Parameters:

LeastSquaresSettings!double settings

Levenberg-Marquart data structure