multi_roots.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_MULTI_ROOTS_H
#define THEORETICA_MULTI_ROOTS_H
 
#include "../autodiff/autodiff.h"
#include "../core/iter_result.h"
 
 
namespace theoretica {
 
 
    template <
        typename Vector = vec<real>,
        typename ReturnVector = Vector,
        typename DualObjectiveFunction,
        autodiff::enable_vector_field<DualObjectiveFunction> = true
    >
    inline iter_result<ReturnVector> multiroot_newton(
        DualObjectiveFunction f,
        Vector guess,
        real tolerance = OPTIMIZATION_MINGRAD_TOLERANCE,
        unsigned int max_iter = OPTIMIZATION_MINGRAD_ITER
    ) {
 
        // Current position
        ReturnVector x = guess;
 
        // Number of iterations
        unsigned int iter = 0;
 
        // Extract the size of the vector type
        constexpr size_t N = ReturnVector::size_argument;
 
        // The current value of f(x)
        ReturnVector f_x = multidual<N>::extract_real(
            f(multidual<N>::make_argument(x))
        );
 
        // Jacobian matrix
        mat<real, N, N> J;
 
        while(f_x.sqr_norm() > square(tolerance) && iter <= max_iter) {
 
            // Compute the function value and Jacobian
            // using automatic differentiation
            multidual<N>::extract(
                f(multidual<N>::make_argument(x)), f_x, J
            );
 
            // Update the current best guess
            x -= algebra::solve(J, f_x);
            iter++;
        }
 
        if(iter > max_iter) {
 
            TH_MATH_ERROR("multiroot_newton", iter, MathError::NoConvergence);
            algebra::vec_error(x);
            return iter_result<ReturnVector>(
                x, ConvergenceStatus::MaxIterations, iter, algebra::norm(f_x)
            );
        }
 
        return iter_result<ReturnVector>(x, iter, algebra::norm(f_x));
    }
 
}
 
#endif