一、算法原理
1、引入
先来回顾一下之前我们介绍的最速下降法;
假设函数为f(x),最速下降法通过给定一个初始点xk,选择xk处的负梯度方向为最速下降方向,然后进行线搜索来确定步长。
其迭代公式为:x(k+1)=xk+a*gradient(f)/norm(gradient(f)),
其中a为步长,∇f(xk)=gradient(f)/norm(gradient(f))为梯度。
步长a要满足条件:f(x(k+1))=min(f(x(k+1))。
2、最速下降法每次都是在局部最优的方向上选取了最优的步长,但是在其后面的迭代过程中可能会出现相同的方向,这样使
我们在迭代过程中多次对一个方向进行迭代,在极值点附近的迭代轨迹为锯齿波,影响了收敛速度。
共轭梯度法解决了这个问题,每一次迭代都将该方向优化为最优方向,理论上对n维问题,只需求出n个方向上的极小值。
3、对于多元二次函数,例如f(x)=x1^2+x2^2-4x1-4x2+8=(x1-2)^2+(x2-2)^2,我们可以将其写成如下形式,
=[x1 [2 [x1 * * [x1 x2]*1/2 + * [-4 -4] +8 |
其中A=[2;2],b=[-4;-4],c=8。为求其解,我们需要假定一个初始点x0,但是我们怎样进行迭代呢?
二、算法步骤
1、初始优化方向为f(x)的负梯度方向,p0=r0=-∇f(x0)=b-Ax0,将该方向记为r0。
2、按照最速下降法进行一次迭代;
3、(1)下一次优化方向pk的选取
每一次的优化方向要与之前的优化方向正交,这里采时线形代数学过的smith正交化就派上用场了;
(2)下一次优化步长ak的选取
f(x(k+1))=f(xk+ak*pk),对其求导,导数为零的点的ak的值即为最优步长(与最速下降法相同),可得。
由于博主水平有限,只是大致讲解一下共轭梯度法的计算原理,随着博主水平的加深后面会不断丰富其内容。
三、matlab代码
%% 共轭梯度
clc
clear
f=@(x1,x2) x1.^2+x2.^2-x1*x2-10*x1-4*x2+60;
[x,result]=Min_gongetidu(f,[0;0])
function [x,result]=Min_gongetidu(f,x0,tol)
if nargin == 2 %输入参数为2,默认误差为le-6
tol=1e-6;
end
%求当前点梯度
%% FF为了计算梯度值 ff为了计算函数值,f_x1 是 带入x1后的表达式
F_td=matlabFunction(gradient(sym(f)));%求梯度并获得梯度函数句柄
x01=x0(1);
x02=x0(2);
f_td_0=F_td(x01,x02);%初始优化方向为负梯度方向
d_0=-f_td_0;%初始优化方向为负梯度方向
if norm(f_td_0) < tol %若初始点为最小点,输出结果
x=x0;
x01=x0(1);
x02=x0(2);
result=f(x01,x02);
return;
end
%%步长alfa
syms alfa %按照最速下降法进行第一次迭代
x1=x0(:)+alfa*d_0;%将这个点带入函数,求alfa
x11=x1(1);
x12=x1(2);
fx1=f(x11,x12);%计算x1带入后的原函数表达式
d_x1=diff(fx1);%对这个表达是求导
d_alfa=double(solve(d_x1));%求解表达式为0时的alfa
x0=x0(:)+d_alfa*d_0;
x01=x0(1);
x02=x0(2);
f_td_1=F_td(x01,x02);%计算当前梯度
if norm(f_td_1) < tol %如果一次迭代达到精度要求,输出结果
x=x0;
x01=x0(1);
x02=x0(2);
result=f(x01,x02);
return;
end
while 1 %共轭梯度法大循环
beta0=norm(f_td_1)^2/norm(f_td_0)^2; %优化方向的确定
d_1=-f_td_1+beta0*d_0; %优化方向的确定
x1=x0(:)+alfa*d_1;%优化步长的确定
x11=x1(1);
x12=x1(2);
fx1=f(x11,x12);%计算x1带入后的原函数表达式
d_x1=diff(fx1);%对这个表达是求导
d_alfa=double(solve(d_x1));%求解表达式为0时的alfa
x0=x0(:)+d_alfa*d_1;%下一个迭代点
x01=x0(1);
x02=x0(2);
d_0=d_1;
f_td_0=f_td_1;
f_td_1=F_td(x01,x02); %计算当前梯度
if norm(f_td_1) < tol %判断
x=x0;
x01=x0(1);
x02=x0(2);
result=f(x01,x02);
return;
end
end
end