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

这是我的课程作业,用了 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();
}

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/194483.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

使用pytorch从零开始实现迷你GPT

生成式建模知识回顾: [1] 生成式建模概述 [2] Transformer I&#xff0c;Transformer II [3] 变分自编码器 [4] 生成对抗网络&#xff0c;高级生成对抗网络 I&#xff0c;高级生成对抗网络 II [5] 自回归模型 [6] 归一化流模型 [7] 基于能量的模型 [8] 扩散模型 I, 扩散模型 II…

wordpress建站优化加速教程-Redis加速

这篇文章适合宝塔面板&#xff0c;在宝塔面板安装 Redis 实现网站加速&#xff08; Redis是一个高性能的key-value数据库(PHP连接redis&#xff0c;需PHP设置中安装redis扩展) &#xff09;。对在word press网站有着明显的加速效果。关于Redis具体说明请自己百度&#xff0c;…

编程中常见的技术难题有哪些?By AI

编程对于现代社会发展的重要性 编程&#xff0c;即按照特定的规则和逻辑&#xff0c;为计算机设计指令的过程&#xff0c;已经深深地融入现代社会的各个角落。它对人们的生活、工作和科技发展产生了深远的影响。 首先&#xff0c;编程改变了人们的生活方式。如今&#xff0c;…

Qt 如何操作SQLite3数据库?数据库创建和表格的增删改查?

# 前言 项目源码下载 https://gitcode.com/m0_45463480/QSQLite3/tree/main # 第一步 项目配置 平台:windows10 Qt版本:Qt 5.14.2 在.pro添加 QT += sql 需要的头文件 #include <QSqlDatabase>#include <QSqlError>#include <QSqlQuery>#include &…

SCAU:数字字符序列2

数字字符序列 Time Limit:1000MS Memory Limit:65535K 题型: 填空题 语言: G;GCC;VC 描述 有一个数字字符序列&#xff0c;它是由平方数1&#xff0c;4&#xff0c;9&#xff0c;16&#xff0c;25&#xff0c;36&#xff0c;49&#xff0c;64&#xff0c;81&#xff0c;…

Pandas进阶:拼接 concat 使用方法

1.处理索引和轴 假设我们有2个关于考试成绩的数据集。 df1 pd.DataFrame&#xff08;{ name&#xff1a;[A&#xff0c;B&#xff0c;C&#xff0c;D]&#xff0c;math&#xff1a;[60,89,82,70]&#xff0c;physics&#xff1a;[66&#xff0c; 95,83,66]&#xff0c;chemi…

SCAU:前一个和后一个字符

前一个和后一个字符 Time Limit:1000MS Memory Limit:65535K 题型: 编程题 语言: G;GCC;VC 描述 编写程序&#xff0c;输入一个数字字符&#xff0c;输出其前一个和后一个的数字字符&#xff0c;如果输入的是0前一个输出 “first”&#xff0c;9后一个则输出“last”&…

springBoot整合task

springBoot整合task 文章目录 springBoot整合task开开关设置任务&#xff0c;并设置执行周期定时任务的相关配置 开开关 设置任务&#xff0c;并设置执行周期 Component public class MyBean {Scheduled(cron "0/1 * * * * ?")public void print(){System.out.prin…

C++进阶篇6---C++11新语法

目录 目录 一、统一的列表初始化 二、声明 1.auto 2.decltype 3.nullptr 三、范围for 四、STL中的变化 五、右值引用和移动语义(重点) 一、统一的列表初始化 在c11之前&#xff0c;我们能用{}初始化数组和结构体 struct Point {int x;int y; }; int main() {int a[] …

机器学习 - 导论

简单了解 机器学习关于数据集的概念 、

堆排序算法

目录 1.基本原理2.例子3.代码实现 本文主要介绍堆排序的原理、例子以及代码实现。 1.基本原理 堆排序&#xff08;Heap Sort&#xff09;是一种基于比较的排序算法&#xff0c;它的工作原理是首先将待排序的序列构造成一个大顶堆或小顶堆&#xff0c;然后交换堆顶元素和最后一…

