Image Component Library (ICL)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Classes | Public Types | Public Member Functions | Static Public Member Functions | Protected Member Functions | Protected Attributes
icl::math::LevenbergMarquardtFitter< Scalar > Class Template Reference

Utility class implementing the multidimensional Levenberg Marquardt Algorithm for non-linear Optimization. More...

#include <LevenbergMarquardtFitter.h>

List of all members.

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< Jacobiancreate_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< Jacobianjs
 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

Detailed Description

template<class Scalar>
class icl::math::LevenbergMarquardtFitter< Scalar >

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 $ f : R^N \rightarrow R^O $, an intial set of parameters $ \beta_0 $ and a set of training data $ \{(x_i,y_i)\}$, where $ f $ depends on model parameters $ \beta $, LMA will find a local minumum of function parameters that minimizes

\[ S(\beta)​ = \sum_{o=1}^O \sum_{i=1}^m [y_i[o] - f(x_i, \ \beta)[o]​ ]​^2 \]

The complete algorithm can be found at at wikipedia (see http://en.wikipedia.org/wiki/Levenberg-Marquardt_algorithm)

What type of Function can be optimized

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)

\[ f(x) = \beta_1 + \beta_2 x + \beta_3 x^2 + beta_4 x^3 \]

Then LMA will find find the solution in a single step (depened on the step parameters lambda)

What is the Jacobian

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 $ f $ is defined by a vector (here, a row-vector) of partial derivations where the ith component is $ \partial f_o / \partial \beta_i$. 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 $ x_i $ and for each output dimension. Therefore, the Jacobian J_o (for output dimension o) becomes a matrix, where each row $ J_{oi] $ contains the partial derivations for all parameters beta for a sigle data point $ x_i $ \

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 $ f $ twice, and is usually not as exact as the analytic Jacobian, which howver can lead to faster or slower convergence.

Examples

Here, a short example is given, that shows how to use the LevenbergMarquardtFitter class.

Example 1

The input data space and the outputspace are 1D, and our function has 4 parameters (a,b,c,d) :

\[ f(x) = a + b*x + c*x^2 + d*x^3 \]

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

\[ ( 1, x, x^2, x^3 ) \]

        #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.

Second Example

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 $ \beta = (a,b,c,d,e,f) $.

\[ f(x,y,z;\beta) = xyab + yzd^2c - zxfeb + (a + b + c + d)^2; \]

Even though the function is way more complex, the analytic derivative is still quite easy to estimate:

\[ \frac{\partial f}{\partial a} = xyb + 2 (a + b + c + d) \]

\[ \frac{\partial f}{\partial b} = xya - zxfe + 2 (a + b + c + d); \]

\[ \frac{\partial f}{\partial c} = yzd^2 + 2 (a + b + c + d)\]

\[ \frac{\partial f}{\partial d} = 2dyzc + 2 (a + b + c + d)\]

\[ \frac{\partial f}{\partial e} = -zxfb \]

\[ \frac{\partial f}{\partial f} = zxeb \]

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

Multidimensional Output

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


Member Typedef Documentation

template<class Scalar >
typedef icl::utils::Function<void,const Result&> icl::math::LevenbergMarquardtFitter< Scalar >::DebugCallback

Optionally given debug callback, that is called in every iterations.

template<class Scalar >
typedef icl::utils::Function<Vector,const Params&,const Vector&> icl::math::LevenbergMarquardtFitter< Scalar >::Function

to-be-optimized function type y = f(params, x)

template<class Scalar >
typedef icl::utils::Function<void,const Params&, const Vector&, Vector &> icl::math::LevenbergMarquardtFitter< Scalar >::Jacobian

jacobian of F

See also:
What is the Jacobian
template<class Scalar >
typedef DynMatrix<Scalar> icl::math::LevenbergMarquardtFitter< Scalar >::Matrix

matrix type (used for input data)

template<class Scalar >
typedef DynColVector<Scalar> icl::math::LevenbergMarquardtFitter< Scalar >::Params

parameter vector type

template<class Scalar >
typedef DynColVector<Scalar> icl::math::LevenbergMarquardtFitter< Scalar >::Vector

vector type


Constructor & Destructor Documentation

template<class Scalar >
icl::math::LevenbergMarquardtFitter< Scalar >::LevenbergMarquardtFitter ( )

