具体数学实现原理可参考这篇文章:最速下降法/steepest descent,牛顿法/newton,共轭方向法/conjugate direction,共轭梯度法/conjugate gradient 及其他
对于正定二次函数
f
(
x
)
=
1
2
x
T
G
x
+
b
T
x
+
c
f(x) = \frac{1}{2} x^T G x + b^T x + c
f(x)=21xTGx+bTx+c
我们有线性共轭梯度算法伪码如下:
Step1:给定初始点
x
0
x_0
x0 和容许误差
ϵ
>
0
\epsilon > 0
ϵ>0 ,令
k
=
0
k = 0
k=0 ;
Step2:计算
g
k
=
G
x
k
+
b
g_k = G x_k + b
gk=Gxk+b ,若
∣
∣
g
k
∣
∣
<
ϵ
||g_k|| < \epsilon
∣∣gk∣∣<ϵ ,则迭代停止 ;
Step3:若
k
=
0
k = 0
k=0 ,令
d
k
=
−
g
k
d_k = -g_k
dk=−gk ,否则计算
β
k
−
1
=
g
k
T
g
k
g
k
−
1
T
g
k
−
1
,
d
k
=
−
g
k
+
β
k
−
1
d
k
−
1
\beta_{k-1} = \dfrac{g_k^T g_k}{g_{k-1}^T g_{k-1}}, d_k = -g_k + \beta_{k-1} d_{k-1}
βk−1=gk−1Tgk−1gkTgk,dk=−gk+βk−1dk−1 ;
Step4:计算步长
α
k
=
g
k
T
g
k
d
k
T
G
d
k
\alpha_k = \dfrac{g_k^T g_k}{d_k^T G d_k}
αk=dkTGdkgkTgk ;
Step5:令
x
k
+
1
=
x
k
+
α
k
d
k
,
k
=
k
+
1
x_{k+1} = x_k + \alpha_k d_k, k = k + 1
xk+1=xk+αkdk,k=k+1 ,go to Step2.
根据此算法可求解出 f(x) 的极小点,即令 f(x) 的梯度函数 g(x) 值为 0 的点值。
对于线性方程组 A x = b ,我们记 A 为方程组的系数矩阵,x 为方程组的根向量,b 为方程组的值向量,我们可以发现,对于正定二次函数
f
(
x
)
=
1
2
x
T
G
x
+
b
T
x
+
c
f(x) = \frac{1}{2} x^T G x + b^T x + c
f(x)=21xTGx+bTx+c
它的梯度函数为
g
(
x
)
=
G
x
+
b
T
g(x) = G x + b^T
g(x)=Gx+bT
若我们令 g(x) = 0,则可得到如下表达式
G
x
=
−
b
T
G x = - b^T
Gx=−bT
因此,若我们做出如下构建,
令线性方程组的 系数矩阵A 为正定二次函数的 Hessen矩阵G,以线性方程组的 值向量b 构建 -b 作为正定二次函数中一次项前面的系数,则正定二次函数的极小点即为线性方程组的解。
利用线性共轭梯度法求解二次函数极小值,该极小值即为构建线性方程组的解:
function [x, k] = Conjudate(G, b, x0, eps, kmax) % % function [x, k] = Conjudate(G, b, x0, kmax, eps) % 线性共轭梯度法求解正定二次函数 f(x) = 1/2 x' G x + b' x % 的极小点,即求解线性方程 G x = b 的过程,故最终输出 x 既为 % 正定二次函数的极小点,也为线性方程的根 % --------------------------------------------------- % 输入: % G n 阶正定对称矩阵,正定二次函数的Hessen矩阵, % 也为线性方程的系数矩阵 % b 正定二次函数中的一次项系数,也为线性方程的值向量乘以负一 % x0 初始点 % eps 精确度 % kmax 函数的最大迭代次数 % % 输出: % x 正定二次函数的极小点,线性方程的根 % k 算法迭代次数 % --------------------------------------------------- % by Zhi Qiangfeng % gk = G * x0 + b; % x0 处的正定二次函数梯度值 dk = -gk; % 初始下降方向 xk = x0; k = 0; while k <= kmax if norm(gk) < eps break end gk_ = gk; % 迭代前的梯度值 gk = G * xk + b; % 迭代后的梯度值 dk_ = dk; % 上一次选取的下降方向 if k == 0 dk = -gk; % 初始点的下降方向为负梯度方向 else beta = (gk' * gk) / (gk_' * gk_); dk = -gk + beta * dk_; % 共轭梯度法迭代方向 end % 正定二次函数步长公式 alpha = (gk' * gk) / (dk' * G * dk); xk = xk + alpha * dk; % 更新迭代点 k = k + 1; end x = xk; end
求解线性方程组 Ax = b,其中 A 为 n 阶希尔伯特矩阵,对于希尔伯特矩阵可参考百度百科:希尔伯特矩阵
b为系数矩阵A每行的和构成的列向量,其中我们初始阶数选取为40,即未知数个数为40。
根据线性代数知识,我们有
对于值向量b,我们有
对比 Ax 与 b,我们不难得出线性方程组 Ax = b 的解为
编写脚本文件 linear_equ.m,内容如下
clear; clc; A = hilb(40); % 我们选取希尔伯特矩阵作为系数矩阵 b = sum(A, 2); % b为nx1矩阵,为希尔伯特矩阵每行的和构成的向量 x0 = zeros(40, 1); % 初始点我们选为0 kmax = 1000; % 最大迭代次数 eps = 1e-6; % 精度 [x, k] = Conjudate(A, -b, x0, kmax, eps)
其中,A 和 b 分别为线性方程组的系数矩阵和值向量,x0为正定二次函数迭代的初始点,x为所构建正定二次函数的极小点,即线性方程组 Ax = b 的解,得到输出如下:
>> linear_equ x = 1.0001 0.9990 1.0035 0.9982 0.9971 0.9987 1.0006 1.0019 1.0023 1.0021 1.0015 1.0007 0.9998 0.9991 0.9985 0.9981 0.9980 0.9980 0.9982 0.9985 0.9989 0.9993 0.9998 1.0003 1.0008 1.0012 1.0016 1.0019 1.0021 1.0022 1.0021 1.0020 1.0017 1.0012 1.0007 1.0000 0.9992 0.9982 0.9971 0.9959 k = 10 >>
可以看到,与精确解十分接近!
最优化相关算法设计数学原理:最优化/Optimization文章合集
最优化相关MATLAB实现程序:Z.Q.Feng的最优化笔记
有帮助可以点赞哦,谢谢大家的支持~