D:/Programmierung/C/TRTK/include/TRTK/Optimization.hpp

Go to the documentation of this file.

/*
    Copyright (C) 2010 - 2014 Christoph Haenisch

    Chair of Medical Engineering (mediTEC)
    RWTH Aachen University
    Pauwelsstr. 20
    52074 Aachen
    Germany

    See license.txt for more information.

    Version 0.1.0 (2012-06-10)
*/

#ifndef OPTIMIZATION_HPP_4562354833
#define OPTIMIZATION_HPP_4562354833


#include <stdexcept>

#include <Eigen/Core>
#include <Eigen/Dense>

#include "Coordinate.hpp"
#include "Definitions.hpp"


namespace TRTK
{


namespace Optimization
{


template <class ValueType>
struct Options
{
    typedef ValueType value_type;

    value_type   error_tolerance;               
    value_type   lambda;                        
    unsigned     max_number_iterations;         
    value_type   nu;                            
    value_type   spacing;                       

    enum Algorithm
    {
        LEVENBERG_MARQUARDT,                    
        NEWTON_RAPHSON                          
    } algorithm;                                

    Options() :
        error_tolerance(1e-15),
        lambda(1e-2),
        max_number_iterations(25),
        nu(10),
        spacing(1e-6),
        algorithm(LEVENBERG_MARQUARDT)
    {
    }
};


template <class ValueType>
struct Result
{
    typedef Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic> matrix_type;
    typedef TRTK::Coordinate<ValueType> Coordinate;
    typedef ValueType value_type;

    value_type   error;
    matrix_type  jacobian;
    unsigned     number_of_iterations;
    Coordinate   root;

    Result() : error(0), number_of_iterations(0)
    {
    }
};


template <class Function, class ValueType>
class Gradient
{
public:
    Gradient(Function & function, Options<ValueType> options = Options<ValueType>());
    virtual ~Gradient();

