This namespace contains various classes and functions to find the roots (zeros) or the local minima and maxima of a function. More...

Classes
struct	Options
	Structure used by several functions to control their operation. More...
struct	Result
	Structure returned by several functions. More...
class	Gradient
	Generates a functor which is the gradient of the given function. More...
Functions
template<class Function , class ValueType >
Eigen::Matrix< ValueType, Eigen::Dynamic, Eigen::Dynamic >	jacobian (Function f, const Coordinate< ValueType > &x, Options< ValueType > options=Options< ValueType >())
	Computes the Jacobian of f at x.
template<class ValueType , class Function >
Gradient< Function, ValueType >	makeGradient (Function &function, Options< ValueType > options=Options< ValueType >())
	Returns an instance of Gradient. The template argument types are automatically deduced.
template<class Function , class ValueType >
Result< ValueType >	solve (Function function, const Coordinate< ValueType > &start_value, Options< ValueType > options=Options< ValueType >())
	This function solves for the root (zero) of a given function.
template<class UnaryFunction , class ValueType >
ValueType	solve1D (UnaryFunction f, ValueType start_value, unsigned max_number_iterations=25)
	This function solves for the root (zero) of a given function.
template<class Function , class ValueType >
Result< ValueType >	solveLevenbergMarquardt (Function function, const Coordinate< ValueType > &start_value, Options< ValueType > options=Options< ValueType >())
	This function solves for the root (zero) of a given function.
template<class Function , class ValueType >
Result< ValueType >	solveNewtonRaphsonMethod (Function function, const Coordinate< ValueType > &start_value, Options< ValueType > options=Options< ValueType >())
	This function solves for the root (zero) of a given function.

Detailed Description

This namespace contains various classes and functions to find the roots (zeros) or the local minima and maxima of a function.

List of classes and functions:

Macros:

If TRTK_PARALLELIZE is defined some algorithms will make use of OpenMP.

Example:

The following example shows how these functions can be used to estimate an unknown transformation between two point sets:

 #include <cassert>
 #include <cmath>
 #include <iomanip>
 #include <iostream>
 #include <vector>

 #include <TRTK/Clock.hpp>
 #include <TRTK/Coordinate.hpp>
 #include <TRTK/Optimization.hpp>
 #include <TRTK/Tools.hpp>
 #include <TRTK/Transform2D.hpp>

 using namespace std;
 using namespace TRTK;
 using namespace TRTK::Optimization;
 using namespace TRTK::Tools;


 class EstimateTransformation
 {
 public:
     EstimateTransformation(const vector<Coordinate<double> > & source_points,
                            const vector<Coordinate<double> > & target_points);

     // Returns the error of the transformation which is defined by the given coefficients.
     double operator()(const Coordinate<double> & coefficients) const;

 private:
     const vector<Coordinate<double> > & source_points;
     const vector<Coordinate<double> > & target_points;
 };


 EstimateTransformation::EstimateTransformation(const vector<Coordinate<double> > & source_points,
                                                const vector<Coordinate<double> > & target_points) :
     source_points(source_points),
     target_points(target_points)
 {
 }


 double EstimateTransformation::operator()(const Coordinate<double> & coefficients) const
 {
     // Returns the Mean Square Error (MSE) of the currently estimated transformation.

     Transform2D<double> transform(coefficients[0], coefficients[1], coefficients[2],
                                   coefficients[3], coefficients[4], coefficients[5],
                                   0,               0,               1);

     double error = 0;

     const unsigned N = source_points.size();

     for (unsigned i = 0; i < N; ++i)
     {
         double residual = (transform * source_points[i] - target_points[i]).norm();
         error += residual * residual;
     }

     return error / N;
 }


 int main()
 {
     // Generate an abitrary transformation.

     Transform2D<double> transform;

     transform.rotate(0.927295218).shear(1, 0).scale(2, 2).translate(1, -2);

     // Construct some source points and some target points which
     // are the transformed source points.

     const unsigned N = 200;

     vector<Coordinate<double> > source_points;
     vector<Coordinate<double> > target_points;

     for (unsigned i = 0; i < N; ++i)
     {
         double x = rand(-100.0, 100.0);
         double y = rand(-100.0, 100.0);

         Coordinate<double> source_point(x, y);
         Coordinate<double> target_point = transform * source_point;

         source_points.push_back(source_point);
         target_points.push_back(target_point);
     }

     // Estimate the transformation.

     Clock clock;

     EstimateTransformation estimateTransformation(source_points, target_points);

     Coordinate<double> coefficients; // = {a11, a12, a13, a21, a22, a23} in Transform2D
     coefficients.resize(6, 1);

     Options<double> options;
     options.error_tolerance = 1e-10;

     Result<double> result = solve(makeGradient<double>(estimateTransformation),
                                   coefficients,
                                   options);

     // Print the result.

     cout << "Original transformation matrix:" << endl << endl;

     cout << setw(6) << transform.a11() << setw(6) << transform.a12() << setw(6) << transform.a13() << endl
          << setw(6) << transform.a21() << setw(6) << transform.a22() << setw(6) << transform.a23() << endl
          << setw(6) << transform.a31() << setw(6) << transform.a32() << setw(6) << transform.a33() << endl
          << endl;

     ios_base::fmtflags flags = cout.flags();

     cout << "Result:    " << endl << endl
          << "Error      " << result.error << endl
          << "Iterations " << result.number_of_iterations << endl
          << "Root       " << fixed << setprecision(3) << result.root << endl << endl;

     cout.flags(flags);

     cout << clock << endl;

     return 0;
 }

