DFP算法-MATLAB

背景

DFP算法是在20世纪60年代初期由Davidon、Fletcher和Powell共同开发的，是拟牛顿法（Quasi-Newton Methods）的一种重要实现。拟牛顿法的基本思想是利用目标函数在当前点附近的二次近似来构造搜索方向，并通过迭代更新这一近似来逼近真实的解。DFP算法通过不断更新Hessian矩阵的逆矩阵（或其近似）来实现这一目标，从而避免了直接计算Hessian矩阵及其逆矩阵的高昂代价。

优点

1. 收敛速度快：DFP算法结合了梯度法和牛顿法的优点，能够在保持梯度法简单性的同时，实现更快的收敛速度。特别是在高维问题中，其收敛性能尤为突出。
2. 计算效率高：与牛顿法相比，DFP算法无需直接计算Hessian矩阵及其逆矩阵，从而大大减少了计算量。它只需计算目标函数的一阶导数（梯度），并通过迭代更新Hessian矩阵的逆矩阵的近似值。
3. 对初始点选择无严格要求：DFP算法对初始点的选择没有严格要求，这使得它在实际应用中更加灵活和方便。
4. 稳定性好：通过适当的迭代更新策略和步长选择方法，DFP算法能够保持较好的稳定性，避免算法在迭代过程中发散或陷入局部最优解。

缺点

1. 可能不收敛：对于一般问题，DFP算法所得的点列未必收敛。这主要是由于算法在迭代过程中可能受到多种因素的影响，如步长选择不当、Hessian矩阵逆矩阵的近似误差等。
2. 计算复杂度高：虽然DFP算法避免了直接计算Hessian矩阵及其逆矩阵的高昂代价，但其迭代更新过程中仍需进行大量的矩阵运算和向量运算，这在一定程度上增加了算法的计算复杂度。特别是在高维问题中，这种计算复杂度可能会成为限制算法性能的一个重要因素。
3. 对参数敏感：DFP算法的性能可能受到多个参数的影响，如步长选择参数、迭代终止条件等。这些参数的选取需要根据具体问题的特点进行仔细调整和优化，否则可能会导致算法性能下降或无法收敛。

function [Xmin] = Davidon_Fletcher_Powell(f, X, tol)
% 项目名称：DFP算法  
% 更新时间：2024/09/21  
% 背景： DFP算法是一种秩2拟牛顿法，DFP算法（Davidon-Fletcher-Powell algorithm）是由戴维登（Davidon）于1959年首先提出的一种求解非线性优化问题的算法，后经弗莱彻（Fletcher）和鲍威尔（Powell）的进一步发展和完善。
% 作者：月白风清江有声 

    %定义梯度函数
    syms x1 x2;
    grad_sym = gradient(f(x1, x2), [x1, x2]);
    grad = matlabFunction(grad_sym, 'vars', [x1, x2]);
    
    %初始化
    x = X;  
    H = eye(length(X)); 
    k = 0;

    while norm(grad(x(1), x(2))) >= tol && k < 10 % 添加迭代次数上限防止死循环

        grad_old = grad(x(1), x(2));
        d = -H * grad_old;
        
        %求解步长ak
        syms c;
        x_bar = x + c * d;   
        f_c = f(x_bar(1), x_bar(2));   
        df_dc = diff(f_c, c);
        sol = solve(df_dc == 0, c);
        ak = double(sol);      %   1.不能eval，不是字符串2.不能sol.c，这个不能引用
        
        x_new = x + ak * d;
        grad_new = grad(x_new(1), x_new(2));
        
        %防止除小数
        delta_x = x_new - x;% 2行1列
        if norm(delta_x) < 1e-10  
            warning('Delta X too small.');  
            break;  
        end  
        delta_y = grad_new - grad_old;% 2行1列
        if norm(delta_y) < 1e-10  
            warning('Delta Y too small.');  
            break;  
        end  
        
        % 确保矩阵乘法的维度匹配，一步一步来
        temp1 = delta_x * delta_x'; %  2行2列
        temp2 = H * delta_y;    %  2行1列
        hess_sim = delta_y'*H*delta_y; % 标量
        if hess_sim == 0
            error('Hess_sim too small.');
        else
            num1 = temp1 /(delta_x' *delta_y) ; %2行2列
            num2 = (temp2 * temp2') / (hess_sim); % 2行2列
            H_new = H + num1 - num2;
        end
        
        %更新迭代
        x=x_new;
        H = H_new;
        k = k + 1;
    
    end
    
    Xmin = x_new;

    fprintf('Final min result: k = %.4f, Xmin=[%.4f;%.4f]\n', k, Xmin(1), Xmin(2));
end


%f=@(x1,x2)x1*x1-2*x1*x2+2*x2*x2-2*x1+x2;X=[0;0];tol=10e-3;Davidon_Fletcher_Powell(f, X,tol);