/home/hauberg/Dokumenter/Capture/humim-tracker-0.1/src/OpenTissue/OpenTissue/core/math/optimization/optimization_bfgs.h

Go to the documentation of this file.

#ifndef OPENTISSUE_CORE_MATH_OPTIMIZATION_BFGS_H
#define OPENTISSUE_CORE_MATH_OPTIMIZATION_BFGS_H
//
// OpenTissue Template Library
// - A generic toolbox for physics-based modeling and simulation.
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
//
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
//
#include <OpenTissue/configuration.h>

#include <OpenTissue/core/math/big/big_types.h>
#include <OpenTissue/core/math/optimization/optimization_constants.h>
#include <OpenTissue/core/math/optimization/optimization_armijo_backtracking.h>
#include <OpenTissue/core/math/optimization/optimization_stationary_point.h>

#include <OpenTissue/core/math/math_value_traits.h>
#include <OpenTissue/core/math/math_is_number.h>

#include <cmath>
#include <stdexcept>
#include <cassert>

namespace OpenTissue
{
  namespace math
  {
    namespace optimization
    {

      namespace detail
      {

        template<typename vector_type,typename matrix_type>
        inline void bfgs_update_inverse_hessian(vector_type const & y,vector_type const & s, matrix_type & H)
        {
          using std::fabs;

          typedef typename matrix_type::value_type                     real_type;
          typedef typename ublas::identity_matrix<real_type>           identity_matrix_type;
          typedef typename OpenTissue::math::ValueTraits<real_type>    value_traits;

          if(H.size1() <= 0 || H.size2() <= 0)
            throw std::invalid_argument("bfgs_update_inverse_hessian(): Hessian was empty");
          if(y.size() <= 0 )
            throw std::invalid_argument("bfgs_update_inverse_hessian(): y_i was empty");
          if(s.size() <= 0 )
            throw std::invalid_argument("bfgs_update_inverse_hessian(): s_i was empty");

          assert(y.size()==H.size2() || "bfgs_update_inverse_hessian(): inconsistent matrix vector dimensions");

          assert(H.size1()==H.size2() || "bfgs_update_inverse_hessian(): Incompatible dimensions of Hessian matrix");

          size_t const N = H.size1();

          real_type dot = inner_prod(y,s);
          if(  fabs(dot) <= value_traits::zero() ) // Make sure we do not divide by zero
          {            
            dot = value_traits::one();
          }

          real_type rho = value_traits::one() / dot;
          
          //matrix_type A,B,C;          
          //A.resize(N,N,false);
          //B.resize(N,N,false);
          //C.resize(N,N,false);
          //identity_matrix_type I(N,N);
          //A = I - outer_prod(s,y)*rho;
          //B = I - outer_prod(y,s)*rho;
          //C = outer_prod(s,s)*rho;
          //matrix_type tmp;
          //tmp.resize(N,N,false);
          //ublas::noalias(tmp) = ublas::sparse_prod<matrix_type>(H,B);
          //ublas::noalias(H)   = ublas::sparse_prod<matrix_type>(A,tmp) + C;

          //
          //The update formula looks like this
          //
          //\f[H' = (I - rho s y^T) H(I - rho  y s^T) + rho s s^T\f]
          //
          //doing some manipulation one gets
          //
          //\f[H' = H + rho^2  (y^T H y) s s^T + rho s s^T - rho s y^T H   - rho H  y s^T\f]
          //
          //\f$H\f$ is symmetric postive definite so \f$(H y)^T = y^Y H^T = y^T H\f$ and
          //
          //\f[H' = H + (\rho^2  (y^T H y)+\rho) s s^T  - \rho ( s (H y)^T   +  (H  y)  s^T )\f]
          //
          //Introducing \f$v = H y\f$ and \f$\gamma = (\rho^2  (y^T H y)+rho)\f$
          //(note \f$y^T H y = y^T v\f$) this reads
          //
          //\f[H' = H + \gamma (s s^T)  - \rho ( s v^T   + v s^T )\f]
          //
          //Doing one more re-write (just a matter of parentheses) one has
          //
          //\f[H' = H + (\gamma s - \rho v) s^T - (\rho s) v^T \f]
          //
          //This can be efficiently implemented as a succesion of two rank-1 updates
          //
          vector_type v;
          v.resize(N,false);
          ublas::axpy_prod(H,y,v,true); // v = H y
          real_type gamma = rho*rho*ublas::inner_prod(y,v) + rho;
          for (typename vector_type::size_type i = 0; i < N; ++ i)
          {
            ublas::row(H,i) -= rho*s(i) * v;
            ublas::row(H,i) += (gamma*s(i) - rho * v(i)) * s;
          } 

        }

      }