The output is:

 Original transformation matrix:

    2.8  -0.4     1
    1.6   1.2    -2
      0     0     1

 Result:

 Error      8.0176e-14
 Iterations 3
 Root       (2.800, -0.400, 1.000, 1.600, 1.200, -2.000)

 Elapsed time: 0.11 seconds.

Function Documentation

template<class Function , class ValueType >

Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic> TRTK::Optimization::jacobian	(	Function	f,
		const Coordinate< ValueType > &	x,
		Options< ValueType >	options = `Options<ValueType>()`
	)

Computes the Jacobian of f at x.

Parameters:

[in]	f	Vector-valued multivariate function or functor.
[in]	x	Vector-valued input argument.
[in]	options	Options to control the operation of this function.

The Jacobian of function is computed numerically using a higher order method, namely the five-point stencil method. The spacing used when computing the finit difference can be set in options.

Note:: The type of the input argument as well as the type of the return value of the given function must be Coordinate.

Example:

 #include <cassert>
 #include <cmath>
 #include <iostream>

 #include <TRTK/Coordinate.hpp>
 #include <TRTK/Optimization.hpp>

 using namespace std;
 using namespace TRTK;
 using namespace TRTK::Optimization;

 Coordinate<double> f(const Coordinate<double> & arg)
 {
     assert(arg.size() == 2);

     const double x = arg.x();
     const double y = arg.y();

     Coordinate<double> result(0, 0);

     result.x() = (x - 2) * (y - 1);
     result.y() = -3 * y;

     return result;
 }

 int main()
 {
     Coordinate<double> arg(5, 5);
     cout << jacobian(f, arg);
     return 0;
 }

Output:

 4  3
 0 -3

Macros:: If TRTK_PARALLELIZE_JACOBIAN is defined jacobian will make use of OpenMP. Be aware that f needs to be reentrant (i.e. it can be safely executed concurrently).

Author:: Christoph Haenisch

Version:: 0.1.1

Date:: last changed on 2013-08-15

Definition at line 503 of file Optimization.hpp.

template<class ValueType , class Function >

Gradient<Function, ValueType> TRTK::Optimization::makeGradient	(	Function &	function,
		Options< ValueType >	options = `Options<ValueType>()`
	)

Returns an instance of Gradient. The template argument types are automatically deduced.

Please, have a look at Gradient for more details.

Example:

 #include <cassert>
 #include <cmath>
 #include <iostream>

 #include <TRTK/Coordinate.hpp>
 #include <TRTK/Optimization.hpp>

 using namespace std;
 using namespace TRTK;
 using namespace TRTK::Optimization;

 double f(const TRTK::Coordinate<double> & arg)
 {
     assert(arg.size() == 2);

     const double x = arg.x();
     const double y = arg.y();

     return (x - 2) * (x - 2)  - (y - 1) * (y - 1) - 9;
 }

 int main()
 {
     // Find the minimum of f (i.e. the zero of grad f).

     Coordinate<double> arg(100, -2); // initial guess

     Result<double> result = solveNewtonRaphsonMethod(makeGradient<double>(f), arg);

     cout << "Result:    " << endl
          << "Error      " << result.error << endl
          << "Iterations " << result.number_of_iterations << endl
          << "Root       " << result.root << endl;

     return 0;
 }

