[最优化理论] 梯度下降法 + 精确线搜索（单峰区间搜索 + 黄金分割）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
#endif

using 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 search
        if (phi1 < phi0)
        {
            t      = 2.0 * t;
            alpha  = alpha0;
            alpha0 = alpha1;
            phi0   = phi1;
        }
        else
        {
            // backward search
            if (k == 0)
            {
                t     = -t;
                alpha = alpha1;
            }
            // find another end
            else
            {
                break;
            }
        }
        ++k;
    }

#if DEBUG
    std::cout << "ForwardBackward total iteration = " << std::endl;
    std::cout << k << std::endl;
#endif

    XClass 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 DEBUG
        std::cout << "GoldenSectionSearch total it eration = " << std::endl;
        std::cout << k << std::endl;
#endif

        return {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 DEBUG
    std::cout << "GoldenSectionSearch total iteration = " << std::endl;
    std::cout << k << std::endl;
#endif

    if (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 pointer
    std::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 file
    fout << x0[0] << ", " << x0[1] << ", " << phi_min << "\n";
#endif

    while (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 DEBUG
        std::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 file
        fout << 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.mfbz.cn/a/213589.html

如若内容造成侵权/违法违规/事实不符，请联系我们进行投诉反馈qq邮箱809451989@qq.com，一经查实，立即删除！

相关文章

【Tkinter 入门教程】

【Tkinter 入门教程】

【Tkinter 入门教程】 1. Tkinter库的简介&#xff1a;1.1 GUI编程1.2 Tkinter的定位 2. Hello word! 程序起飞2.1 第⼀个程序2.2 字体颜色主题 3. 组件讲解3.1 tkinter 的核⼼组件3.2 组件的使⽤3.3 标签Label3.3.1 标签显示内容3.3.2 多标签的应⽤程序3.3.3 总结 3.4 按钮but…
阅读更多...
leetcode 201 数字范围按位与

leetcode 201 数字范围按位与

leetcode 201 题目题解代码 题目 给你两个整数 left 和 right &#xff0c;表示区间 [left, right] &#xff0c;返回此区间内所有数字 按位与 的结果&#xff08;包含 left、right 端点&#xff09;。 具体示例如下&#xff1a; 题解 本题是一个在思维上的方法&#xff0c;不…
阅读更多...
Ant Design Pro 框架设置API Token拦截器的功能

Ant Design Pro 框架设置API Token拦截器的功能

分享记录一个解决方法&#xff0c;希望对大家有帮助。 找到文件&#xff0c;然后定义一个方法。最后调用一下即可。 代码我也给你贴上了。 // 获取token 拦截方法 const setTokenRequest (config: any) > {const token 30|eh5GNXWRe5rO4XLjbbnqy132RABfiKqI338EoIhqc790a…
阅读更多...
sourceTree的下载和安装

sourceTree的下载和安装

sourceTree的下载和安装 一、概述 SourceTree 是一款免费的 Git 和 Hg 客户端管理工具&#xff0c;支持 Git 项目的创建、克隆、提交、push、pull 和合并等操作。它拥有一个精美简洁的界面&#xff0c;大大简化了开发者与代码库之间的 Git 操作方式&#xff0c;这对于不熟悉 …
阅读更多...
java学习part32StringBuffer和StringBuilder

java学习part32StringBuffer和StringBuilder

Java中的值传递和引用传递&#xff08;详解&#xff09; - 知乎 (zhihu.com) 146-常用类与基础API-StringBuffer与StringBuilder的源码分析、常用方法_哔哩哔哩_bilibili 1. 2.扩容机制 不够用&#xff1a;长度为 原长度*22&#xff1b;如果还不够&#xff0c;那么就扩容到目…
阅读更多...
STM32踩坑--串口发送乱码

STM32踩坑--串口发送乱码

一、发现问题 今天在STM32F407新板子上测试串口时&#xff0c;发现发送数据一直乱码。 二、解决问题 针对STM32F407系列校准PLLCLK时钟&#xff1a; ①由 时钟树 可以看出PLLCLKHSE&#xff08;高速外部时钟&#xff09;*N/(M*P)。因为SYSTICK一般取最高的时钟168M&#xff…
阅读更多...
【Linux】第二十五站：深入理解文件系统

【Linux】第二十五站：深入理解文件系统

文章目录 一、前言二、认识硬件----磁盘1.基本介绍2.磁盘的存储构成3.磁盘的逻辑结构4.回归到硬件 三、文件系统1.划分2.Block group(1)Data blocks(2)inode Table(3)Block Bitmap(4)inode Bitmap(5)Group Descriptor Table(GDT)(6)Super Block 3.总结4.一些其他问题5.如何理解…
阅读更多...
C语言--求一个十进制整数中1的个数

C语言--求一个十进制整数中1的个数

一.题目描述⭐ 求一个十进制整数中1的个数 比如&#xff1a; 输入:10201 输出&#xff1a;2 &#xff08;这个数字中1的个数是2&#xff09; 二.思路分析⭐ 数字类的问题我们可以用取模&#xff0c;或者取余运算。 首先定义一个计数器&#xff0c;用来统计1的个数。 输入数字…
阅读更多...
《管家婆》辉煌2005+（V4.0）简单教程

《管家婆》辉煌2005+（V4.0）简单教程

《管家婆》辉煌2005&#xff08;V4.0&#xff09;简单教程 呉師傅 运行环境&#xff1a; 　　操作系统推荐使用Win2000&#xff08;32位&#xff09;、WinXP&#xff08;32位&#xff09;、Win7&#xff08;位&#xff09; 兼容&#xff1a;Win7&#xff08;64位&#xff09…
阅读更多...
【毕业设计】基于雷达与深度学习的摔倒检测——微多普勒效应

【毕业设计】基于雷达与深度学习的摔倒检测——微多普勒效应

运动物体的微多普勒效应为人体动作识别提供了可能&#xff0c;基于雷达的居家检测具有良好的隐私保护性&#xff0c;且不易受环境因素影响&#xff08;如光照、温度等&#xff09;&#xff0c;近年来已受到国内外学者的广泛关注。由于雷达信号的非平稳特性&#xff0c;通过短时…
阅读更多...
html电子签名

html电子签名

html电子签名 html5实现手写签名板&#xff0c;点击保存即可生成base64格式的图片 使用H5自带的canvas&#xff0c;无需引入js无需引入别的js 效果图 html代码 <!DOCTYPE html> <html> <head><meta http-equiv"Content-Type" content"…
阅读更多...
最大乘积分解（动态规划）

最大乘积分解（动态规划）

相较于我上一题写的动态规划&#xff0c;这一题比较简单 代码如下&#xff1a; #include<stdio.h>int main(void) {long long n, max[101] {0, 1};scanf("%lld", &n);for(int i 1; i < n; i)max[i] i;for(int i 1; i < n; i)for(int j 1; j &…
阅读更多...
java蚁群算法的物流管理系统eclipse定制开发mysql数据库BS模式java编程百度地图

java蚁群算法的物流管理系统eclipse定制开发mysql数据库BS模式java编程百度地图

一、源码特点 java 基于蚁群算法的物流管理系统是一套完善的web设计系统 &#xff0c;对理解JSP java编程开发语言有帮助&#xff0c;系统具有完整的源代码和数据库&#xff0c;系统主要采用B/S模式开发。开发环境为 TOMCAT7.0,eclipse开发&#xff0c;数据库为Mysql5.0&a…
阅读更多...
[linux进程控制]进程替换

[linux进程控制]进程替换

文章目录 1.进程替换的概念和原理2.如何完成进程替换2.1exec系列函数加载器的底层系统调用接口基于execve的封装--库函数 2.2 int execl(const char *path, const char *arg, ...);1.在当前进程进行替换2.在子进程进行替换 2.3 int execv(const char *path, char *const argv[]…
阅读更多...
基于JNI实现调用C++ SDK

基于JNI实现调用C++ SDK

基于JNI实现调用C SDK 背景分析解决实践 背景 上篇文章总结了几种Java项目调用C/C SDK项目方法&#xff0c;在逐一实践、踩坑后&#xff0c;最终还是敲定采用 JNI 方式进行实现。在文章开始的过程&#xff0c;会先大概讲讲笔者遇到的情况&#xff0c;因为封装方式需要根据实际…
阅读更多...
Debian12配置ssh服务器

Debian12配置ssh服务器

Debian12配置ssh服务器 安装ssh-server sudo apt install openssh-server启动ssh sudo systemctl start ssh启用ssh sudo systemctl enable ssh查看ssh状态 sudo systemctl status ssh可以看到有enabled和running字样 说明ssh启用成功 连接到服务器 # username是你的用…
阅读更多...
React立即更新DOM

React立即更新DOM

正常情况下&#xff0c;react会等待set完毕后再进行页面渲染&#xff0c;所以在set时无法拿到更新后的dom import { useRef, useState } from "react"export default () > {const div useRef(null)const [count, setCount] useState(0)const btnClick () >…
阅读更多...
备战春招——12.3 算法

备战春招——12.3 算法

哈希表 哈希表主要是使用 map、unordered_map、set、unorerdered_set、multi_&#xff0c;完成映射操作&#xff0c;主要是相应的函数。map和set是有序的&#xff0c;使用的是树的形式&#xff0c;unordered_map和unordered_set使用的是散列比表的&#xff0c;无序。 相应函数…
阅读更多...
【PUSDN】java中easyexcel导入导出带有图片的Excel（main方法方式）

【PUSDN】java中easyexcel导入导出带有图片的Excel（main方法方式）

简述 java中easyexcel导入导出带有图片的Excel&#xff08;main方法方式&#xff09;&#xff0c;web方式详见另一篇 由于电脑音频问题&#xff0c;视频暂时没有解说声音&#xff0c; 回头重新补上 前情提示 如果有任何疑问、需求、技术支持&#xff0c;欢迎点赞&#xff0…
阅读更多...
微机原理——定时器8253(8254)学习2应用与设计

微机原理——定时器8253(8254)学习2应用与设计

目录 简要说明 用户扩展的定时计数器应用举例 1 8254作测量脉冲宽度 2 8254作定时 3 8254作分频 4 8254同时用作计数与定时 硬件设计 ​编辑软件设计 微机系统中定时计数器应用举例 5 计时器设计 硬件设计 软件设计 6 发生器设计 硬件设计 软件设计 简要说明 定…
阅读更多...
最新文章

Copyright @ 2022~2023

友情链接： 我的编程经验分享 一条长河网 尧图网站开发 尧图建网站