优化方法求二次函数极小点MATLAB实现

问题描述

二次函数的形式如下

设问题的维度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的不同。