Rosenbrock.hpp

Go to the documentation of this file.

#ifndef HOPS_ROSENBROCK_HPP
#define HOPS_ROSENBROCK_HPP
 
#include <cmath>
#include <stdexcept>
#include <unsupported/Eigen/MatrixFunctions>
#include <utility>
 
#include <hops/Model/Model.hpp>
#include <hops/Utility/MatrixType.hpp>
#include <hops/Utility/VectorType.hpp>
 
namespace hops {
 
    class Rosenbrock : public Model {
    public:
        Rosenbrock(double scaleParameter, VectorType shiftParameter);
 
        [[nodiscard]] typename MatrixType::Scalar computeNegativeLogLikelihood(const VectorType &x) const override;
 
        [[nodiscard]] MatrixType computeHessian(const VectorType &x) const;
 
        [[nodiscard]] std::optional<VectorType> computeLogLikelihoodGradient(const VectorType &x) const override;
 
        [[nodiscard]] std::optional<MatrixType> computeExpectedFisherInformation(const VectorType &x) const override;
 
    private:
        typename MatrixType::Scalar scaleParameter;
        VectorType shiftParameter;
        long numberOfDimensions;
    };
 
    Rosenbrock::Rosenbrock(double scaleParameter,
                           VectorType shiftParameter) :
            scaleParameter(scaleParameter),
            shiftParameter(std::move(shiftParameter)) {
        numberOfDimensions = this->shiftParameter.rows() * 2;
    }
 
    typename MatrixType::Scalar
    Rosenbrock::computeNegativeLogLikelihood(const VectorType &x) const {
        if (x.rows() != numberOfDimensions) {
            throw std::runtime_error("Input x has wrong number of rows.");
        }
        typename MatrixType::Scalar result = 0;
        for (long i = 0; i < shiftParameter.rows(); ++i) {
            result += scaleParameter *
                      (100 * std::pow(std::pow(x(2 * i), 2) - x(2 * i + 1), 2) +
                       std::pow(x(2 * i) - shiftParameter(i), 2));
        }
        return result;
    }
 
    MatrixType
    Rosenbrock::computeHessian(const VectorType &x) const {
        MatrixType observedFisherInformation(x.rows(), x.rows());
 
        for (long i = 0; i < shiftParameter.rows(); ++i) {
            observedFisherInformation(2 * i, 2 * i) =
                    scaleParameter * (1200 * std::pow(x(2 * i), 2) - 400 * x(2 * i + 1) + 2);
            observedFisherInformation(2 * i + 1, 2 * i) = scaleParameter * -400 * x(2 * i);
            observedFisherInformation(2 * i, 2 * i + 1) = scaleParameter * -400 * x(2 * i);
            observedFisherInformation(2 * i + 1, 2 * i + 1) = scaleParameter * 200;
        }
 
        return observedFisherInformation;
    }
 
    std::optional<VectorType>
    Rosenbrock::computeLogLikelihoodGradient(const VectorType &x) const {
        VectorType gradient(x.rows());
        for (long i = 0; i < shiftParameter.rows(); ++i) {
            gradient(2 * i) = scaleParameter * (4 * 100 * (x(2 * i + 1) - std::pow(x(2 * i), 2)) * (-2 * x(2 * i)) +
                                                2 * (x(2 * i) - shiftParameter(i)));
 
            gradient(2 * i + 1) = scaleParameter * 2 * 100 * (x(2 * i + 1) - std::pow(x(2 * i), 2));
        }
        return gradient;
    }
 
    std::optional<MatrixType>
    Rosenbrock::computeExpectedFisherInformation(const VectorType &x) const {
        MatrixType hessian = computeHessian(x);
        // regularization should be between 0 and infinity. This value is guessed for now.
        double regularization = 1. / hessian.maxCoeff();
        MatrixType expPositive = (regularization * hessian).exp();
        MatrixType expNegative = (-regularization * hessian).exp();
        hessian = (expPositive + expNegative) * hessian * (expPositive - expNegative).inverse();
        return
                hessian;
    }
}
 
#endif //HOPS_ROSENBROCK_HPP