共轭梯度法是非常重要的一种优化算法，仅仅需要函数的一阶导数，步收敛性，稳定性高，而且不需要任何外来参数。共轭梯度法的本质是选择一个矩阵，使之尽量接近海森矩阵的逆。共轭梯度法经常用于求解多元函数的极值，今天主要是以实际的例子讲解MATLAB的共轭梯度法求多元函数的极值。共轭梯度法的计算步骤如下：
1.实例
主程序
clc;
clear all;
close all;
syms x y;%定义函数变量 x y
x0 = [1 1]';
[x,val,k]=frcg(@fun111,@gfun222,x0)
x = -5:0.01:5;
y = x;
[x,y] = meshgrid(x,y);
F = x.^6+2*x.^3+24*x.^2+y.^6+12*y.^2;
figure;
mesh(x,y,F);
xlabel('x');
ylabel('y');
zlabel('f');目标函数自定义函数
function f = fun111(X)
x = X(1);
y = X(2);
f = x^6+2*x^3+24*x^2+y^6+12*y^2;
end梯度函数自定义函数
function Y = gfun222(X)
x = X(1);
y = X(2);
Y = [6*x^5 + 6*x^2 + 48*x;
6*y^5 + 24*y];
end共轭梯度法函数
function [x,val,k]=frcg(fun,gfun,x0)
%功能：用FR共辄梯度法求解无约束问题：min f(x)
%输入：x0是初始点，fun,gfun分别是目标函数和梯度
%输出：x,val分别是近似最优点和最优值，k是迭代次数.
maxk=5000; %最大迭代次数
rho=0.6;
sigma=0.4;
k=0;
epsilon=1e-6;
n=length(x0);
while(k<maxk)
g=feval(gfun,x0);%计算梯度
itern=k-(n+1)*floor(k/(n+1));
itern=itern+1;
%计算搜索方向
if(itern==1)
d=-g;
else
beta=(g'*g)/(g0'*g0);
d=-g+beta*d0;gd=g'*d;
if(gd>=0.0)
d=-g;
end
end
if(norm(g)<epsilon), break; end %检验终止条件
m=0; mk=0;
while(m<20) %Armijo 搜索
if(feval(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d)
mk=m; break;
end
m=m+1;
end
x0=x0+rho^mk*d;
val=feval(fun,x0);
g0=g; d0=d;
fprintf('迭代次数：%d 误差：%.20f 极值点：(x,y) = (%f,%f) 极值：f(x,y) = %.20f\n',k,vpa(norm(g)),x0(1),x0(2),vpa(val));
k=k+1;
end
x=x0;
val=feval(fun,x);
end
运行结果
迭代次数：0 误差：67.08203932499368704612 极值点：(x,y) = (-0.007770,0.496115) 极值：f(x,y) = 2.96992196514529815943
迭代次数：1 误差：12.09283402266935070202 极值点：(x,y) = (-0.051922,0.130463) 极值：f(x,y) = 0.26867338102736537664
迭代次数：2 误差：3.99202899953111067788 极值点：(x,y) = (0.012581,0.002958) 极值：f(x,y) = 0.00390773472646663463
迭代次数：3 误差：0.60898932445470077557 极值点：(x,y) = (0.002422,0.001766) 极值：f(x,y) = 0.00017823014834137801
迭代次数：4 误差：0.12377341781553205524 极值点：(x,y) = (0.000049,0.001005) 极值：f(x,y) = 0.00001216966817470481
迭代次数：5 误差：0.02422649965973680888 极值点：(x,y) = (-0.000168,0.000281) 极值：f(x,y) = 0.00000162869643094180
迭代次数：6 误差：0.01052822247166376388 极值点：(x,y) = (0.000058,0.000092) 极值：f(x,y) = 0.00000018247934170596
迭代次数：7 误差：0.00355209861285059932 极值点：(x,y) = (0.000006,0.000009) 极值：f(x,y) = 0.00000000174914745315
迭代次数：8 误差：0.00035156057294959651 极值点：(x,y) = (-0.000003,0.000002) 极值：f(x,y) = 0.00000000020403920753
迭代次数：9 误差：0.00013092395308230562 极值点：(x,y) = (-0.000000,0.000001) 极值：f(x,y) = 0.00000000002390850274
迭代次数：10 误差：0.00003772948173371275 极值点：(x,y) = (0.000000,0.000001) 极值：f(x,y) = 0.00000000000544374612
迭代次数：11 误差：0.00001636004270145506 极值点：(x,y) = (0.000000,0.000000) 极值：f(x,y) = 0.00000000000056847182
迭代次数：12 误差：0.00000732057953279620 极值点：(x,y) = (0.000000,0.000000) 极值：f(x,y) = 0.00000000000002787083
迭代次数：13 误差：0.00000152498720182152 极值点：(x,y) = (0.000000,0.000000) 极值：f(x,y) = 0.00000000000000235274
x =
1.0e-07 *
0.0039
0.1399
val =
2.3527e-15
k =
14主程序
clc;
clear all;
close all;
% syms x y;%定义函数变量 x y
% f = (1-x)^2+100*(y-x^2)^2;
% dy1 = diff(f,x,1)
% dy2 = diff(f,y,1)
x0 = [0 0]';
[x,val,k]=frcg(@f11,@f22,x0)
x = -2:0.1:2;
y = x;
[X,Y] = meshgrid(x,y);
F =(1-X).^2+100.*(Y-X.^2).^2;
figure;
mesh(X,Y,F);
xlabel('x');
ylabel('y');
zlabel('z');
目标函数自定义函数
function f = f11(X)
x = X(1);
y = X(2);
f = (1-x)^2+100*(y-x^2)^2;
end梯度函数自定义函数
function f = f22(X)
x = X(1);
y = X(2);
f = [2*x - 400*x*(- x^2 + y) - 2 ;
- 200*x^2 + 200*y];
end共轭梯度函数同上
运行结果
x =
1.0000
1.0000
val =
9.4130e-13
k =
892.网上其他的共轭梯度法程序
在网上看到一个写的很好的共轭梯度法的程序，结果是动态的，分享一下。
程序
%来源于知乎作者：达瓦里希也喝脉动的文章
%《最速下降法、牛顿法、共轭梯度法》
%https://zhuanlan.zhihu.com/p/510421800
% 共轭梯度法
clear; clc; close;
%%
% 绘图
a=1; b=100;
x = -1.2:0.03:1.2; y = -0.1:0.03:1.2;
[X,Y] = meshgrid(x,y); F = (a-X).^2 + b*(Y-X.^2).^2;
surf(X, Y, F);
title('CG Method for Rosenbrock'); xlabel('X'); ylabel('Y'); zlabel('Z');
axis square; hold on;
% 初始值
x=0.2; y=0.8; XY = [x; y];
f=func(x, y); % 目标函数
x_last = x; y_last = y; f_last=func(x_last, y_last);
% 计算一阶导(梯度)和二阶导(海森矩阵)
[fx, fy] = grad(x, y); [fxx, fxy, fyy] = hess(x, y);
G = [fx; fy];
H = [fxx, fxy; fxy, fyy];
for n = 0:100
n = n + 1;
plot3(x,y,func(x,y),'-ro','markersize',8); hold on; % 绘图
pause(1);
% 更新方向d
if n == 1
f=func(x, y);
d = [fx; fy];
else
% 计算一阶导(梯度)和二阶导(海森矩阵)
f = func(x, y); [fx, fy] = grad(x, y); [fxx, fxy, fyy] = hess(x, y);
G = [fx; fy];
H = [fxx, fxy; fxy, fyy];
d = -G + ((d'*H*G) / (d'*H*d)) * d; % 更新方向d, [num, num]
end
% 更新步长
lamda = - (d'*G) / (d'*H*d); % 更新步长lamda, num
XY = XY + lamda * d; x = XY(1); y = XY(2); % 更新x和y, [num, num]
% 画线条
plot3([x_last,x], [y_last,y], [func(x_last,y_last),func(x,y)], 'Linewidth', 1); hold on;
% 判断是否终结循环
delta = abs(func(x,y) - f);
if delta < 1*10^(-10)
break;
end
% 记录这一步(x,y)
x_last = x;
y_last = y;
f_last = func(x_last,y_last);
fprintf('迭代次数:%d, f(%f,%f) = %f.\n',n,x,y,func(x,y));
end
plot(x,y,'-ro','markersize',8,'MarkerFaceColor','r'); % 绘图
text(-1, 1, 200, ['迭代次数:',num2str(n),', 最小f(',num2str(x),',',num2str(y),')=',num2str(func(x,y))],...
'backgroundcolor', [0.95,0.95,0.95]);
% 保存
saveas(gcf, 'CG.png');
%% 子函数
% func1: function 函数
function [f] = func(x, y)
a = 1; b = 100;
f = (a-x)^2 + b*(y-x^2)^2;
end
% func2: gradient of func1 函数一阶导-梯度
function [fx, fy] = grad(x, y)
a = 1; b = 100;
% 梯度
fx = 2*(x-a) + 4*b*x*(x^2-y);
fy = -2*b*(x^2-y);
end
% func3: Hessian of func1 函数二阶导-海森矩阵
function [fxx, fxy, fyy] = hess(x, y)
a = 1; b = 100;
% 梯度
fx = 2*(x-a) + 4*b*x*(x^2-y);
fy = -2*b*(x^2-y);
% 梯度的梯度
fxx = 2 + 4*b*(x^2-y) + 8*b*x^2;
fxy = -4*b*x;
fyy = 2*b;
end运行结果
迭代次数:1, f(0.531166,-0.006687) = 8.561771.
迭代次数:2, f(0.617559,0.099603) = 8.086014.
迭代次数:3, f(0.619917,0.384313) = 0.144463.
迭代次数:4, f(1.001153,0.856950) = 2.112874.
迭代次数:5, f(1.126771,1.158013) = 1.261516.
迭代次数:6, f(1.063372,1.118446) = 0.019181.
迭代次数:7, f(1.053420,1.097051) = 0.018838.
迭代次数:8, f(1.038044,1.077288) = 0.001453.
迭代次数:9, f(1.002693,1.003429) = 0.000394.
迭代次数:10, f(0.994178,0.987765) = 0.000073.
迭代次数:11, f(1.000396,1.000598) = 0.000004.
迭代次数:12, f(0.999438,0.998840) = 0.000000.
迭代次数:13, f(1.000008,1.000014) = 0.000000.
迭代次数:14, f(0.999996,0.999991) = 0.000000.
>>3.参考内容
[1] 知乎作者达瓦里希也喝脉动 的文章《最速下降法、牛顿法、共轭梯度法》：https://zhuanlan.zhihu.com/p/510421800
本文内容来源于网络，仅供参考学习，如内容、图片有任何版权问题，请联系处理，24小时内删除。
作 者 | 郭志龙
编 辑 | 郭志龙
校 对 | 郭志龙