Image Component Library (ICL)
|
Utility class implementing the multidimensional Levenberg Marquardt Algorithm for non-linear Optimization. More...
#include <LevenbergMarquardtFitter.h>
Classes | |
struct | Data |
utility structure that is used in the static create_data utlity method More... | |
struct | Result |
Utility structure, that represents a fitting result. More... | |
Public Types | |
typedef DynColVector< Scalar > | Vector |
vector type | |
typedef DynColVector< Scalar > | Params |
parameter vector type | |
typedef DynMatrix< Scalar > | Matrix |
matrix type (used for input data) | |
typedef icl::utils::Function < Vector, const Params &, const Vector & > | Function |
to-be-optimized function type y = f(params, x) | |
typedef icl::utils::Function < void, const Params &, const Vector &, Vector & > | Jacobian |
jacobian of F | |
typedef icl::utils::Function < void, const Result & > | DebugCallback |
Optionally given debug callback, that is called in every iterations. | |
Public Member Functions | |
LevenbergMarquardtFitter () | |
creates a dummy (null instance) | |
LevenbergMarquardtFitter (Function f, int outputDim, const std::vector< Jacobian > &js=std::vector< Jacobian >(), Scalar initialLambda=1.e-8, int maxIterations=200, Scalar minError=1.e-6, Scalar lambdaMultiplier=10, const std::string &linSolver="svd") | |
create an instance with given parameters | |
void | setUseMultiThreading (bool enable) |
enables openmp based multithreading | |
void | init (Function f, int outputDim, const std::vector< Jacobian > &js=std::vector< Jacobian >(), Scalar initialLambda=1.e-8, int maxIterations=1000, Scalar minError=1.e-6, Scalar lambdaMultiplier=10, const std::string &linSolver="svd") |
(re)-initialization method | |
Result | fit (const Matrix &xs, const Matrix &ys, Params initParams) |
actual parameter fitting with given data and start parameters | |
void | setDebugCallback (DebugCallback dbg=default_debug_callback) |
sets a debug callback method, which is called automatically in every interation | |
Static Public Member Functions | |
static Jacobian | create_numerical_jacobian (int o, Function f, float delta=1.E-5) |
creates a single numerical Jacobian for a given function f and output dim | |
static std::vector< Jacobian > | create_numerical_jacobians (int n, Function f, float delta=1.e-5) |
creates a set of numerical jacobians for output dimension n | |
static Data | create_data (const Params &p, Function f, int xDim, int yDim, int num=1000, Scalar minX=-5, Scalar maxX=5) |
creates test data using a given function | |
static void | default_debug_callback (const Result &r) |
default debug callback that simply streams r to std::cout | |
Protected Member Functions | |
Scalar | error (const Matrix &ys, const Matrix &y_est) const |
returns the complete error | |
Protected Attributes | |
Function | f |
Function f. | |
bool | useMultiThreading |
flag whether multithreading is enabled | |
Scalar | initialLambda |
initial damping parameter lambda | |
int | maxIterations |
maximum number of iterations | |
Scalar | minError |
minimum error threshold | |
Scalar | lambdaMultiplier |
mulitplier that is used to increase/decreas lambda | |
std::string | linSolver |
linear solver that is used | |
Vector | dst |
b in Mx=b of linear system solved internally | |
Matrix | J |
Jacobian Matrix (rows are Ji) | |
Matrix | H |
Hessian Matrix (J^T J) | |
std::vector< Jacobian > | js |
output buffers | |
DebugCallback | dbg |
debug callback | |
Params | params_new |
new parameters (after update step) | |
Matrix | y_est |
current estimated outputs | |
Vector | dy |
buffer for current delta |
Utility class implementing the multidimensional Levenberg Marquardt Algorithm for non-linear Optimization.
The well known Levenberg Marquardt algorithms (LMA) is used for non-linear parameter optimization. Given a function , an intial set of parameters and a set of training data , where depends on model parameters , LMA will find a local minumum of function parameters that minimizes
The complete algorithm can be found at at wikipedia (see http://en.wikipedia.org/wiki/Levenberg-Marquardt_algorithm)
The implemented Function type, is of type icl::LevenbergMarquardtFitter::Function. It gets the parameters and input value and returns a vector output. Basically all functions that match this signature can be minimized. But actually it is quite important, that dependent on the function complexity, LMA will be slower of faster, or it will perhaps even get stuck in a local minimum. As long as the function is purely linear to the parameters (e.g)
Then LMA will find find the solution in a single step (depened on the step parameters lambda)
Due to the fact, that the LevenbergMarquardtFitter can optimized functions, with vector output, a set of jacobians js is used (of type std::vector<LevenbergMarquardtFitter::Jacobian>). js[o] defines the Jacobian for the o-th output dimension of f. Each Jacobian of is defined by a vector (here, a row-vector) of partial derivations where the ith component is . The Jacobian is always evaluated at a given position in the input data space. For the LMA, the Jacobian is evaluated in each step at each given input data point and for each output dimension. Therefore, the Jacobian J_o (for output dimension o) becomes a matrix, where each row contains the partial derivations for all parameters beta for a sigle data point \
In the code, the Jacobian is of type icl::LevenbergMarquardtFitter::Jacobian. It defines how to compute the lines of J_o. The LevenbergMarquardtFitter class supports analytic and numeric Jacobians. In case of being not able to derive an analytic Jacobian, a numerical Jacobian can automatically be generated. This will estimate the real partial derivations by an numerical approximation. However, it needs to evaluate twice, and is usually not as exact as the analytic Jacobian, which howver can lead to faster or slower convergence.
Here, a short example is given, that shows how to use the LevenbergMarquardtFitter class.
The input data space and the outputspace are 1D, and our function has 4 parameters (a,b,c,d) :
Obviously, f is linear to it's parameters. Therefore, providing an analytical Jacobian leads to a one-step optimiziation.
The analytical Jacobian for f is
#include <ICLMath/LevenbergMarquardtFitter.h> using namespace icl::utils; using namespace icl::math; typedef float real; typedef icl::math::LevenbergMarquardtFitter<real> LM; LM::Vector f(const LM::Params &p, const LM::Vector &vx){ const real x = vx[0]; return LM::Vector(1, p[0] + p[1]*x + p[2]*x*x + p[3] * x*x*x); } void j(const LM::Params &p, const LM::Vector &vx, LM::Vector &dst){ real x = vx[0]; const real Ji[] = { 1, x, x*x, x*x*x }; std::copy(Ji,Ji+4,dst.begin()); } int main(){ LM lm(f,1); // LM lm(f,1,std::vector<LM::Jacobian>(1,j)); const float p[4] = { 1,2,3,4 }; LM::Params realParams(4,p), startParams(4,1); LM::Data data = LM::create_data(realParams,f,1,1); LM::Result result = lm.fit(data.x,data.y,startParams); SHOW(result); }
In this example, the optimization takes 7 iterations if the numerical Jacobian is used. The Analytic Jacobian leads to a much better convergence and leads to a one-step fitting. If doubles instead of floats are used, the numerical Jacobian is also accurate enough for one-step fitting.
Now we want to optimize a more complex function that has a 3D input data space (x,y,z) and which is no longer linear to all of it's parameters: The function uses 6 parameters .
Even though the function is way more complex, the analytic derivative is still quite easy to estimate:
In case of this function using real parameters (1,2,3,4,5,6), and starting parameters (3,3,3,3,3,3), we see, that the numerical Jacbobian does only work with double precision. If we use float precision, the analytic Jacobian needs to be given in order to get convergence:
#include <ICLMath/LevenbergMarquardtFitter.h> using namespace icl::utils; using namespace icl::math; typedef float real; typedef icl::math::LevenbergMarquardtFitter<real> LM; LM::Vector f3(const LM::Params &p, const LM::Vector &vx){ real x = vx[0], y=vx[1], z=vx[2]; real a = p[0], b=p[1], c=p[2], d=p[3], e=p[4], f=p[5]; return LM::Vector(1,x*y*a*b + y*z*d*d*c - z*x*f*e*b + sqr(a + b + c + d)); } void j3(const LM::Params &p, const LM::Vector &vx, LM::Vector &dst){ real x = vx[0], y=vx[1], z=vx[2]; real a = p[0], b=p[1], c=p[2], d=p[3], e=p[4], f=p[5]; real k = 2 * (a+b+c+d); dst[0] = x*y*b + k; dst[1] = x*y*a - z*x*f*e + k; dst[2] = y*z*d*d + k; dst[3] = 2*d*y*z*c + k; dst[4] = -z*x*f*b; dst[5] = -z*x*e*b; } int main(){ LM lm(f3,1,std::vector<LM::Jacobian>(1,j3)); // use numerical jacobian LM lm(f); const real p[6] = { 1,2,3,4,5,6 }; LM::Params realParams(6,p), startParams(6,3); LM::Data data = LM::create_data(realParams,f3,3,1); LM::Result result = lm.fit(data.x,data.y,startParams); SHOW(result); }
The LevenbergMarquardtFitter class can also optimize functions with multi-dinensional output. In this case the internal optimization is performed for each output dimension seperately in each step. A common multidimensional optimization problem is 6D pose optimization. Please note, that there are many local minima (due to the complex relationship between params and output). Therefore, we need to start somewhere close to the original parameters. For the rotation parameters, it is always ususal to find different solutions that are higher or lower by n*pi or n*2pi. In this example, we did not try to derive an analytical jacobian so we use the numerical default jacobian here
#include <ICLMath/LevenbergMarquardtFitter.h> #include <ICLMath/FixedVector.h> using namespace icl::utils; using namespace icl::math; typedef float real; typedef icl::math::LevenbergMarquardtFitter<real> LM; LM::Vector f4(const LM::Params &p, const LM::Vector &vx){ FixedColVector<real,4> x(vx[0],vx[1], vx[2], 1); FixedMatrix<real,4,4> T = create_hom_4x4<float>(p[0],p[1],p[2],p[3],p[4],p[5]); FixedColVector<real,4> y = T*x; return LM::Vector(3,y.begin()); } int main(){ LM lm(f4,3); const real r[6] = { 3, 1, 2, 100, 100, 100}; const real s[6] = { 0.1, 0.1, 0.1, 10, 10, 10 }; LM::Params realParams(6,r), startParams(6,s); LM::Data data = LM::create_data(realParams,f4,3,3); LM::Result result = lm.fit(data.x,data.y,startParams); SHOW(result); }
typedef icl::utils::Function<void,const Result&> icl::math::LevenbergMarquardtFitter< Scalar >::DebugCallback |
Optionally given debug callback, that is called in every iterations.
typedef icl::utils::Function<Vector,const Params&,const Vector&> icl::math::LevenbergMarquardtFitter< Scalar >::Function |
to-be-optimized function type y = f(params, x)
typedef icl::utils::Function<void,const Params&, const Vector&, Vector &> icl::math::LevenbergMarquardtFitter< Scalar >::Jacobian |
jacobian of F
typedef DynMatrix<Scalar> icl::math::LevenbergMarquardtFitter< Scalar >::Matrix |
matrix type (used for input data)
typedef DynColVector<Scalar> icl::math::LevenbergMarquardtFitter< Scalar >::Params |
parameter vector type
typedef DynColVector<Scalar> icl::math::LevenbergMarquardtFitter< Scalar >::Vector |
vector type
icl::math::LevenbergMarquardtFitter< Scalar >::LevenbergMarquardtFitter | ( | ) |
creates a dummy (null instance)
icl::math::LevenbergMarquardtFitter< Scalar >::LevenbergMarquardtFitter | ( | Function | f, |
int | outputDim, | ||
const std::vector< Jacobian > & | js = std::vector< Jacobian >() , |
||
Scalar | initialLambda = 1.e-8 , |
||
int | maxIterations = 200 , |
||
Scalar | minError = 1.e-6 , |
||
Scalar | lambdaMultiplier = 10 , |
||
const std::string & | linSolver = "svd" |
||
) |
create an instance with given parameters
f | function to be optimized |
j | optionally given Jacobian of f. If j is null (e.g. an empty Jacobian() is passed), a numerical jacobian is created automatically. This will use a numerical delta of 1e-5, which is the default parameter for the static LevenbergMarquardtFitter::create_numerical_jacobian method. If a numerical Jacobian with another delta shall be used, create_numerical_jacobian can be called manually to obtain another automatically created numerical jacobian. |
initialLambda | initial damping parameter (usually 1e-6 ist not the worst choice) |
maxIterations | maximum number of iterations (usually, LevenbergMarquardtFitter will converge fast, or not at all. Therefore, the default argument of 200 iterations somehow assumes, that the minError criterion is met much earlier. |
minError | if the current error gets less than this threshold, the optimization is finished |
lambdaMultiplier | mulitiplyer used for increasing/decreasing the damping parameter lambda |
linSolver | linear solver that is used to estimate a local step internally possible values are documented in icl::DynMatrix::solve (we recommend the most-stable method svd) |
static Data icl::math::LevenbergMarquardtFitter< Scalar >::create_data | ( | const Params & | p, |
Function | f, | ||
int | xDim, | ||
int | yDim, | ||
int | num = 1000 , |
||
Scalar | minX = -5 , |
||
Scalar | maxX = 5 |
||
) | [static] |
creates test data using a given function
p | real function parameters |
f | function |
xDim | input dimension of function f |
num | number of samples |
minX | the function input values are randomly drawn from a uniform distribution in range [minX,maxX] |
maxX | see minX |
static Jacobian icl::math::LevenbergMarquardtFitter< Scalar >::create_numerical_jacobian | ( | int | o, |
Function | f, | ||
float | delta = 1.E-5 |
||
) | [static] |
creates a single numerical Jacobian for a given function f and output dim
The Jacobian will evaluate f twice for each parameter in . Here is the internal implementation snippet:
static std::vector<Jacobian> icl::math::LevenbergMarquardtFitter< Scalar >::create_numerical_jacobians | ( | int | n, |
Function | f, | ||
float | delta = 1.e-5 |
||
) | [static] |
creates a set of numerical jacobians for output dimension n
static void icl::math::LevenbergMarquardtFitter< Scalar >::default_debug_callback | ( | const Result & | r | ) | [static] |
default debug callback that simply streams r to std::cout
Scalar icl::math::LevenbergMarquardtFitter< Scalar >::error | ( | const Matrix & | ys, |
const Matrix & | y_est | ||
) | const [protected] |
returns the complete error
Result icl::math::LevenbergMarquardtFitter< Scalar >::fit | ( | const Matrix & | xs, |
const Matrix & | ys, | ||
Params | initParams | ||
) |
actual parameter fitting with given data and start parameters
This method actually fits the internal function model to the given data. The fitting procedure starts at the given set of initial parameters.
xs | input data matrix, where each data point is one row. Please note, that you can shallowly wrap a Matrix instance around existing other data types. |
ys | output data (Scalar) ys[i] belongs to the ith row of xs |
initParams | initial parameters for starting optimizaiton |
void icl::math::LevenbergMarquardtFitter< Scalar >::init | ( | Function | f, |
int | outputDim, | ||
const std::vector< Jacobian > & | js = std::vector< Jacobian >() , |
||
Scalar | initialLambda = 1.e-8 , |
||
int | maxIterations = 1000 , |
||
Scalar | minError = 1.e-6 , |
||
Scalar | lambdaMultiplier = 10 , |
||
const std::string & | linSolver = "svd" |
||
) |
(re)-initialization method
void icl::math::LevenbergMarquardtFitter< Scalar >::setDebugCallback | ( | DebugCallback | dbg = default_debug_callback | ) |
sets a debug callback method, which is called automatically in every interation
void icl::math::LevenbergMarquardtFitter< Scalar >::setUseMultiThreading | ( | bool | enable | ) |
enables openmp based multithreading
multithreading using 2 instead of 1 core provides a processing time reduction of about 35%. Please note, that openmp must be enabled using the compiler optionn "-fopenmp" as well
DebugCallback icl::math::LevenbergMarquardtFitter< Scalar >::dbg [protected] |
debug callback
Vector icl::math::LevenbergMarquardtFitter< Scalar >::dst [protected] |
b in Mx=b of linear system solved internally
Vector icl::math::LevenbergMarquardtFitter< Scalar >::dy [protected] |
buffer for current delta
Function icl::math::LevenbergMarquardtFitter< Scalar >::f [protected] |
Function f.
Matrix icl::math::LevenbergMarquardtFitter< Scalar >::H [protected] |
Hessian Matrix (J^T J)
Scalar icl::math::LevenbergMarquardtFitter< Scalar >::initialLambda [protected] |
initial damping parameter lambda
Matrix icl::math::LevenbergMarquardtFitter< Scalar >::J [protected] |
Jacobian Matrix (rows are Ji)
std::vector<Jacobian> icl::math::LevenbergMarquardtFitter< Scalar >::js [protected] |
output buffers
Scalar icl::math::LevenbergMarquardtFitter< Scalar >::lambdaMultiplier [protected] |
mulitplier that is used to increase/decreas lambda
std::string icl::math::LevenbergMarquardtFitter< Scalar >::linSolver [protected] |
linear solver that is used
int icl::math::LevenbergMarquardtFitter< Scalar >::maxIterations [protected] |
maximum number of iterations
Scalar icl::math::LevenbergMarquardtFitter< Scalar >::minError [protected] |
minimum error threshold
Params icl::math::LevenbergMarquardtFitter< Scalar >::params_new [protected] |
new parameters (after update step)
bool icl::math::LevenbergMarquardtFitter< Scalar >::useMultiThreading [protected] |
flag whether multithreading is enabled
Matrix icl::math::LevenbergMarquardtFitter< Scalar >::y_est [protected] |
current estimated outputs