Output:

 Result:
 Error      0
 Iterations 3
 Root       (2, 1)

Example:

 class A
 {
 public:
     double f(const Coordinate<double> & arg) const
     {
         // ...
     }

     double operator()(const Coordinate<double> & arg) const
     {
         return f(arg);
     }
 };

 int main()
 {
     // ...
     Result<double> result = solveNewtonRaphsonMethod(makeGradient<double>(a), Coordinate<double> arg(0, 0));
     // ...
 }

Author:: Christoph Haenisch

Version:: 0.1.1

Date:: last changed on 2012-11-07

Definition at line 632 of file Optimization.hpp.

template<class Function , class ValueType >

Result<ValueType> TRTK::Optimization::solve	(	Function	function,
		const Coordinate< ValueType > &	start_value,
		Options< ValueType >	options = `Options<ValueType>()`
	)

This function solves for the root (zero) of a given function.

Parameters:

[in]	function	Vector-valued multivariate function or functor.
[in]	start_value	Value in the neighborhood of a root.
[in]	options	Options to control the operation of this function.

This functions serves as a uniform interface for all solver functions.

By default the solver uses the Levenberg-Marquardt algorithm. This can be changed with the options parameter. Also, all parameters relevant for a particular algorithm can be set with an option object.

For further details please have a look at the respective function documentation.

Example:

 #include <cassert>
 #include <iostream>

 #include <TRTK/Coordinate.hpp>
 #include <TRTK/Optimization.hpp>

 using namespace std;
 using namespace TRTK;
 using namespace TRTK::Optimization;

 Coordinate<double> f(const Coordinate<double> & arg)
 {
     assert(arg.size() == 2);

     const double x = arg.x();
     const double y = arg.y();

     Coordinate<double> result(0, 0);

     result.x() = (x - 2) * (y - 1);
     result.y() = -3 * y;

     return result;
 }

 int main()
 {
     Coordinate<double> arg(100, -2);

     Result<double> result = solve(f, arg);

     cout << "Result:    " << endl
          << "Iterations " << result.number_of_iterations << endl
          << "Error      " << result.error << endl
          << "Root       " << result.root << endl;

     return 0;
 }

Output:

 Result: 
 Iterations 5
 Error      2.2523e-023
 Root       (2, 1.06174e-023)

See also:: solveLevenbergMarquardt(), solveNewtonRaphsonMethod()

Author:: Christoph Haenisch

Version:: 0.1.0

Date:: last changed on 2012-11-09

Definition at line 712 of file Optimization.hpp.

template<class UnaryFunction , class ValueType >

ValueType TRTK::Optimization::solve1D	(	UnaryFunction	f,
		ValueType	start_value,
		unsigned	max_number_iterations = `25`
	)

This function solves for the root (zero) of a given function.

Parameters:

[in]	f	Unary function or functor (e.g. `double f(double)`).
[in]	start_value	Value in the neighborhood of a root.
[in]	max_number_iterations	Maximum number of iterations to perform.

This function implements the Newton-Raphson method. In order to converge, the initial estimate (start_value) must be sufficient close to the root. At most max_number_iterations iterations are performed. The derivative of f is computed numerically using a higher order method (five-point stencil).

Example:

 #include <TRTK/Coordinate.hpp>
 #include <TRTK/Optimization.hpp>

 using namespace TRTK;
 using namespace TRTK::Optimization;

 double f(double x)
 {
     return x * x - 2;
 }

 int main()
 {
     solve1D(f, 1.0); // returns sqrt(2)
     return 0;
 }

Author:: Christoph Haenisch

Version:: 0.1.0

Date:: last changed on 2012-11-07

Definition at line 767 of file Optimization.hpp.

template<class Function , class ValueType >

Result<ValueType> TRTK::Optimization::solveLevenbergMarquardt	(	Function	function,
		const Coordinate< ValueType > &	start_value,
		Options< ValueType >	options = `Options<ValueType>()`
	)

This function solves for the root (zero) of a given function.

Parameters:

[in]	function	Vector-valued multivariate function or functor.
[in]	start_value	Value in the neighborhood of a root.
[in]	options	Options to control the operation of this function.

This function implements the Levenberg-Marquardt algorithm [1] which aims to minimize the error criterion

\[ \Phi = \| f(x) \|^2 = \sum_{i=0}^{n} f_i(x)^2 \]

by iteratively solving the equation

\[ (JJ^T + \lambda \text{diag}(JJ^T)) \Delta x = -J^T f(x) \]

