共轭梯度法 Conjugate Gradient Method (线性及非线性)

1. 线性共轭梯度法

共轭梯度法(英语:Conjugate gradient method),是求解系数矩阵为对称正定矩阵的线性方程组的数值解的方法。

共轭梯度法是一个迭代方法,它适用于

1. 求解线性方程组,

2. 共轭梯度法也可以用于求解无约束的最优化问题

我们想最小化目标函数f(x),假设其拥有二次形式:

最优化问题表示如下:

上式可以等价于求解线性方程组 Ax = b,因为

目标方程的梯度为:

此外,双共轭梯度法(英语:BiConjugate gradient method)提供了一种处理非对称矩阵情况的推广。

共轭梯度法中,搜索方向p,是关于A共轭的,即

因此也称为共轭方向(conjugate directions)。

示例:

代码:

import numpy as np# Define the objective function
def f(x): return x[0]**2/2 + x[0]*x[1] + x[1]**2 - 2*x[1]# Define A and b
A = np.array(([1/2, 1/2], [1/2, 1]), dtype=float)
b = np.array([0., 2.])# Make sure A is a symmetric positive definite matrix.
if (A.T==A).all()==True: print("A is symmetric")
eigs = np.linalg.eigvals(A)
print("The eigenvalues of A:", eigs)
if (np.all(eigs>0)):print("A is positive definite")
elif (np.all(eigs>=0)):print("A is positive semi-definite")
else:print("A is negative definite")# Implements the linear conjugate gradient algorithm
def linear_CG(x, A, b, epsilon):res = A.dot(x) - b # Initialize the residualdelta = -res # Initialize the descent directionwhile True:if np.linalg.norm(res) <= epsilon:return x, f(x) # Return the minimizer x* and the function value f(x*)D = A.dot(delta)beta = -(res.dot(delta))/(delta.dot(D)) # Line (11) in the algorithmx = x + beta*delta # Generate the new iterateres = A.dot(x) - b # generate the new residualchi = res.dot(D)/(delta.dot(D)) # Line (14) in the algorithm delta = chi*delta -  res # Generate the new descent direction# Solve the equations
sol, funValue = linear_CG(np.array([2.3, -2.2]), A, b, 1e-5)# Check the result
if ( np.linalg.norm(A@sol - b) < 1e-5):print("solution verified")

2. 非线性共轭梯度法 Nonlinear Conjugate Gradient method

NCG方法被用来解非线性优化问题,常见下面几种算法:

  • Fletcher-Reeves algorithm,
  • Polak-Ribiere algorithm,
  • Hestenes-Stiefel algorithm,
  • Dai-Yuan algorithm, and
  • Hager-Zhang algorithm.

CG方法第一次被用来解非线性优化问题,是由Fletcher和Reeves提出的。搜索方向delta关于A共轭。

 chi计算式如下:

NCG的迭代公式为:

其中delta为方向,beta为步长。

示例:

代码:

首先安装自动微分工具 

pip install autograd

使用了autograd里面的梯度求解函数

此外,用到了 scipy中的line_search函数

scipy.optimize.line_search — SciPy v1.13.0 Manual

import numpy as np
from autograd import grad# Define the objective function
def func(x): # Objective functionreturn x[0]**4 - 2*x[0]**2*x[1] + x[0]**2 + x[1]**2 - 2*x[0] + 1Df = grad(func) # Gradient of the objective function# Next we define the function Fletcher_Reeves()
from scipy.optimize import line_search
NORM = np.linalg.normdef Fletcher_Reeves(Xj, tol, alpha_1, alpha_2):x1 = [Xj[0]]x2 = [Xj[1]]D = Df(Xj)delta = -D # Initialize the descent directionwhile True:start_point = Xj # Start point for step length selection beta = line_search(f=func, myfprime=Df, xk=start_point, pk=delta, c1=alpha_1, c2=alpha_2)[0] # Selecting the step lengthif beta!=None:X = Xj+ beta*delta #Newly updated experimental pointif NORM(Df(X)) < tol:x1 += [X[0], ]x2 += [X[1], ]return X, func(X) # Return the resultselse:Xj = Xd = D # Gradient at the preceding experimental pointD = Df(Xj) # Gradient at the current experimental pointchi = NORM(D)**2/NORM(d)**2 # Line (16) of the Fletcher-Reeves algorithmdelta = -D + chi*delta # Newly updated descent directionx1 += [Xj[0], ]x2 += [Xj[1], ]sol, funValue = Fletcher_Reeves(np.array([2., -1.8]), 10**-5, 10**-4, 0.38)
print(sol)