creates a dummy (null instance)

template<class Scalar >
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

Parameters:
ffunction to be optimized
joptionally 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.
initialLambdainitial damping parameter (usually 1e-6 ist not the worst choice)
maxIterationsmaximum 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.
minErrorif the current error gets less than this threshold, the optimization is finished
lambdaMultipliermulitiplyer used for increasing/decreasing the damping parameter lambda
linSolverlinear 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)

Member Function Documentation

template<class Scalar >
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

Parameters:
preal function parameters
ffunction
xDiminput dimension of function f
numnumber of samples
minXthe function input values are randomly drawn from a uniform distribution in range [minX,maxX]
maxXsee minX
template<class Scalar >
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 $ \beta $. Here is the internal implementation snippet:

          void jacobian(const Params &params, const Vector &x,Vector &target) const{
            Vector p = params;
            for(unsigned int i=0;i<params.dim();++i){
              p[i] = params[i] + delta/2;
              Scalar f1 = f(p,x)[o];
              p[i] = params[i] - delta/2;
              Scalar f2 = f(p,x)[o];
              p[i] = params[i];
              target[i] = ( f1 - f2 ) / delta;
            }
          }
template<class Scalar >
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

template<class Scalar >
static void icl::math::LevenbergMarquardtFitter< Scalar >::default_debug_callback ( const Result r) [static]

default debug callback that simply streams r to std::cout

template<class Scalar >
Scalar icl::math::LevenbergMarquardtFitter< Scalar >::error ( const Matrix ys,
const Matrix y_est 
) const [protected]

returns the complete error

template<class Scalar >
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.

Parameters:
xsinput data matrix, where each data point is one row. Please note, that you can shallowly wrap a Matrix instance around existing other data types.
ysoutput data (Scalar) ys[i] belongs to the ith row of xs
initParamsinitial parameters for starting optimizaiton
template<class Scalar >
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

sets a debug callback method, which is called automatically in every interation

template<class Scalar >
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


Member Data Documentation

template<class Scalar >
DebugCallback icl::math::LevenbergMarquardtFitter< Scalar >::dbg [protected]

debug callback

template<class Scalar >
Vector icl::math::LevenbergMarquardtFitter< Scalar >::dst [protected]

b in Mx=b of linear system solved internally

template<class Scalar >
Vector icl::math::LevenbergMarquardtFitter< Scalar >::dy [protected]

buffer for current delta

template<class Scalar >
Function icl::math::LevenbergMarquardtFitter< Scalar >::f [protected]

Function f.

template<class Scalar >
Matrix icl::math::LevenbergMarquardtFitter< Scalar >::H [protected]

Hessian Matrix (J^T J)

template<class Scalar >
Scalar icl::math::LevenbergMarquardtFitter< Scalar >::initialLambda [protected]

initial damping parameter lambda

template<class Scalar >
Matrix icl::math::LevenbergMarquardtFitter< Scalar >::J [protected]

Jacobian Matrix (rows are Ji)

template<class Scalar >
std::vector<Jacobian> icl::math::LevenbergMarquardtFitter< Scalar >::js [protected]

output buffers

template<class Scalar >
Scalar icl::math::LevenbergMarquardtFitter< Scalar >::lambdaMultiplier [protected]

mulitplier that is used to increase/decreas lambda

template<class Scalar >
std::string icl::math::LevenbergMarquardtFitter< Scalar >::linSolver [protected]

linear solver that is used

See also:
icl::DynMatrix::solve
template<class Scalar >
int icl::math::LevenbergMarquardtFitter< Scalar >::maxIterations [protected]

maximum number of iterations

template<class Scalar >
Scalar icl::math::LevenbergMarquardtFitter< Scalar >::minError [protected]

minimum error threshold

template<class Scalar >
Params icl::math::LevenbergMarquardtFitter< Scalar >::params_new [protected]

new parameters (after update step)

template<class Scalar >
bool icl::math::LevenbergMarquardtFitter< Scalar >::useMultiThreading [protected]

flag whether multithreading is enabled

template<class Scalar >
Matrix icl::math::LevenbergMarquardtFitter< Scalar >::y_est [protected]

current estimated outputs


The documentation for this class was generated from the following file:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines