Реализация алгоритма наивысшего спуска, изменяемый размер шага

Question

Я пытаюсь реализовать алгоритм наименьшего спуска в языках программирования (C/C++/fortran).

К примеру минимизации F (x1, x2) = x1^3 + x2^3 - 2 * х1 * х2

1. Расчетный начиная расчетной точке x0, итерационный счетчик k0, параметр tolerence сходимости = 0,1 , Скажем, что эта точка зрения (1,0)

2. Вычислить градиент f (x1, x2) в текущей точке x (k) как grad (f). Здесь я буду использовать численное дифференцирование.

   д/дх1 (е) = Пт (h-> 0) (F (x 1 + ч, х2) - F (x1, x2))/ч

   Это градиента (е) = (3 * x1^2 - 2 * х2, 3 * х2^2 - 2 * x1)

   град (е) в точке (0,1) с0 = (3, -2)

3. , поскольку норма L2 с0> tolerence, мы переходим к следующему шагу

4. направление d0 = с0 = (-3,2)

5. Вычислить размер шага a. Минимизировать f (a) = f (x0 + a d0) = (1-3a, 2a) = (1-3a)^3 + (2a)^3 - 2 (1-3a) * (2a). Я не придерживаюсь постоянного размера шага.

6. обновление: новый [x1, x2] = старый [x1, x2] x + a * d0.

Я не понимаю, как сделать шаг 5. У меня есть программа минимизации 1D с методом бисекций, и это выглядит следующим образом:

program main() 
    ... 
    ... 
    define upper, lower interval 
    call function value 
    ...calculations 
    ... 
    ... 


function value (input x1in) (output xout) 
    ...function is x^4 - 2x^2 + x + 10 
    xout = (xin)^4 - 2*(xin)^2 + (xin) + 10 

В этом случае, глядя на шаге 5, я не может проходить символом a. Любые идеи о том, как реализовать алгоритм в языке программирования, особенно на этапе 5? Пожалуйста, предложите, если есть совсем другой способ программирования этого. Я видел много программ с постоянным размером шага, но я хочу вычислить его на каждом шагу. Этот алгоритм может быть легко реализован в MATLAB ot python sympy с использованием символики, но я не хочу использовать символику. Любые предложения оценены. Благодарю.

Answer 1

Если C++ - это опция, вы можете использовать functors и lambdas.

Давайте рассмотрим функцию, мы хотим, чтобы свести к минимуму, например, у = х - х + 2. Он может быть представлен в виде функции объекта, который является классом с перегруженной operator():

struct MyFunc { 
    double operator()(double x) const { 
     return x * x - x + 2.0; 
    } 
}; 

Теперь мы можем объявить объект этого типа, использовать его как функцию и передать его в другую шаблонную функцию как шаблонный параметр.Функция

// given this templated function: 
template < typename F > 
void tabulate_function(F func, double a, double b, int steps) { 
    // the functor  ^^^^^^ is passed to the templated function 
    double step = (b - a)/(steps - 1); 

    std::cout << " x   f(x)\n------------------------\n"; 
    for (int i = 0; i < steps; ++i) { 
     double x = a + i * step, 
       fx = func(x); 
     //   ^^^^^^^ call the operator() of the functor 
     std::cout << std::fixed << std::setw(8) << std::setprecision(3) << x 
        << std::scientific << std::setw(16) << std::setprecision(5) 
        << fx << '\n'; 
    } 
} 

// we can use the previous functor like this: 
MyFunc example; 
tabulate_function(example, 0.0, 2.0, 21); 

OP может быть реализована (дана вспомогательный класс для представления 2D точек) аналогичным образом:

struct MyFuncVec { 
    double operator()(const Point &p) const { 
     return p.x * p.x * p.x + p.y * p.y * p.y - 2.0 * p.x * p.y; 
    } 
}; 

градиента этой функции можно представить (данный класс, который реализовать 2D вектор) по:

struct MyFuncGradient { 
    Vector operator()(const Point &p) { 
     return Vector(3.0 * p.x * p.x - 2.0 * p.y, 3.0 * p.y * p.y - 2.0 * p.x); 
    } 
};  

Теперь, пятый шаг запросов вопрос OP, чтобы свести к минимуму первой функции вдоль направления градиента с использованием monodimensional оптимизации Алгоритм Построения m, который требует передачи мономерной функции. Мы можем решить эту проблему с помощью лямбда:

MyFuncVec funcOP; 
    MyFuncGradient grad_funcOP; 
    Point p0(0.2, 0.8); 
    Vector g = grad_funcOP(p0); 

    // use a lambda to transform the OP function to 1D 
    auto sliced_func = [&funcOP, &p0, &g] (double t) -> double { 
     // those variables ^^^ ^^^ ^^ are captured and used 
     return funcOP(p0 - t * g); 
    }; 

    tabulate_function(sliced_func, 0, 0.5, 21); 

Живой пример HERE.