问题描述
二次函数的形式如下
设问题的维度n=158,取初始点为全为0的n维向量。
由于问题的形式特殊,所以步长α采用精确线搜索的显示表示。
对于函数的参数G,b,采用随机生成。
函数文件
function [fun,grad,Hess,b]=f(x)
n=158;
a=unidrnd(10,n,1);
G=a*a'+unidrnd(2)*eye(n);
b=0.5*G*ones(n,1);
grad=G*x-b; %由于函数的系数是随机生成的,所以在下面程序中的迭代都是采用该方法
Hess=G;
fun=0.5*x'*G*x-b'*x;
end
最速下降法
function [x_star]=steepest(x0,eps)
k=0;
xk=x0; %选择初始点
maxit=100000; %最大迭代次数
[~,gradk,G,b]=f(xk); %设定相关参数
while 1
res=norm(gradk); %gradk的二范数
k=k+1; %第k次迭代
if res<=eps||k>maxit %算法停止条件
if k>maxit
fprintf('');
end
break;
end
dk=-gradk; %确定下降方向
alphak=(-dk'*gradk)/(dk'*G*dk); %计算步长
xk=xk+alphak*dk; %迭代
gradk=G*xk-b;
fprintf('At the %d-th iteration,the residual is %f\n',k,res);
end
x_star=xk;
牛顿法
function [x_star]=newton(x0,eps)
k=0;
xk=x0; %选择初始点
maxit=1000000; %最大迭代次数
[~,gradk,G,b]=f(xk); %设定相关参数
while 1
res=norm(gradk); %gradk的二范数
k=k+1; %第k次迭代
if res<=eps||k>maxit %算法停止条件
break;
end
dk=-inv(G)*gradk; %牛顿法确定下降方向
alphak=(-dk'*gradk)/(dk'*G*dk); %计算步长
xk=xk+alphak*dk; %迭代
gradk=G*xk-b;
fprintf('At the %d-th iteration,the residual is %f\n',k,res);
end
x_star=xk;
拟牛顿法
function [x_star]=bfgs(x0,eps)
k=0;
xk=x0; %选择初始点
maxit=100000; %最大迭代次数
Hk=eye(length(x0)); %Hk矩阵初始化
I=eye(length(x0));
[~,gradk,G,b]=f(xk); %设定相关参数
while 1
res=norm(gradk); %gradk的二范数
k=k+1; %第k次迭代
if res<=eps||k>maxit %算法停止条件
break;
end
dk=-Hk*gradk; %确定下降方向
alphak=-dk'*gradk/(dk'*G*dk);%计算步长
xk=xk+alphak*dk; %迭代核心
%更新Hk矩阵,利用Sherman-Morrison-Woodbury公式
grad_temp=gradk; %计算yk
gradk=G*xk-b;
yk=gradk-grad_temp;
deltak=G\yk; %计算deltak
temp1=I-(deltak*yk')/(deltak'*yk); %更新Hk矩阵
temp2=(deltak*deltak')/(deltak'*yk);
Hk=temp1*Hk*temp1+temp2;
fprintf('At the %d-th iteration,the residual is %f\n',k,res);
end
x_star=xk;
共轭梯度法
function [x_star]=fr(x0,eps)
k=0;
xk=x0; %选择初始点
maxit=10000; %最大迭代次数
[~,gradk,G,b]=f(xk); %设定相关参数
while 1
res=norm(gradk); %gradk的二范数
k=k+1; %第k次迭代
if res<=eps||k>maxit %算法停止条件
break;
end
dk=-gradk; %确定下降方向
alphak=-gradk'*dk/(dk'*G*dk); %计算步长
xk=xk+alphak*dk; %迭代核心
%计算Beta
grad_temp=gradk;
gradk=G*xk-b;
Betak=gradk'*gradk/(grad_temp'*grad_temp);
dk=-gradk+Betak*dk;
fprintf('At the %d-th iteration,the residual is %f\n',k,res);
end
x_star=xk;
总结
运行结果如下:
>> steepest(zeros(158,1),1e-4)
At the 1-th iteration,the residual is 34754.134481
At the 2-th iteration,the residual is 3.002778
At the 3-th iteration,the residual is 1.633110
At the 4-th iteration,the residual is 0.000141
>> newton(zeros(158,1),1e-4)
At the 1-th iteration,the residual is 37774.400535
>> bfgs(zeros(158,1),1e-4)
At the 1-th iteration,the residual is 34511.098765
At the 2-th iteration,the residual is 5.758404
>> fr(zeros(158,1),1e-4)
At the 1-th iteration,the residual is 27961.323766
At the 2-th iteration,the residual is 6.037534
At the 3-th iteration,the residual is 3.306320
At the 4-th iteration,the residual is 0.000714
At the 5-th iteration,the residual is 0.000391
最后极小点为一个全是0.5的n维向量。
其实根据函数的构造来看,尽管系数是随机的,但是对于极小点有Gx=b,解就很明显了。
牛顿法一次迭代就能得到结果,得益于使用精确线搜索,二次函数具有显示表达。
最速下降法相比其他方法,显然不是“最速”的。。。
拟牛顿法与牛顿法的区别在于迭代矩阵Hk的不同。