ICode9

精准搜索请尝试: 精确搜索
首页 > 编程语言> 文章详细

线性共轭梯度法求解正定二次函数极小点以及线性方程组的解--MATLAB源程序

2021-06-11 21:01:25  阅读:279  来源: 互联网

标签:函数 求解 -- 线性方程组 矩阵 正定 二次 MATLAB 源程序


目录

实现原理

具体数学实现原理可参考这篇文章:最速下降法/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)=21​xTGx+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−1T​gk−1​gkT​gk​​,dk​=−gk​+βk−1​dk−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​=dkT​Gdk​gkT​gk​​ ;

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​+αk​dk​,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)=21​xTGx+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的最优化笔记


有帮助可以点赞哦,谢谢大家的支持~

标签:函数,求解,--,线性方程组,矩阵,正定,二次,MATLAB,源程序
来源: https://blog.csdn.net/weixin_46584887/article/details/117810760

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有