實現最速下降算法，可變步長

Question

我想在編程語言（C/C++/fortran）中實現最速下降算法。

與F的例如最小化（X1，X2）= X1^3 + X 2^3 - 2 * X1 * X2

1. 估計開始設計點X0，迭代計數器K0，收斂參數耐性= 0.1 。 （1,0）

2. 將當前點x（k）處f（x1，x2）的梯度計算爲grad（f）。我將在這裏使用數字區分。

   d/DX1（F）= LIM（H-> 0）（F（X1 + H，X 2） - F（X1，X2））/ h的

   這是研究所（F）=（3 * X1^2 - 2 * X 2，3 * X 2^2 - 2 * X1）

   在（0,1研究所（F））是C0 =（3，-2）

3. 因爲L2範數C0>耐性的，我們繼續進行下一步驟

4. 方向D0 = -c0 =（-3,2）

5. 計算步長a。最小化f（a）= f（x0 + a d0）=（1-3a，2a）=（1-3a）^ 3 +（2a）^ 3 - 2（1-3a）*（2a）。我不保持步長不變。

6. update：new [x1，x2] = old [x1，x2] x + a * d0。

我不明白怎麼做步驟5 我有二分法一維最小化程序，它看起來像：

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中使用符號實現，但我不想使用符號。 任何建議表示讚賞。謝謝。

Answer 1

如果C++是一個選項，您可以利用functors和lambdas。

讓我們考慮我們要最小化的功能，例如Y = X - X + 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問題請求的第五步驟使用單維優化algorith以最小化沿着梯度的方向上的第一功能這需要傳遞一個單維函數。我們可以解決使用Lambda這個問題：

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。