当前位置：首页 > news >正文

[最优化理论] 梯度下降法 + 精确线搜索（单峰区间搜索 + 黄金分割）C++ 代码

news 来源：原创 2024/5/20 9:07:03

这是我的课程作业，用了 Eigen 库，最后的输出是 latex 的表格的一部分

具体内容就是梯度下降法 + 精确线搜索（单峰区间搜索 + 黄金分割）

从书本的 Matlab 代码转译过来的其实，所以应该是一看就懂了

这里定义了两个测试函数 fun 和 fun2

整个最优化方法包装在 SteepestDescent 类里面

用了模板封装类，这样应该是 double 和 Eigne 的 Vector 都可以支持的

用了 tuple 返回值，用了 functional 接受函数形参，所以应该要 C++11 以上进行编译

#include "3rdparty/Eigen/Eigen/Dense"#include <cstdint>
#include <fstream>
#include <functional>
#include <iostream>
#include <string>
#include <tuple>#ifndef DEBUG
#    define DEBUG 0
#endifusing namespace Eigen;template<class YClass, class XClass>
class SteepestDescent
{
public:SteepestDescent(std::function<YClass(XClass)> const& fun,std::function<XClass(XClass)> const& gfun,double                               delta,double                               epsilon): m_fun(fun), m_gfun(gfun), m_delta(delta), m_epsilon(epsilon) {};/*** @brief Find single peak interval.** It will stop if the number of iterations exceeds the given upper limit.** @param fun Target function.* @param alpha0 Start point.* @param h Search direction.** @return XClass Left end of single peak interval.* @return XClass Right end of single peak interval.* @return XClass Inner point of single peak interval.* 1 represents same direction w.r.t. h, -1 represents reversed direction w.r.t. h.*/std::tuple<XClass, XClass, XClass> ForwardBackward(XClass alpha0, XClass h);/*** @brief Find a minimum of a function inside a specified interval.** @param fun Target function.* @param a Left end of interval.* @param b Right end of interval.* @param delta Tolerable error of input variable.* @param epsilon Tolerable error of target function value.** @return bool Is early stop. Let interpolation points to be p, q, if fun(a) < fun(p) and fun(q) > fun(b)* @return XClass Minimum point.* @return YClass Function value of minimum point.*/std::tuple<bool, XClass, YClass> GoldenSectionSearch(XClass a, XClass b);/*** @brief Run Forward Backward and Golden Section Search** @param fun Target function.* @param gfun Gredient of target function.* @param x0 Start point.* @param h Search direction.* @param delta Tolerable error of input variable.* @param epsilon Tolerable error of target function value.* @return std::tuple<YClass, YClass, uint32_t>*/std::tuple<XClass, YClass, uint32_t> ForwardBackwardAndGoldenSectionSearch(XClass x0);/*** @brief Run Armijo Search** @param fun Target function.* @param gfun Gredient of target function.* @param x0 Start point.* @param h Search direction.* @param delta Tolerable error of input variable.* @param epsilon Tolerable error of target function value.* @return std::tuple<YClass, YClass, uint32_t>*/std::tuple<XClass, YClass, uint32_t> ArmijoSearch(XClass x0);private:std::function<YClass(XClass)> m_fun;std::function<XClass(XClass)> m_gfun;double                        m_delta;double                        m_epsilon;
};template<class YClass, class XClass>
std::tuple<XClass, XClass, XClass> SteepestDescent<YClass, XClass>::ForwardBackward(XClass alpha0, XClass h)
{uint32_t k = 0, max_k = 500;bool     reversed = false;XClass alpha1 = alpha0, alpha = alpha0;YClass phi0 = m_fun(alpha0), phi1 = m_fun(alpha0);double t = 1e-2;while (k < max_k){alpha1 = alpha0 + t * h;phi1   = m_fun(alpha1);// forward searchif (phi1 < phi0){t      = 2.0 * t;alpha  = alpha0;alpha0 = alpha1;phi0   = phi1;}else{// backward searchif (k == 0){t     = -t;alpha = alpha1;}// find another endelse{break;}}++k;}#if DEBUGstd::cout << "ForwardBackward total iteration = " << std::endl;std::cout << k << std::endl;
#endifXClass left  = t > 0.0 ? alpha : alpha1;XClass right = t < 0.0 ? alpha : alpha1;return {left, right, alpha0};
}template<class YClass, class XClass>
std::tuple<bool, XClass, YClass> SteepestDescent<YClass, XClass>::GoldenSectionSearch(XClass a, XClass b)
{uint32_t k = 0, max_k = 500;double t = (sqrt(5) - 1.0) / 2.0;XClass h = b - a;XClass p = a + (1 - t) * h, q = a + t * h;YClass phia = m_fun(a), phib = m_fun(b);YClass phip = m_fun(p), phiq = m_fun(q);bool is_early_stop = false;if (phia < phip && phiq > phib){is_early_stop = true;#if DEBUGstd::cout << "GoldenSectionSearch total it eration = " << std::endl;std::cout << k << std::endl;
#endifreturn {is_early_stop, a, phia};}while (((abs(phip - phia) > m_epsilon) || (h.norm() > m_delta)) && k < max_k){if (phip < phiq){b = q;q = p;phib = phiq;phiq = phip;h = b - a;p = a + (1 - t) * h;phip = m_fun(p);}else{a = p;p = q;phia = phip;phip = phiq;h = b - a;q = a + t * h;phiq = m_fun(q);}++k;}#if DEBUGstd::cout << "GoldenSectionSearch total iteration = " << std::endl;std::cout << k << std::endl;
#endifif (phip <= phiq){return {is_early_stop, p, phip};}else{return {is_early_stop, q, phiq};}
}template<class YClass, class XClass>
std::tuple<XClass, YClass, uint32_t> SteepestDescent<YClass, XClass>::ForwardBackwardAndGoldenSectionSearch(XClass x0)
{uint32_t k = 0, max_k = 5000;YClass phi_min = m_fun(x0);#if DEBUG// file pointerstd::fstream fout;// opens an existing csv file or creates a new file.fout.open("SteepestDescent.csv", std::ios::out | std::ios::trunc);// Insert the data to filefout << x0[0] << ", " << x0[1] << ", " << phi_min << "\n";
#endifwhile (k < max_k){Vector2d h = -m_gfun(x0);if (h.norm() < m_epsilon){return {x0, phi_min, k};}auto [left, right, inner] = ForwardBackward(x0, h);auto [is_early_stop, x1, phix1] = GoldenSectionSearch(left, right);if (is_early_stop){x1    = inner;phix1 = m_fun(x1);}x0      = x1;phi_min = phix1;++k;#if DEBUGstd::cout << "iteration " << k << ":" << std::endl;std::cout << "h = " << std::endl;std::cout << h << std::endl;std::cout << "left pointer = " << std::endl;std::cout << left << std::endl;std::cout << "right pointer = " << std::endl;std::cout << right << std::endl;std::cout << "inner pointer = " << std::endl;std::cout << inner << std::endl;std::cout << "current point = " << std::endl;std::cout << x1 << std::endl;std::cout << "current evaluation = " << std::endl;std::cout << phix1 << std::endl;// Insert the data to filefout << x0[0] << ", " << x0[1] << ", " << phi_min << "\n";
#endif}return {x0, phi_min, k};
}template<class YClass, class XClass>
std::tuple<XClass, YClass, uint32_t> SteepestDescent<YClass, XClass>::ArmijoSearch(XClass x0)
{uint32_t k = 0, max_k = 5000;YClass phi_min = m_fun(x0);double rho   = 0.5;double sigma = 0.4;while (k < max_k){Vector2d h = -m_gfun(x0);if (h.norm() < m_epsilon){return {x0, phi_min, k};}uint32_t m  = 0;uint32_t mk = 0;while (m < 20) // Armijo Search{phi_min = m_fun(x0 + pow(rho, m) * h);if (phi_min < m_fun(x0) + sigma * pow(rho, m) * (-pow(h.norm(), 2.0))){mk = m;break;}m = m + 1;}x0 = x0 + pow(rho, mk) * h;++k;}return {x0, phi_min, k};
}double fun(Vector2d x) { return 100.0 * pow(pow(x[0], 2.0) - x[1], 2.0) + pow(x[0] - 1, 2.0); }Vector2d gfun(Vector2d x)
{return Vector2d(400.0 * x[0] * (pow(x[0], 2.0) - x[1]) + 2.0 * (x[0] - 1.0), -200.0 * (pow(x[0], 2.0) - x[1]));
}double fun2(Vector2d x) { return 3.0 * pow(x[0], 2.0) + 2.0 * pow(x[1], 2.0) - 4.0 * x[0] - 6.0 * x[1]; }Vector2d gfun2(Vector2d x) { return Vector2d(6.0 * x[0] - 4.0, 4.0 * x[1] - 6.0); }int main()
{std::vector<Vector2d> points {Vector2d(0.0, 0.0),Vector2d(2.0, 1.0),Vector2d(1.0, -1.0),Vector2d(-1.0, -1.0),Vector2d(-1.2, 1.0),Vector2d(10.0, 10.0)};SteepestDescent<double, Vector2d> sd(fun, gfun, 1e-4, 1e-5);std::fstream fout_result_1, fout_result_2;fout_result_1.open("ForwardBackwardAndGoldenSectionSearch_Result.csv", std::ios::out | std::ios::trunc);fout_result_2.open("ArmijoSearch_Result.csv", std::ios::out | std::ios::trunc);fout_result_1 << "初始点 ($x_0$) & 目标函数值 ($f(x_k)$) & 迭代次数 ($k$) \\\\"<< "\n";fout_result_1 << "\\midrule"<< "\n";fout_result_2 << "初始点 ($x_0$) & 目标函数值 ($f(x_k)$) & 迭代次数 ($k$) \\\\"<< "\n";fout_result_2 << "\\midrule"<< "\n";for (size_t i = 0; i < points.size(); ++i){auto [x, val, k] = sd.ForwardBackwardAndGoldenSectionSearch(points[i]);fout_result_1 << "$(" << points[i][0] << ", " << points[i][1] << ")^T$ & " << val << " & " << k << " \\\\"<< "\n";auto [x2, val2, k2] = sd.ArmijoSearch(points[i]);fout_result_2 << "$(" << points[i][0] << ", " << points[i][1] << ")^T$ & " << val2 << " & " << k2 << " \\\\"<< "\n";}fout_result_1.close();fout_result_2.close();
}