最后的解为 [1, 1],f(x)最小值为0

迭代过程如下:

## +----+-----------+------------+--------------+--------------+
## |    |       x_1 |        x_2 |         f(X) |     ||grad|| |
## |----+-----------+------------+--------------+--------------|
## |  0 |  2        | -1.8       | 34.64        | 49.7707      |
## |  1 | -0.98032  | -1.08571   |  8.1108      | 12.6662      |
## |  2 |  1.08966  |  0.0472277 |  1.30794     |  5.6311      |
## |  3 |  0.642619 |  0.473047  |  0.131332    |  0.877485    |
## |  4 |  0.766371 |  0.46651   |  0.0691785   |  0.260336    |
## |  5 |  0.932517 |  0.704482  |  0.0318138   |  0.583346    |
## |  6 |  1.0149   |  1.06008   |  0.00112543  |  0.110081    |
## |  7 |  1.02357  |  1.0596    |  0.000697231 |  0.0238509   |
## |  8 |  1.02489  |  1.05473   |  0.000638128 |  0.0331525   |
## |  9 |  1.00544  |  0.999549  |  0.000158528 |  0.0609372   |
## | 10 |  0.996075 |  0.987011  |  4.19723e-05 |  0.016347    |
## | 11 |  0.994792 |  0.986923  |  3.43476e-05 |  0.00538401  |
## | 12 |  0.994466 |  0.987575  |  3.25511e-05 |  0.00620548  |
## | 13 |  0.9956   |  0.992867  |  2.20695e-05 |  0.015708    |
## | 14 |  0.999909 |  1.00171   |  3.59093e-06 |  0.008628    |
## | 15 |  1.00088  |  1.00254   |  1.3779e-06  |  0.00206337  |
## | 16 |  1.00102  |  1.00249   |  1.24228e-06 |  0.000925229 |
## | 17 |  1.00106  |  1.00226   |  1.14704e-06 |  0.00161353  |
## | 18 |  1.00056  |  1.00065   |  5.3011e-07  |  0.00313135  |
## | 19 |  0.999916 |  0.99956   |  8.14653e-08 |  0.00107299  |
## | 20 |  0.999816 |  0.999511  |  4.85294e-08 |  0.000269684 |
## | 21 |  0.999798 |  0.999526  |  4.57054e-08 |  0.000185146 |
## | 22 |  0.999803 |  0.999615  |  3.90603e-08 |  0.000435884 |
## | 23 |  0.99995  |  0.999991  |  1.08357e-08 |  0.000499645 |
## | 24 |  1.00003  |  1.00009   |  2.25348e-09 |  0.000130632 |
## | 25 |  1.00004  |  1.00009   |  1.75917e-09 |  3.97529e-05 |
## | 26 |  1.00004  |  1.00009   |  1.66947e-09 |  4.22905e-05 |
## | 27 |  1.00003  |  1.00006   |  1.1931e-09  |  0.000108964 |
## | 28 |  1        |  0.999989  |  2.11734e-10 |  6.79786e-05 |
## | 29 |  0.999994 |  0.999982  |  7.24881e-11 |  1.61034e-05 |
## | 30 |  0.999993 |  0.999982  |  6.4458e-11  |  6.72611e-06 |
## +----+-----------+------------+--------------+--------------+

参考链接:

Chapter 5 Conjugate Gradient Methods | Introduction to Mathematical Optimization

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

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

相关文章

学习基于pytorch的VGG图像分类 day5