HCIP —— 双点重发布 + 路由策略 实验

目录 实验拓扑&#xff1a; 实验要求&#xff1a; 实验配置&#xff1a; 1.配置IP地址 2.配置动态路由协议 —— RIP 、 OSPF R1 RIP R4 OSPF R2 配置RIP、OSPF 双向重发布 R3配置RIP、OSPF 双向重发布 3.查询路由表学习情况 4.使用路由策略控制选路 R2 R3 5.检…

C/C++ 快速排序

个人主页&#xff1a;仍有未知等待探索_C语言疑难,数据结构,小项目-CSDN博客 专题分栏&#xff1a;算法_仍有未知等待探索的博客-CSDN博客 快速排序的思想——分治 目录 一、引言 二、讲解 1、步骤 2、代码 1.以左边界作为基准 2.以右边界作为基准 3.以中心点作为基准 …

Linux shell编程学习笔记32:declare 命令

0 前言 在 Linux shell编程学习笔记16&#xff1a;bash中的关联数组https://blog.csdn.net/Purpleendurer/article/details/134053506?spm1001.2014.3001.5501 中&#xff0c;我们在定义关联数组时使用了declare命令。 其实&#xff0c;declare命令的功能不只是定义定义关…

排序算法介绍(三)选择排序

0. 简介 选择排序&#xff08;Selection Sort&#xff09;是一种简单直观的排序算法。它的工作原理是每一次从待排序的数据元素中选出最小&#xff08;或最大&#xff09;的一个元素&#xff0c;存放在序列的起始位置&#xff0c;直到全部待排序的数据元素排完。选择排序是不稳…

超大规模集成电路设计----学习框架(一)

本文仅供学习&#xff0c;不作任何商业用途&#xff0c;严禁转载。绝大部分资料来自----数字集成电路——电路、系统与设计(第二版)及中国科学院段成华教授PPT 超大规模集成电路设计----学习框架&#xff08;一&#xff09; 这门课在学什么&#xff1f;这门课该怎么学&#xf…

PostgreSQL中常用的几种连接池总结及更新

前言 PostgreSQL的多进程结构&#xff0c;使得在支持大规模连接的时候&#xff0c;服务器端显得比较吃亏。一般上了1000个连接以上的时候&#xff0c;系统就会受到很大影响。这个时候&#xff0c;使用连接池&#xff0c;优势就会突显出来了。 在云环境下&#xff0c;一个JAVA…

Redis quicklist源码+listpack源码(5.0版本以上的优化)

ziplist设计上的问题&#xff0c;每一次增删改都需要计算前面元素的空间和长度&#xff08;prevlen&#xff09;&#xff0c;这种设计缺陷非常明显&#xff0c;因此引入了quicklist的设计。 quicklist quicklist实际就是双端链表&#xff0c;链表里的每一个节点都是ziplist&a…

Python---函数递归---练习:猴子吃桃问题(本文以递归算法 解法为主)

相关链接&#xff1a;Python---函数递归---练习&#xff1a;斐波那契数列&#xff08;本文以递归算法为主&#xff09;-CSDN博客 案例&#xff1a;猴子吃桃问题 猴子吃桃问题。猴子第1天摘下若干个桃子&#xff0c;当即吃了一半&#xff0c;还不过瘾&#xff0c;又多吃了一个。…

类和对象——(5)定义对象数组

归纳编程学习的感悟&#xff0c; 记录奋斗路上的点滴&#xff0c; 希望能帮到一样刻苦的你&#xff01; 如有不足欢迎指正&#xff01; 共同学习交流&#xff01; &#x1f30e;欢迎各位→点赞 &#x1f44d; 收藏⭐ 留言​&#x1f4dd; 芳华没有草稿纸&#xff0c;我们永久不…