实现最速下降算法，可变步长

Implementing steepest descent algorithm, variable step size

本文关键字：算法 实现      更新时间：2023-10-16

我正在尝试用编程语言(C/C++/fortran)实现最速下降算法。

例如，f(x1，x2)的最小化=x1^3+x2^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，x2)-f(x1，x2))/h

    这是梯度(f)=(3*x1^2-2*x2，3*x2^2-2*x1)

    (0,1)处的grad(f)为c0=(3，-2)

  3. 由于c0>公差的L2范数，我们进行下一步

  4. 方向d0=-c0=(-3,2)

  5. 计算步长a。最小化f(a)=f(x0+ad0)=(1-3a，2a)=(1-3 a)^3+(2a)^3-2(1-3a)*(2a)。我没有保持恒定的步长。

  6. update:new[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中使用符号来实现，但我不想使用符号。欢迎提出任何建议。谢谢

如果C++是一个选项，则可以利用函子和lambda。

让我们考虑一个我们想要最小化的函数，例如y=x2-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问题的第五步要求使用需要传递一维函数的一维优化算法来最小化沿梯度方向的第一个函数。我们可以使用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);

请在这里举个例子。

相关文章：