注&#xff1a;本系列博客在于汇总CSDN的精华帖&#xff0c;类似自用笔记&#xff0c;不做学习交流&#xff0c;方便以后的复习回顾&#xff0c;博文中的引用都注明出处&#xff0c;并点赞收藏原博主. 目录 VGG的数据集处理 1.数据的分类 2.对数据集的处理 VGG的分类标签设置 …

2款Notepad++平替工具(实用、跨平台的文本编辑器)

前言 今天大姚给大家分享2款Notepad平替工具&#xff0c;实用、跨平台&#xff08;支持Window/MacOS/Linux操作系统平台&#xff09;的文本编辑器。 NotepadNext NotepadNext是一个跨平台的 Notepad 的重新实现。开发是使用 QtCreator 和 Microsft Visual C (msvc) 编译器完…

python输入某年某月某日判断这一天是这一年的第几天

如何使用python实现输入某年某月某日判断这一天是这一年的第几天 from datetime import datetime #引入日期类 def is_leap_year(year):"""判断是否为闰年"""return (year % 4 0 and year % 100 ! 0) or (year % 400 0)# 根据年份和月份返回当…

深度解析 Spark(进阶):架构、集群运行机理与核心组件详解

关联阅读博客文章&#xff1a;深度解析SPARK的基本概念 引言&#xff1a; Apache Spark作为一种快速、通用、可扩展的大数据处理引擎&#xff0c;在大数据领域中备受关注和应用。本文将深入探讨Spark的集群运行原理、核心组件、工作原理以及分布式计算模型&#xff0c;带领读者…

test4141

欢迎关注博主 Mindtechnist 或加入【Linux C/C/Python社区】一起学习和分享Linux、C、C、Python、Matlab&#xff0c;机器人运动控制、多机器人协作&#xff0c;智能优化算法&#xff0c;滤波估计、多传感器信息融合&#xff0c;机器学习&#xff0c;人工智能等相关领域的知识和…

ES6: set和map数据结构以及使用场景

ES6:set和map数据结构 一、Set 数据结构&#xff1a;二、使用场景&#xff1a;使用Set 进行去重三、Map 数据结构四、使用场景&#xff1a;使用Map进行树型数据懒加载刷新五、Set和Map的区别六、Map、Set的实际使用场景 Set 和 Map 是 ES6 中引入的两种新的数据结构&#xff0c…

JavaScript中的Blob、Buffer、ArrayBuffer和TypedArray详解

文章的更新路线&#xff1a;JavaScript基础知识-Vue2基础知识-Vue3基础知识-TypeScript基础知识-网络基础知识-浏览器基础知识-项目优化知识-项目实战经验-前端温习题&#xff08;HTML基础知识和CSS基础知识已经更新完毕&#xff09; 正文 摘要&#xff1a;本文详细介绍了JavaS…

(UDP)其他信息: 通常每个套接字地址(协议/网络地址/端口)只允许使用一次。

“System.Net.Sockets.SocketException”类型的异常在 mscorlib.dll 中发生&#xff0c;但未在用户代码中进行处理其他信息: 通常每个套接字地址(协议/网络地址/端口)只允许使用一次。这个异常表示端口已经被占用了&#xff0c;需要释放端口或者使用其他端口来建立连接。您可以…

vite+react+ts+scss 创建项目

npm create vitelatest输入项目名称选择react选择typescript swc WC 通过利用 Rust 编写的编译器&#xff0c;使用了更先进的优化技术&#xff0c;使得它在处理 TypeScript 代码时能够更快地进行转换和编译。特别是在大型项目中&#xff0c;SWC 相对于传统的 TypeScript 编译器…

springboot火车运输物资定额领用管理系统java

本系统设计的现状和趋势&#xff0c;从需求、结构、数据库等方面的设计到系统的实现&#xff0c;分别为管理员和库房部、财务部、用户、生产部的实现。论文的内容从系统的设计、描述、实现、分析、测试方面来表明开发的过程。本系统根据现实情况来选择一种可行的开发方案&#…

2023数据要素白皮书(免费下载)

