newton_raphson.h

Go to the documentation of this file.

/*
 * Copyright 2021 MusicScience37 (Kenta Kabashima)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#pragma once
 
#include <cmath>
#include <limits>
#include <type_traits>
 
#include <Eigen/LU>
 
#include "num_collect/base/exception.h"
#include "num_collect/base/index_type.h"
#include "num_collect/base/precondition.h"
#include "num_collect/logging/iterations/iteration_logger.h"
#include "num_collect/logging/log_tag_view.h"
#include "num_collect/logging/logging_macros.h"
#include "num_collect/roots/concepts/differentiable_function.h"
#include "num_collect/roots/concepts/multi_variate_differentiable_function.h"
#include "num_collect/roots/concepts/single_variate_differentiable_function.h"
#include "num_collect/roots/function_root_finder_base.h"
 
namespace num_collect::roots {

constexpr auto newton_raphson_tag =
    logging::log_tag_view("num_collect::roots::newton_raphson");

template <concepts::differentiable_function Function>
class newton_raphson;

template <concepts::single_variate_differentiable_function Function>
class newton_raphson<Function>
    : public function_root_finder_base<newton_raphson<Function>, Function> {
public:
    using this_type = newton_raphson<Function>;

    using base_type =
        function_root_finder_base<newton_raphson<Function>, Function>;
 
    using typename base_type::function_type;
    using typename base_type::variable_type;
 
    static_assert(
        std::is_same_v<variable_type, typename function_type::jacobian_type>);

    explicit newton_raphson(const function_type& function = function_type())
        : base_type(newton_raphson_tag, function) {}

    void init(const variable_type& variable) {
        variable_ = variable;
        last_change_ = std::numeric_limits<variable_type>::infinity();
        iterations_ = 0;
        evaluations_ = 0;
 
        function().evaluate_on(variable_);
        ++evaluations_;
        using std::abs;
        value_norm_ = abs(function().value());
    }

    void iterate() {
        change_ = -function().value() / function().jacobian();
        variable_ += change_;
 
        function().evaluate_on(variable_);
        ++evaluations_;
        ++iterations_;
        using std::abs;
        last_change_ = abs(change_);
        value_norm_ = abs(function().value());
    }

    [[nodiscard]] auto is_stop_criteria_satisfied() const -> bool {
        return (iterations() > max_iterations_) ||
            (last_change() < tol_last_change_) ||
            (value_norm() < tol_value_norm_);
    }

    void configure_iteration_logger(
        logging::iterations::iteration_logger<this_type>& iteration_logger)
        const {
        iteration_logger.template append<index_type>(
            "Iter.", &this_type::iterations);
        iteration_logger.template append<index_type>(
            "Eval.", &this_type::evaluations);
        iteration_logger.template append<variable_type>(
            "Value", &this_type::value_norm);
        iteration_logger.template append<variable_type>(
            "Change", &this_type::last_change);
    }
 
    using base_type::function;

    [[nodiscard]] auto variable() const -> const variable_type& {
        return variable_;
    }

    [[nodiscard]] auto value() const
        -> std::invoke_result_t<decltype(&function_type::value),
            const function_type> {
        return function().value();
    }

    [[nodiscard]] auto jacobian() const
        -> std::invoke_result_t<decltype(&function_type::jacobian),
            const function_type> {
        return function().jacobian();
    }

    [[nodiscard]] auto iterations() const noexcept -> index_type {
        return iterations_;
    }

    [[nodiscard]] auto evaluations() const noexcept -> index_type {
        return evaluations_;
    }

    [[nodiscard]] auto last_change() const noexcept -> variable_type {
        return last_change_;
    }

    [[nodiscard]] auto value_norm() const noexcept -> variable_type {
        return value_norm_;
    }

    auto max_iterations(index_type val) -> this_type& {
        NUM_COLLECT_PRECONDITION(val > 0, this->logger(),
            "Maximum number of iterations must be a positive integer.");
        max_iterations_ = val;
        return *this;
    }

    auto tol_last_change(const variable_type& val) -> this_type& {
        NUM_COLLECT_PRECONDITION(val >= static_cast<variable_type>(0),
            this->logger(),
            "Tolerance of last change of the variable must be a "
            "non-negative value.");
        tol_last_change_ = val;
        return *this;
    }

    auto tol_value_norm(const variable_type& val) -> this_type& {
        NUM_COLLECT_PRECONDITION(val >= static_cast<variable_type>(0),
            this->logger(),
            "Tolerance of the norm of function value must be a "
            "non-negative value.");
        tol_value_norm_ = val;
        return *this;
    }
 
private:
    variable_type variable_{};

    variable_type change_{};

    index_type iterations_{};

    index_type evaluations_{};

    variable_type last_change_{};

    variable_type value_norm_{};

    static constexpr index_type default_max_iterations = 1000;

    index_type max_iterations_{default_max_iterations};

    static constexpr auto default_tol_last_change =
        static_cast<variable_type>(1e-6);

    variable_type tol_last_change_{default_tol_last_change};

    static constexpr auto default_tol_value_norm =
        static_cast<variable_type>(1e-6);

    variable_type tol_value_norm_{default_tol_value_norm};
};

template <concepts::multi_variate_differentiable_function Function>
class newton_raphson<Function>
    : public function_root_finder_base<newton_raphson<Function>, Function> {
public:
    using this_type = newton_raphson<Function>;

    using base_type =
        function_root_finder_base<newton_raphson<Function>, Function>;
 
    using typename base_type::function_type;
    using typename base_type::variable_type;

    using scalar_type = typename variable_type::Scalar;

    using jacobian_type = typename function_type::jacobian_type;

    using jacobian_solver_type =
        Eigen::PartialPivLU<typename Function::jacobian_type>;

    explicit newton_raphson(const function_type& function = function_type())
        : base_type(newton_raphson_tag, function) {}

    void init(const variable_type& variable) {
        variable_ = variable;
        last_change_ = std::numeric_limits<scalar_type>::infinity();
        iterations_ = 0;
        evaluations_ = 0;
 
        function().evaluate_on(variable_);
        ++evaluations_;
        value_norm_ = function().value().norm();
    }

    void iterate() {
        jacobian_solver_.compute(function().jacobian());
        change_ = -jacobian_solver_.solve(function().value());
        variable_ += change_;
 
        function().evaluate_on(variable_);
        ++evaluations_;
        ++iterations_;
        last_change_ = change_.norm();
        value_norm_ = function().value().norm();
    }

    [[nodiscard]] auto is_stop_criteria_satisfied() const -> bool {
        return (iterations() > max_iterations_) ||
            (last_change() < tol_last_change_) ||
            (value_norm() < tol_value_norm_);
    }

    void configure_iteration_logger(
        logging::iterations::iteration_logger<this_type>& iteration_logger)
        const {
        iteration_logger.template append<index_type>(
            "Iter.", &this_type::iterations);
        iteration_logger.template append<index_type>(
            "Eval.", &this_type::evaluations);
        iteration_logger.template append<scalar_type>(
            "Value", &this_type::value_norm);
        iteration_logger.template append<scalar_type>(
            "Change", &this_type::last_change);
    }
 
    using base_type::function;

    [[nodiscard]] auto variable() const -> const variable_type& {
        return variable_;
    }

    [[nodiscard]] auto value() const
        -> std::invoke_result_t<decltype(&function_type::value),
            const function_type> {
        return function().value();
    }

    [[nodiscard]] auto jacobian() const
        -> std::invoke_result_t<decltype(&function_type::jacobian),
            const function_type> {
        return function().jacobian();
    }

    [[nodiscard]] auto iterations() const noexcept -> index_type {
        return iterations_;
    }

    [[nodiscard]] auto evaluations() const noexcept -> index_type {
        return evaluations_;
    }

    [[nodiscard]] auto last_change() const noexcept -> scalar_type {
        return last_change_;
    }

    [[nodiscard]] auto value_norm() const noexcept -> scalar_type {
        return value_norm_;
    }

    auto max_iterations(index_type val) -> this_type& {
        NUM_COLLECT_PRECONDITION(val > 0, this->logger(),
            "Maximum number of iterations must be a positive integer.");
        max_iterations_ = val;
        return *this;
    }

    auto tol_last_change(const scalar_type& val) -> this_type& {
        NUM_COLLECT_PRECONDITION(val >= static_cast<scalar_type>(0),
            this->logger(),
            "Tolerance of last change of the variable must be a "
            "non-negative value.");
        tol_last_change_ = val;
        return *this;
    }

    auto tol_value_norm(const scalar_type& val) -> this_type& {
        NUM_COLLECT_PRECONDITION(val >= static_cast<scalar_type>(0),
            this->logger(),
            "Tolerance of the norm of function value must be a "
            "non-negative value.");
        tol_value_norm_ = val;
        return *this;
    }
 
private:
    variable_type variable_{};

    variable_type change_{};

    jacobian_solver_type jacobian_solver_{};

    index_type iterations_{};

    index_type evaluations_{};

    scalar_type last_change_{};

    scalar_type value_norm_{};

    static constexpr index_type default_max_iterations = 1000;

    index_type max_iterations_{default_max_iterations};

    static constexpr auto default_tol_last_change =
        static_cast<scalar_type>(1e-6);

    scalar_type tol_last_change_{default_tol_last_change};

    static constexpr auto default_tol_value_norm =
        static_cast<scalar_type>(1e-6);

    scalar_type tol_value_norm_{default_tol_value_norm};
};
 
}  // namespace num_collect::roots