where \( f(x + \Delta x) \approx f(x) + J \Delta x \).

The Levenberg-Marquardt algorithm can be interpreted as interpolating between the Gradient Descent and the Newton-Raphson method depending on the damping parameter \( \lambda \). A high value of \( \lambda \) leads to a more Gradient Descent like characteristic which is advantageous in the case of a suboptimal start value. The parameter is adjusted at each iteration. The damping can be controlled with the paramters options.lambda and options.nu.

In order to converge, the initial estimate (start_value) must be sufficient close to the root. The algorithm is terminated if the residual is smaller than the given error tolerance or if the maximum number of iterations is reached.

The Jacobian of function is computed numerically using a higher order method (i.e. five-point stencil method). The spacing used when computing the finit difference can be set in options.

Note:: The type of the input argument as well as the type of the return value of the given function must be Coordinate.

Example:

 #include <cassert>
 #include <iostream>

 #include <TRTK/Coordinate.hpp>
 #include <TRTK/Optimization.hpp>

 using namespace std;
 using namespace TRTK;
 using namespace TRTK::Optimization;

 Coordinate<double> f(const Coordinate<double> & arg)
 {
     assert(arg.size() == 2);

     const double x = arg.x();
     const double y = arg.y();

     Coordinate<double> result(0, 0);

     result.x() = (x - 2) * (y - 1);
     result.y() = -3 * y;

     return result;
 }

 int main()
 {
     Coordinate<double> arg(100, -2);

     Result<double> result = solveLevenbergMarquardt(f, arg);

     cout << "Result:    " << endl
          << "Iterations " << result.number_of_iterations << endl
          << "Error      " << result.error << endl
          << "Root       " << result.root << endl;

     return 0;
 }

Output:

 Result: 
 Iterations 5
 Error      2.2523e-023
 Root       (2, 1.06174e-023)

References:

[1] Donald W. Marquardt, "An Algorithm for Least Squares Estimation of Nonlinear Parameters", Journal of the Society dor Industrial and Applied Mathematics, Vol. 11, No. 2, 1963

Author:: Christoph Haenisch

Version:: 0.1.0

Date:: last changed on 2012-11-07

Bug:: There is a known bug in the implementation. Do not use this algorithm until the bug is fixed.

Definition at line 894 of file Optimization.hpp.

template<class Function , class ValueType >

Result<ValueType> TRTK::Optimization::solveNewtonRaphsonMethod	(	Function	function,
		const Coordinate< ValueType > &	start_value,
		Options< ValueType >	options = `Options<ValueType>()`
	)

This function solves for the root (zero) of a given function.

Parameters:

[in]	function	Vector-valued multivariate function or functor.
[in]	start_value	Value in the neighborhood of a root.
[in]	options	Options to control the operation of this function.

This function implements the Newton-Raphson method. In order to converge, the initial estimate (start_value) must be sufficient close to the root. The algorithm is terminated if the residual is smaller than the given error tolerance or if the maximum number of iterations is reached.

The Jacobian of function is computed numerically using a higher order method (i.e. five-point stencil method). The spacing used when computing the finit difference can be set in options.

Note:: The type of the input argument as well as the type of the return value of the given function must be Coordinate.

Example:

 #include <cassert>
 #include <iostream>

 #include <TRTK/Coordinate.hpp>
 #include <TRTK/Optimization.hpp>

 using namespace std;
 using namespace TRTK;
 using namespace TRTK::Optimization;

 Coordinate<double> f(const Coordinate<double> & arg)
 {
     assert(arg.size() == 2);

     const double x = arg.x();
     const double y = arg.y();

     Coordinate<double> result(0, 0);

     result.x() = (x - 2) * (y - 1);
     result.y() = -3 * y;

     return result;
 }

 int main()
 {
     Coordinate<double> arg(100, -2);

     Result<double> result = solveNewtonRaphsonMethod(f, arg);

     cout << "Result:    " << endl
          << "Iterations " << result.number_of_iterations << endl
          << "Error      " << result.error << endl
          << "Root       " << result.root << endl;

     return 0;
 }

Output:

 Result:
 Iterations 3
 Error      1.04589e-31
 Root       (2, -4.93038e-32)

Author:: Christoph Haenisch

Version:: 0.1.0

Date:: last changed on 2012-06-10

Definition at line 1118 of file Optimization.hpp.

TRTK::Optimization Namespace Reference

Classes

Functions

Detailed Description

Function Documentation