    Coordinate<ValueType> operator()(const Coordinate<ValueType> & value);

private:
    Function & function;
    Options<ValueType> options;
};


template<class Function, class ValueType>
Gradient<Function, ValueType>::Gradient(Function & function, Options<ValueType> options) :
    function(function),
    options(options)
{
}


template<class Function, class ValueType>
Gradient<Function, ValueType>::~Gradient()
{
}


template<class Function, class ValueType>
Coordinate<ValueType> Gradient<Function, ValueType>::operator()(const Coordinate<ValueType> & value)
{
    // function : R^n --> R

    const unsigned n = value.size();

    Coordinate<ValueType> gradient;
    gradient.resize(n);

    #ifdef TRTK_PARALLELIZE_GRADIENT
    #pragma omp parallel for
    #endif // TRTK_PARALLELIZE_GRADIENT

    for (int i = 0; i < int(n); ++i)
    {
        // Use the five-point stencil method to compute the derivatives.

        const ValueType spacing = options.spacing;

        Coordinate<ValueType> h;
        h.resize(n);
        h[i] = spacing; // = (0, ..., 0, spacing, 0, ..., 0)

        ValueType a = function(value - 2.0 * h);
        ValueType b = function(value -       h);
        ValueType c = function(value +       h);
        ValueType d = function(value + 2.0 * h);

        gradient[i] = (a - 8.0 * b + 8.0 * c - d) / (12.0 * spacing);
    }

    return gradient;
}


template <class Function, class ValueType>
Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic> jacobian(Function f, const Coordinate<ValueType> & x, Options<ValueType> options = Options<ValueType>())
{
    using namespace Eigen;

    typedef TRTK::Coordinate<ValueType> Coordinate;

    // function : R^n --> R^m

    const unsigned n = x.size();
    const unsigned m = f(x).size();

    Matrix<ValueType, Dynamic, Dynamic> jacobian(m, n);

    #ifdef TRTK_PARALLELIZE_JACOBIAN
    #pragma omp parallel for
    #endif // TRTK_PARALLELIZE_JACOBIAN

    for (int col = 0; col < int(n); ++col)
    {
        // Use the five-point stencil method to compute the derivatives.

        const ValueType spacing = options.spacing;

        Coordinate h;
        h.resize(n);
        h[col] = spacing; // = (0, ..., 0, spacing, 0, ..., 0)

        Coordinate a = f(x - 2.0 * h);
        Coordinate b = f(x -       h);
        Coordinate c = f(x +       h);
        Coordinate d = f(x + 2.0 * h);

        Coordinate derivatives_in_col = (a - 8.0 * b + 8.0 * c - d) / (12.0 * spacing);

        for (unsigned row = 0; row < m; ++row)
        {
            jacobian(row, col) = derivatives_in_col[row];
        }
    }

    return jacobian;
}


template <class ValueType, class Function>
Gradient<Function, ValueType> makeGradient(Function & function, Options<ValueType> options = Options<ValueType>())
{
    return Gradient<Function, ValueType>(function, options);
}


template <class Function, class ValueType>
Result<ValueType> solve(Function function, const Coordinate<ValueType> & start_value, Options<ValueType> options = Options<ValueType>())
{
    switch(options.algorithm)
    {
        case Options<ValueType>::LEVENBERG_MARQUARDT:
            return solveLevenbergMarquardt(function, start_value, options);

        case Options<ValueType>::NEWTON_RAPHSON:
            return solveLevenbergMarquardt(function, start_value, options);

        default:
            throw std::invalid_argument("TRTK::Optimization::solve(): Unknown solver algorithm.");
    }
}


template<class UnaryFunction, class ValueType>
ValueType solve1D(UnaryFunction f, ValueType start_value, unsigned max_number_iterations = 25)
{
    const ValueType epsilon = 0.000001;

    ValueType x = start_value;

    for (unsigned i = 0; i < max_number_iterations; ++i)
    {
        ValueType a = f(x - 2 * epsilon);
        ValueType b = f(x -     epsilon);
        ValueType c = f(x +     epsilon);
        ValueType d = f(x + 2 * epsilon);

        ValueType f_prime = (a - 8 * b + 8 * c - d) / (12 * epsilon); // five-point stencil method

        if (f_prime == 0) break; // Avoid a division by zero.

        x = x - f(x) / f_prime;
    }

    return x;
}


template <class Function, class ValueType>
Result<ValueType> solveLevenbergMarquardt(Function function, const Coordinate<ValueType> & start_value, Options<ValueType> options = Options<ValueType>())
{
    // The algorithm is implemented as described in the original article
    // of Marquardt. However, instead of scaling the matrix A and the vector
    // g we add lambda * diag(A) to A rather than lambda * I since this is
    // equivalent but more efficient.

    using std::numeric_limits;
    using std::sqrt;

    using namespace Eigen;

    typedef TRTK::Coordinate<ValueType> Coordinate;
    typedef Eigen::Matrix<ValueType, Dynamic, 1> Vector;
    typedef Eigen::Matrix<ValueType, Dynamic, Dynamic> Matrix;

    unsigned iteration = 0;
    const unsigned m = function(start_value).size();    // function : R^n --> R^m
    ValueType error = 0;
    ValueType lambda = options.lambda;                  // damping parameter
    ValueType nu = options.nu;                          // damping parameter
    ValueType Phi = numeric_limits<ValueType>::max();
    Coordinate x = start_value;
    Matrix jacobian;

    assert(lambda > 0);
    assert(nu > 1);

    for (iteration = 0; iteration < options.max_number_iterations; ++iteration)
    {
        // Compute Phi(lambda) and Phi(lambda/nu).

        // Now, compute (A + lambda * diag(A)) * delta = g where A = Jacobian(f)^T *
        // Jacobian(f) and g = Jacobian(f)^T * f. Split the computation due to
        // stability reasons; this also leads to a speed up since no inverse must
        // be computed.

        Vector x1 = x.toArray(); // --> Phi(lambda / nu)
        Vector x2 = x.toArray(); // --> Phi(lambda)
        Vector f  = function(x).toArray();

        jacobian = Optimization::jacobian(function, x, options);

        if (jacobian.norm() == 0)
        {
            error = function(x).norm() / sqrt(ValueType(m)); // RMSE
            break;
        }

        Matrix jj = jacobian.transpose() * jacobian;
        Matrix diagonal = Matrix(jj.diagonal().asDiagonal());
        Matrix left_term1 = jj + lambda / nu * diagonal;
        Matrix left_term2 = jj + lambda * diagonal;
        Vector right_term = jacobian.transpose() * -f;

        // Same as "x1 = x1 - (jj + lambda / nu * diagonal).inverse() * jacobian.transpose() * f;"
        // but faster and more stable:

        Vector delta = left_term1.fullPivLu().solve(right_term);
        x1 = delta + x1;

        delta = left_term2.fullPivLu().solve(right_term);
        x2 = delta + x2;

        ValueType Phi1 = function(Coordinate(x1)).squaredNorm();
        ValueType Phi2 = function(Coordinate(x2)).squaredNorm();

        // Case-by-case analysis

        if (Phi1 <= Phi)
        {
            // case i. in the original article

            lambda = lambda / nu;
            Phi = Phi1;
            x = Coordinate(x1);
        }
        else
        {
            if (Phi2 <= Phi)
            {
                // case ii. in the original article

                lambda = lambda;
                Phi = Phi2;
                x = Coordinate(x2);
            }
            else
            {
                // case iii. in the original article

                Vector x3 = x.toArray();

                ValueType lambda3 = lambda;
                ValueType Phi3 = Phi;

                unsigned w = 0;
                const unsigned MAX_POWER = 100; // this bound is arbitrarily chosen...

                do
                {
                    lambda3 *= nu; // increase lambda by successive multiplication

                    Matrix left_term = jj + lambda3 * diagonal;
                    x3 += left_term.fullPivLu().solve(right_term);

                    Phi3 = function(Coordinate(x3)).squaredNorm();
                }
                while (!(Phi3 <= Phi || w++ >= MAX_POWER));

                if (w > MAX_POWER)
                {
                    // increasing lambda does not seem to work, hence do nothing
                    // and keep the old value (Phi3 might now be even worse...)
                }
                else
                {
                    lambda = lambda3;
                    Phi = Phi3;
                    x = Coordinate(x3);
                }
            }
        }

        // Converged or small enough RMSE?

        Coordinate return_value = function(x);

        // Check whether the error (RMSE) is less than the given error
        // tolerance. Should the situation arise, stop the iteration.

        error = return_value.norm() / sqrt(ValueType(m));

        if (error < options.error_tolerance) break;
    }

    Result<ValueType> result;

    result.error = error;
    result.jacobian = jacobian;
    result.number_of_iterations = iteration + 1;
    result.root = x;

    return result;
}


template <class Function, class ValueType>
Result<ValueType> solveNewtonRaphsonMethod(Function function, const Coordinate<ValueType> & start_value, Options<ValueType> options = Options<ValueType>())
{
    using namespace Eigen;

    using std::sqrt;

    typedef TRTK::Coordinate<ValueType> Coordinate;

    // function : R^n --> R^m

    // const unsigned n = start_value.size();
    const unsigned m = function(start_value).size();

    Coordinate x = start_value;

    Matrix<ValueType, Dynamic, Dynamic> jacobian;

    ValueType error = 0;

    unsigned iteration = 0;

    for (iteration = 0; iteration < options.max_number_iterations; ++iteration)
    {
        Coordinate return_value = function(x);

        // Check whether the error (RMS) is less than the given error
        // tolerance. Should the situation arise, stop the iteration.

        error = return_value.norm() / sqrt(ValueType(m));

        if (error < options.error_tolerance) break;

        jacobian = Optimization::jacobian(function, x, options);

        if (jacobian.norm() == 0) break;

        Matrix<ValueType, Dynamic, 1> x_ = x.toArray();
        Matrix<ValueType, Dynamic, 1> f_ = function(x).toArray();

        // Same as "x_ = x_ - Jacobian.inverse() * f_;" but faster and more stable:

        Matrix<ValueType, Dynamic, 1> delta = jacobian.fullPivLu().solve(-f_);
        x_ = delta + x_;

        x = Coordinate(x_);
    }

    Result<ValueType> result;

    result.error = error;
    result.jacobian = jacobian;
    result.number_of_iterations = iteration;
    result.root = x;

    return result;
}


} // namespace Optimization


} // namespace TRTK


#endif // OPTIMIZATION_HPP_4562354833