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

2024-04-14 08:44

本文主要是介绍共轭梯度法 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

这篇关于共轭梯度法 Conjugate Gradient Method (线性及非线性)的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/902590

相关文章

BD错误集锦7——在集成Spring MVC + MyBtis时使用c3p0作为数据库时报错Method com/mchange/v2/c3p0/impl/NewProxyPreparedStatem

异常信息如下: Type Exception ReportMessage Handler dispatch failed; nested exception is java.lang.AbstractMethodError: Method com/mchange/v2/c3p0/impl/NewProxyPreparedStatement.isClosed()Z is abstractDescr

线性回归(Linear Regression)原理详解及Python代码示例

一、线性回归原理详解         线性回归是一种基本的统计方法,用于预测因变量(目标变量)与一个或多个自变量(特征变量)之间的线性关系。线性回归模型通过拟合一条直线(在多变量情况下是一条超平面)来最小化预测值与真实值之间的误差。 1. 线性回归模型         对于单变量线性回归,模型的表达式为:         其中: y是目标变量。x是特征变量。β0是截距项(偏置)。β1

DataWhale机器学习——第三章线性模型笔记

读书笔记: 《机器学习》第三章 线性模型 3.1 基本形式 3.1.1 线性模型的定义 3.1.2 线性模型的优点 简单易理解计算效率高容易实现和解释 3.1.3 线性模型的局限性 只能表达线性关系对于复杂的非线性关系,表现较差 3.2 线性回归 3.2.1 基本概念 3.2.2 最小二乘法 正规方程 3.2.3 正则化 为防止过拟合,可以在损失函数中加入正

让IE8支持CSS3属性(border-radius、box-shadow、linear-gradient)

下载 PIE-1.0.0.zip解压后,将文件夹重命名为PIE,放到项目目录下在CSS3文件中添加一行代码 behavior: url(PIE/PIE.htc); 例如: .form__input{border-radius: 0.3em;behavior: url(PIE/PIE.htc);} 参考: TYStudio-专注WEB前端开发 css3pie

智能优化算法改进策略之局部搜索算子(六)--进化梯度搜索

1、原理介绍     进化梯度搜索(Evolutionary Gradient Search, EGS)[1]是兼顾进化计算与梯度搜索的一种混合算法,具有较强的局部搜索能力。在每次迭代过程中,EGS方法首先用受进化启发的形式估计梯度方向,然后以最陡下降的方式执行实际的迭代步骤,其中还包括步长的自适应,这一过程的总体方案如下图所示:     文献[1]

使用matlab的大坑,复数向量转置!!!!!变量区“转置变量“功能(共轭转置)、矩阵转置(默认也是共轭转置)、点转置

近期用verilog去做FFT相关的项目,需要用到matlab进行仿真然后和verilog出来的结果来做对比,然后计算误差。近期使用matlab犯了一个错误,极大的拖慢了项目进展,给我人都整emo了,因为怎么做仿真结果都不对,还好整体的代码都是我来写的,慢慢往下找就找到了问题的来源,全网没有看到多少人把这个愚蠢的错误写出来,我来引入一下。 代码错误的表现:复数向量的虚部被取反,正数变成负数,负数

机器学习算法(二):1 逻辑回归的从零实现(普通实现+多项式特征实现非线性分类+正则化实现三个版本)

文章目录 前言一、普通实现1 数据集准备2 逻辑回归模型3 损失函数4 计算损失函数的梯度5 梯度下降算法6 训练模型 二、多项式特征实现非线性分类1 数据准备与多项式特征构造2 逻辑回归模型 三、逻辑回归 --- 正则化实现1 数据准备2 逻辑回归模型3 正则化损失函数4 计算损失函数的梯度5 梯度下降6 训练模型 总结 前言 今天我们开始介绍逻辑回归的从零开始实现代码了,

「Debug R」报错unable to find an inherited method for function是如何产生的

在一个群里看到这样一条报错,截图如下: 报错信息 当然这种问题解决起来也很快,无非就是把报错信息复制出来放在搜索引擎上,只不过你要挑选合适的搜索引擎。 百度 谷歌 必应 解决方案就是用dplyr::select。 虽然报错解决了,但是我还想着要重复出这个报错。因为只有能重复出报错,才能证明你不是运气好才解

智能优化算法改进策略之局部搜索算子(四)--梯度搜索法

2、仿真实验 以海洋捕食者算法(MPA)为基本算法。考察基于梯度搜索的改进海洋捕食者算法(命名为GBSMPA) vs. 海洋捕食者算法(MPA)  在Sphere函数上的比较      在Penalized1函数上的比较    在CEC2017-1上的比较    在CEC2017-3上的比较 在CEC2017-4上的比较 代码获取:

abstract的method是否可同时是static,是否可同时是native,是否可同时是synchronized

1,abstract的method是否可同时是static,是否可同时是native,是否可同时是synchronized 都不可以,因为abstract申明的方法是要求子类去实现的,abstract只是告诉你有这样一个接口,你要去实现,至于你的具体实现可以是native和synchronized,也可以不是,抽象方法是不关心这些事的,所以写这两个是没有意义的。然后,static方法是不会被覆