文章目录
- 写在前面
- 共轭梯度法
- 代码实现
- 标题函数定义部分
- 函数测试
- 画图显示效果
- 显示方法
- 运行方法
写在前面
写这篇博客是为了增加对共轭梯度的理解。最近最优化课一直在讲共轭梯度,雅克比,梯度下降啊,课上没听太懂,就课下花了点时间好好学一下,在此记录一下。
这篇博客主要介绍了共轭梯度的matlab实现,代码可以运行看到效果。代码分为2个部分,第一部分代码是函数定义,第二部分代码是测试代码。
共轭梯度法
我是参考了下面这篇博客,以及B站的一个视频
共轭梯度法的简单直观理解https://www.bilibili.com/video/BV16a4y1t76z/
如果你了解梯度下降法的话,理解共轭梯度就很容易了。概念性的东西和数学公式我就不说了,用大白话直接说说吧,有需要数学公式推导的可以去检索一下。
首先是初始点,你在初始点想去下一个点该怎么走,可以把它想象成一个人下山,第一步是不是先确定方向,第二步是不是确定走多远。共轭梯度法就是一直迭代上面2个步骤。
确定方向该怎么确定,最速下降法是直接是梯度的方向,共轭梯度法在初始点时是梯度的方向,接着下一个点就不是梯度的方向了,它要一直更新了。
r0=b-A*x0 %起始点的梯度方向
p0=r0 %方向
新的方向用p1表示
p1=r1-bat*p0
可以把bat就是一个系数
r1= b-A*x1
r1是新的残差 ,b是已知的,
x1是下一个点的位置
x1=x0+alp*p0
x0是初始点的位置
现在可以连起来,假设你站在x0的位置上方向p0,走多远是alp,到达下一个点x1。
到达下一个点之后从新计算方向,新的方向就是上面的p1,然后已知循环迭代就行了。
上方的bat和alp系数我就不介绍了可以直接按照公式来。
代码实现
标题函数定义部分
function[x,obj]=cg_solver(A,b,x0,epsilon,max_iter)
x=x0; %初始点
i=0;
r=b-A*x; % 计算残差
d=r; %方向
delta_new=r'*r; %残差的内积
delta0=delta_new;
obj=0.5*x'*A*x-b'*x; %目标函数
while(i<max_iter)&&(delta_new/delta0>epsilon^2)
q=A*d;
alpha=delta_new/(d'*q); %计算系数 alpha
x=x+alpha*d; %迭代一步(向下走一步)
obj=[obj,0.5*x'*A*x-b'*x];
r=b-A*x; % 更新计算残差
delta_old=delta_new;
delta_new=r'*r; %更新残差的内积
beta=delta_new/delta_old; %更新下一步方向的系数
d=r+beta*d; %更新下一步的方向
i=i+1;
end
end
函数测试
% Ax=b
x_true=ones(1000,1);%x_true是向量
A=diag(1:1000); %矩阵对角元素的提取和创建对角阵
b=A*x_true;
x0=zeros(1000,1); %x0是向量
[x,obj]=cg_solver(A,b,x0,1e-8,1000);
画图显示效果
显示方法
步骤1
右键点击obj,选择plot(obj),就会出现上面的效果图
运行方法
一个是函数定义部分也就是你实现共轭梯度法,下图中的文件2,文件1是测试文件。点击文件1运行,然后最右边回显示工作区,按显示方法进行操作显示。