【1】关注本公众号&#xff0c;转发当前文章到微信朋友圈 【2】私信发送 【2023年数据资源入表白皮书】 【3】获取本方案PDF下载链接&#xff0c;直接下载即可。 如需下载本方案PPT原格式&#xff0c;请加入微信扫描以下方案驿站知识星球&#xff0c;获取上万份PPT解决方案&a…

STM32-看门狗

1、看门狗是什么&#xff1a;就是一个向下定时器&#xff0c;定时时间一到&#xff0c;就会触发一个向下的复位的中断&#xff0c;使单片机开始工作 2、作用&#xff1a;MCU微控制器构成的微型计算机系统中&#xff0c;由于微控制器的工作常常会受到来自外界电磁场的干 扰,造成…

MySQL进阶-合

目录 1.使用环境 2.条件判断 2.1.case when 2.2.if 3.窗口函数 3.1.排序函数 3.2.聚合函数 3.3.partiton by ​​​​​​​3.4.order by 4.排序窗口函数 5.聚合窗口函数 1.使用环境 数据库&#xff1a;MySQL 8.0.30 客户端&#xff1a;Navicat 15.0.12 MySQL进阶…

openjudge_2.5基本算法之搜索_166:The Castle

题目 166:The Castle 总时间限制: 1000ms 内存限制: 65536kB 描述 Figure 1 shows the map of a castle.Write a program that calculates how many rooms the castle hashow big the largest room is The castle is divided into m * n (m<50, n<50) square modules.…

React 组件生命周期对比:Class vs. 函数式

在 React 中&#xff0c;Class 组件和函数式组件的生命周期存在一些差异。通过对 React 中 Class 组件和函数式组件的生命周期进行对比&#xff0c;详细探讨了它们在设计哲学、生命周期管理和开发技巧上的异同。全面了解 React 中两种组件类型的生命周期特点&#xff0c;以及如…

MySQL-触发器:触发器概述、触发器的创建、查看删除触发器、 触发器的优缺点

触发器 触发器1. 触发器概述2. 触发器的创建2.1 创建触发器语法2.2 代码举例 3. 查看、删除触发器3.1 查看触发器3.2 删除触发器 4. 触发器的优缺点4.1 优点4.2 缺点4.3 注意点 注&#xff1a;此为笔者学习尚硅谷-宋红康MySQL的笔记&#xff0c;其中包含个人的笔记和理解&#…

Matlab与ROS(1/2)---Simulink(二)

0. 简介 在上一章中我们详细介绍了ROS与Matlab链接的基础用法。这一章我们将来学习如何使用Matlab当中的Simulink来完成。Simulink对机器人操作系统(ROS)的支持使我们能够创建与ROS网络一起工作的Simulink模型。ROS是一个通信层&#xff0c;允许机器人系统的不同组件以消息的形…

前端HTML入门基础5(表单)

前端HTML入门基础5&#xff08;表单&#xff09; 表单-基本结构文本框和密码框单选框和复选框隐藏域确认按钮重置按钮普通按钮文本域下拉框禁用表单控件label标签fieldset与legend的使用 表单-基本结构 在HTML中&#xff0c;表单&#xff08;Forms&#xff09;是一种允许用户与…

嵌入式实时操作系统的调度机制与优化

大家好&#xff0c;今天给大家介绍嵌入式实时操作系统的调度机制与优化&#xff0c;文章末尾附有分享大家一个资料包&#xff0c;差不多150多G。里面学习内容、面经、项目都比较新也比较全&#xff01;可进群免费领取。 嵌入式实时操作系统的调度机制与优化 一、引言 嵌入式实…

最前沿・量子退火建模方法(1) : subQUBO讲解和python实现

前言 量子退火机在小规模问题上的效果得到了有效验证&#xff0c;但是由于物理量子比特的大规模制备以及噪声的影响&#xff0c;还没有办法再大规模的场景下应用。 这时候就需要我们思考&#xff0c;如何通过软件的方法怎么样把大的问题分解成小的问题&#xff0c;以便通过现在…