multi_extrema.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_MULTI_EXTREMA_H
#define THEORETICA_MULTI_EXTREMA_H
 
#include "../core/constants.h"
#include "../autodiff/autodiff.h"
#include "./extrema.h"
 
#include <functional>
 
 
namespace theoretica {
 
 
    template<unsigned int N>
    inline vec<real, N> multi_minimize_grad(
        multidual<N>(*f)(vec<multidual<N>, N>),
        vec<real, N> guess = vec<real, N>(0),
        real gamma = OPTIMIZATION_MINGRAD_GAMMA,
        real tolerance = OPTIMIZATION_MINGRAD_TOLERANCE,
        unsigned int max_iter = OPTIMIZATION_MINGRAD_ITER) {
 
        if(gamma >= 0) {
            TH_MATH_ERROR("multi_minimize_grad", gamma, INVALID_ARGUMENT);
            return vec<real, N>(nan());
        }
 
        vec<real, N> x = guess;
        vec<real, N> grad;
        unsigned int iter = 0;
 
        do {
 
            grad = gradient(f, x);
            x += gamma * grad;
            iter++;
 
        } while(grad.norm() > OPTIMIZATION_MINGRAD_TOLERANCE && iter <= max_iter);
 
        if(iter > max_iter) {
            TH_MATH_ERROR("multi_minimize_grad", iter, NO_ALGO_CONVERGENCE);
            return vec<real, N>(nan());
        }
 
        return x;
    }
 
 
    template<unsigned int N>
    inline vec<real, N> multi_maximize_grad(
        multidual<N>(*f)(vec<multidual<N>, N>),
        vec<real, N> guess = vec<real, N>(0),
        real gamma = OPTIMIZATION_MINGRAD_GAMMA,
        real tolerance = OPTIMIZATION_MINGRAD_TOLERANCE,
        unsigned int max_iter = OPTIMIZATION_MINGRAD_ITER) {
 
        return multi_minimize_grad(f, guess, -gamma, tolerance, max_iter);
    }
 
 
    template<unsigned int N>
    inline vec<real, N> multi_minimize_lingrad(
        multidual<N>(*f)(vec<multidual<N>, N>),
        vec<real, N> guess = vec<real, N>(0),
        real tolerance = OPTIMIZATION_MINGRAD_TOLERANCE,
        unsigned int max_iter = OPTIMIZATION_MINGRAD_ITER) {
 
        vec<real, N> x = guess;
        vec<real, N> grad;
        unsigned int iter = 0;
 
        do {
 
            grad = autodiff::gradient(f, x);
 
            // Minimize f(x + gamma * gradient) in [-1, 0]
            // using Golden Section extrema search
            real gamma = minimize_goldensection(
                [f, x, grad](real gamma){
                    return f(multidual<N>::make_argument(x + gamma * grad)).Re();
                }, -1, 0);
 
            // Fallback value
            if(-gamma <= MACH_EPSILON)
                gamma = OPTIMIZATION_MINGRAD_GAMMA;
 
            x += gamma * grad;
            iter++;
 
        } while(grad.norm() > OPTIMIZATION_MINGRAD_TOLERANCE && iter <= max_iter);
 
        if(iter > max_iter) {
            TH_MATH_ERROR("multi_minimize_lingrad", iter, NO_ALGO_CONVERGENCE);
            return vec<real, N>(nan());
        }
 
        return x;
    }
 
 
    template<unsigned int N>
    inline vec<real, N> multi_maximize_lingrad(
        multidual<N>(*f)(vec<multidual<N>, N>),
        vec<real, N> guess = vec<real, N>(0),
        real tolerance = OPTIMIZATION_MINGRAD_TOLERANCE,
        unsigned int max_iter = OPTIMIZATION_MINGRAD_ITER) {
 
        vec<real, N> x = guess;
        vec<real, N> grad;
        unsigned int iter = 0;
 
        do {
 
            grad = -gradient(f, x);
 
            // Maximize f(x + gamma * gradient) in [-1, 0]
            // using Golden Section extrema search
            real gamma = maximize_goldensection(
                [f, x, grad](real gamma){
                    return f(multidual<N>::make_argument(x + gamma * grad)).Re();
                }, -1, 0);
 
            // Fallback value
            if(-gamma <= MACH_EPSILON)
                gamma = OPTIMIZATION_MINGRAD_GAMMA;
 
            x += gamma * grad;
            iter++;
 
        } while(grad.norm() > OPTIMIZATION_MINGRAD_TOLERANCE && iter <= max_iter);
 
        if(iter > max_iter) {
            TH_MATH_ERROR("multi_maximize_lingrad", iter, NO_ALGO_CONVERGENCE);
            return vec<real, N>(nan());
        }
 
        return x;
    }
 
 
    template<unsigned int N>
    inline vec<real, N> multi_minimize(
        multidual<N>(*f)(vec<multidual<N>, N>),
        vec<real, N> guess = vec<real, N>(0), real tolerance = OPTIMIZATION_MINGRAD_TOLERANCE) {
 
        return multi_minimize_lingrad(f, guess, tolerance);
    }
 
 
    template<unsigned int N>
    inline vec<real, N> multi_maximize(
        multidual<N>(*f)(vec<multidual<N>, N>),
        vec<real, N> guess = vec<real, N>(0), real tolerance = OPTIMIZATION_MINGRAD_TOLERANCE) {
 
        return multi_maximize_lingrad<N>(f, guess, tolerance);
    }
 
}
 
 
#endif