quasi_newton.hpp 2.22 KB
Newer Older
1
// Software License for MTL
2 3
//
// Copyright (c) 2007 The Trustees of Indiana University.
4
//               2008 Dresden University of Technology and the Trustees of Indiana University.
5
//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
6 7
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
8
//
9
// This file is part of the Matrix Template Library
10
//
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
// See also license.mtl.txt in the distribution.

#ifndef ITL_QUASI_NEWTON_INCLUDE
#define ITL_QUASI_NEWTON_INCLUDE

#include <boost/numeric/mtl/concept/collection.hpp>
#include <boost/numeric/mtl/matrix/dense2D.hpp>
#include <boost/numeric/mtl/matrix/operators.hpp>
#include <boost/numeric/mtl/operation/operators.hpp>
#include <boost/numeric/mtl/operation/trans.hpp>
#include <boost/numeric/mtl/utility/gradient.hpp>

// #include <iostream>

namespace itl {

/// Quasi-Newton method
28
template <typename Matrix, typename Vector, typename F, typename Grad,
29
	  typename Step, typename Update, typename Iter>
30 31
Vector quasi_newton(Vector& x, F f, Grad grad_f, Step step, Update update, Iter& iter)
{
32 33 34
    typedef typename mtl::Collection<Vector>::value_type value_type;
    Vector         d, y, x_k, s;
    Matrix         H(size(x), size(x));
35

36 37
    H= 1;
    for (; !iter.finished(two_norm(grad_f(x))); ++iter) {
38
	d= H * -grad_f(x);                                                  // std::cout << "d is " << d << '\n';
39 40 41 42
	value_type alpha= step(x, d, f, grad_f); assert(alpha == alpha);
	x_k= x + alpha * d;                                                 // std::cout << "x_k is " << x_k << '\n';
	s= alpha * d;                                                       // std::cout << "alpha is " << alpha << '\n';
	y= grad_f(x_k) - grad_f(x);
43 44
	update(H, y, s);
	x= x_k;
45 46 47 48 49 50
    }
    return x;
}

/// Quasi-Newton method
template <typename Vector, typename F, typename Grad, typename Step, typename Update, typename Iter>
51
Vector inline quasi_newton(Vector& x, F f, Grad grad_f, Step step, Update update, Iter& iter)
52 53 54 55 56 57 58 59 60 61
{
    typedef typename mtl::traits::gradient<Vector>::type hessian_type;
    // typedef typename mtl::Collection<Vector>::value_type value_type;
    return quasi_newton<hessian_type>(x, f, grad_f, step, update, iter);
}


} // namespace itl

#endif // ITL_QUASI_NEWTON_INCLUDE