      template < 
        typename T
        , typename function_functor
        , typename gradient_functor
      >
      inline void bfgs(
      function_functor  & f
      , gradient_functor & nabla_f
      , ublas::compressed_matrix<T>  & H
      , boost::numeric::ublas::vector<T> & x
      , size_t const & max_iterations
      , T      const & absolute_tolerance
      , T      const & relative_tolerance
      , T      const & stagnation_tolerance
      , size_t       & status
      , size_t       & iteration
      , T            & error
      , T      const & alpha
      , T      const & beta
      , ublas::vector<T> * profiling = 0
      )  
      {
        using std::fabs;
        using std::min;
        using std::max;

        typedef          ublas::compressed_matrix<T>       matrix_type;
        typedef          ublas::vector<T>                  vector_type;
        typedef          T                                 real_type;
        typedef          OpenTissue::math::ValueTraits<T>  value_traits;

        if(max_iterations <= 0)
          throw std::invalid_argument("max_iterations must be larger than zero");
        if(absolute_tolerance < value_traits::zero() )
          throw std::invalid_argument("absolute_tolerance must be non-negative");
        if(relative_tolerance < value_traits::zero() )
          throw std::invalid_argument("relative_tolerance must be non-negative");
        if(stagnation_tolerance < value_traits::zero() )
          throw std::invalid_argument("stagnation_tolerance must be non-negative");
        if (beta >= value_traits::one() )
          throw std::invalid_argument("Illegal beta value");
        if (alpha <= value_traits::zero() )
          throw std::invalid_argument("Illegal alpha value");
        if(beta<=alpha)
          throw std::invalid_argument("beta must be larger than alpha");
        if(profiling == &x)
          throw std::logic_error("profiling must not point to x-vector");

        error     = value_traits::infinity();
        iteration = 0;
        status    = OK;

        size_t const m = x.size();
        if(m==0)
          return;

        status = ITERATING; // Indicate that we are iterating and have not converged

        if(profiling)
        {
          (*profiling).resize( max_iterations );
          (*profiling).clear();
        }

        // Declare temporary storage
        vector_type y_k;
        vector_type s_k;
        vector_type dx;
        vector_type x_old;
        vector_type nabla_f_k;
        vector_type nabla_f_k1;

        // Allocate space for temporaries
        y_k.resize(m);
        s_k.resize(m);
        dx.resize(m);
        x_old.resize(m);
        nabla_f_k.resize(m);
        nabla_f_k1.resize(m);

        // Initialize 
        real_type             f_0   = f(x);
        ublas::noalias( nabla_f_k ) = nabla_f(x);

        // Iterate until convergence
        for (; iteration < max_iterations; ++iteration)
        {
          if(profiling)
            (*profiling)(iteration) = f_0;

          // Check for absolute convergence
          if(stationary_point( nabla_f_k, absolute_tolerance, error ) )
          {
            status = ABSOLUTE_CONVERGENCE;
            return;
          }

          // solve Newton system. Compute Search Direction
          // ublas::noalias( dx ) = - ublas::prod(H, nabla_f_k);
          ublas::axpy_prod(H, -nabla_f_k, dx, true);

          x_old.assign( x );
          real_type f_tau = f_0;
          armijo_backtracking(
            f
            , nabla_f_k
            , x_old
            , x
            , dx
            , relative_tolerance
            , stagnation_tolerance
            , alpha
            , beta
            , f_tau
            , status 
            );

          if(status != OK)
            return;

          // Do the incremental update of the inverse Hessian approximation
          ublas::noalias( nabla_f_k1 ) = nabla_f(x);
          ublas::noalias( s_k )  = x - x_old;
          ublas::noalias( y_k ) = nabla_f_k1 - nabla_f_k;

          detail::bfgs_update_inverse_hessian(y_k,s_k,H);

          // Update values for next iteration
          f_0 = f_tau;
          nabla_f_k.assign( nabla_f_k1 );

        }//end for loop
      }

    } // namespace optimization
  } // namespace math
} // namespace OpenTissue

// OPENTISSUE_CORE_MATH_OPTIMIZATION_BFGS_H
#endif