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

Posted 2021-11-23 Z.Q.Feng

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了线性共轭梯度法求解正定二次函数极小点以及线性方程组的解--MATLAB源程序相关的知识，希望对你有一定的参考价值。

实现原理

具体数学实现原理可参考这篇文章：最速下降法/steepest descent，牛顿法/newton，共轭方向法/conjugate direction，共轭梯度法/conjugate gradient 及其他

拟解决的问题

求解正定二次函数的极小点

对于正定二次函数
$\\frac{1}{2} x^T G x + b^T x + c$

我们有线性共轭梯度算法伪码如下：

Step1：给定初始点 $x_0$ 和容许误差 $\\epsilon > 0$ ，令 $k = 0$ ;
Step2：计算 $g_k = G x_k + b$ ，若 $||g_k|| < \\epsilon$ ，则迭代停止 ;
Step3：若 $k = 0$ ，令 $d_k = -g_k$ ，否则计算

$\\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}$ ;
Step4：计算步长 $\\alpha_k = \\dfrac{g_k^T g_k}{d_k^T G d_k}$ ;
Step5：令 $x_{k+1} = x_k + \\alpha_k d_k, k = k + 1$ ，go to Step2.

根据此算法可求解出 f(x) 的极小点，即令 f(x) 的梯度函数 g(x) 值为 0 的点值。

求解线性方程组的根

对于线性方程组 A x = b ，我们记 A 为方程组的系数矩阵，x 为方程组的根向量，b 为方程组的值向量，我们可以发现，对于正定二次函数
$\\frac{1}{2} x^T G x + b^T x + c$

它的梯度函数为
$g(x) = G x + b^T$

若我们令 g(x) = 0，则可得到如下表达式
$G x = - b^T$

因此，若我们做出如下构建，

令线性方程组的系数矩阵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

>>