diff_hessian.h

Go to the documentation of this file.

/* =========================================================================== *
   |
   |  Copyright (c) 1994-2011 by Kobus Barnard (author)
   |
   |  Personal and educational use of this code is granted, provided that this
   |  header is kept intact, and that the authorship is not misrepresented, that
   |  its use is acknowledged in publications, and relevant papers are cited.
   |
   |  For other use contact the author (kobus AT cs DOT arizona DOT edu).
   |
   |  Please note that the code in this file has not necessarily been adequately
   |  tested. Naturally, there is no guarantee of performance, support, or fitness
   |  for any particular task. Nonetheless, I am interested in hearing about
   |  problems that you encounter.
   |
   |  Author:  Ernesto Brau
 * =========================================================================== */

/* $Id$ */

#ifndef DIFF_HESSIAN_H
#define DIFF_HESSIAN_H

#include <m_cpp/m_matrix.h>
#include <diff_cpp/diff_util.h>
#include <vector>

namespace kjb {

template<class Func, class Model, class Adapter>
Matrix hessian
(
    const Func& f,
    const Model& x,
    const std::vector<double>& dx,
    const Adapter& adapter
)
{
    const size_t D = adapter.size(&x);

    Matrix H(D, D);
    double fx = f(x);
    Model y = x;

    for(size_t i = 0; i < D; i++)
    {
        for(size_t j = 0; j < D; j++)
        {
            if(i == j)
            {
                move_param(y, i, dx[i], adapter);
                double fxp = f(y);

                move_param(y, i, -2*dx[i], adapter);
                double fxm = f(y);

                H(i, i) = (fxp - 2*fx + fxm) / (dx[i]*dx[i]);

                move_param(y, i, dx[i], adapter);
            }
            else
            {
                move_params(y, i, j, dx[i], dx[j], adapter);
                double fxpp = f(y);

                move_param(y, j, -2*dx[j], adapter);
                double fxpm = f(y);

                move_param(y, i, -2*dx[i], adapter);
                double fxmm = f(y);

                move_param(y, j, 2*dx[j], adapter);
                double fxmp = f(y);

                H(i, j) = (fxpp - fxpm - fxmp + fxmm) / (4*dx[i]*dx[j]);

                move_params(y, i, j, dx[i], -dx[j], adapter);
            }
        }
    }

    return H;
}

template<class Func, class Vec>
inline
Matrix hessian
(
    const Func& f,
    const Vec& x,
    const std::vector<double>& dx
)
{
    return hessian(f, x, dx, Vector_adapter<Vec>());
}

template<class Func, class Model, class Adapter>
Matrix hessian_symmetric
(
    const Func& f,
    const Model& x,
    const std::vector<double>& dx,
    const Adapter& adapter
)
{
    const size_t D = adapter.size(&x);

    Matrix H(D, D);
    double fx = f(x);
    Model y = x;

    for(size_t i = 0; i < D; i++)
    {
        for(size_t j = 0; j <= i; j++)
        {
            if(i == j)
            {
                move_param(y, i, dx[i], adapter);
                double fxp = f(y);

                move_param(y, i, -2*dx[i], adapter);
                double fxm = f(y);

                H(i, i) = (fxp - 2*fx + fxm) / (dx[i]*dx[i]);

                move_param(y, i, dx[i], adapter);
            }
            else
            {
                move_params(y, i, j, dx[i], dx[j], adapter);
                double fxpp = f(y);

                move_param(y, j, -2*dx[j], adapter);
                double fxpm = f(y);

                move_param(y, i, -2*dx[i], adapter);
                double fxmm = f(y);

                move_param(y, j, 2*dx[j], adapter);
                double fxmp = f(y);

                // symmetry assumption
                H(i, j) = (fxpp - fxpm - fxmp + fxmm) / (4*dx[i]*dx[j]);
                H(j, i) = H(i, j);

                move_params(y, i, j, dx[i], -dx[j], adapter);
            }
        }
    }

    return H;
}

template<class Func, class Vec>
inline
Matrix hessian_symmetric
(
    const Func& f,
    const Vec& x,
    const std::vector<double>& dx
)
{
    return hessian_symmetric(f, x, dx, Vector_adapter<Vec>());
}

template<class Func, class Model, class Adapter>
Matrix hessian_diagonal
(
    const Func& f,
    const Model& x,
    const std::vector<double>& dx,
    const Adapter& adapter
)
{
    const size_t D = adapter.size(&x);

    Matrix H(D, D);
    double fx = f(x);
    Model y = x;

    for(size_t i = 0; i < D; i++)
    {
        for(size_t j = 0; j < D; j++)
        {
            if(i == j)
            {
                move_param(y, i, dx[i], adapter);
                double fxp = f(y);

                move_param(y, i, -2*dx[i], adapter);
                double fxm = f(y);

                H(i, i) = (fxp - 2*fx + fxm) / (dx[i]*dx[i]);

                move_param(y, i, dx[i], adapter);
            }
        }
    }
    return H;
}

template<class Func, class Vec>
inline
Matrix hessian_diagonal
(
    const Func& f,
    const Vec& x,
    const std::vector<double>& dx
)
{
    return hessian_diagonal(f, x, dx, Vector_adapter<Vec>());
}

} //namespace kjb

#endif /*DIFF_HESSIAN_H */