matlab subs函数_工程优化设计与Matlab实现——无约束问题的导数解法（共轭梯度法）...

“轭”是牛拉车用的木头，同时拉一辆车的两头牛，就是“共轭”关系

共轭梯度法借助了梯度方向来构造共轭方向，综合了最速下降法和牛顿法的的优点

共轭梯度法

Point 1 共轭的概念

“共轭”是一个比较文艺的词，翻译成大白话就是“成对”。准确地说，是按某种关系成对。

其定义是这样的，G为一个对称正定矩阵

如果有两个向量

和

，使

，则称

和

共对矩阵

共轭

如果有一组向量

，满足

，则这组向量共轭。

当G为单位矩阵时，有

，所以正交是共轭的一种特殊情况，共轭可以看作是

经过

线性变换后，与

正交。

Point 2 共轭方向与函数极值的关系

由于目标函数可以使用在任一点

二次泰勒展开式

来近似，

从而转换为求

这一n元二次函数极值点的问题。

所以我们可以先以二元二次函数为例进行说明：

对于二元二次函数

，

任选初始点

和任意搜索方向

，作一维搜索得

，即

如果

为一维搜索的最优步长，那么

对

的一阶导数为0，即：

，这说明了搜索方向

与

的梯度方向垂直，如下图。

假设从

点继续搜索得到的

是目标函数的极小值点，则在极小值点处需要满足的必要条件为：

由此可得：

两边同时乘以

得：

又由

可得：

，即

与

共对矩阵

共轭。

这就说明沿着与

共轭的方向搜索可以更快的找到极小值点。

所以回到在一般的目标函数中，共轭条件是：

Point 3 如何构造共轭方向

在上图中我们可以将

的搜索方向

视为其负梯度方向

与

的搜索方向

的线性组合，即

将目标函数在任意点

进行泰勒二次展开

，并求偏导得

又由

所以

上式左右同乘

，有：

由于我们要求相邻两次搜索方向

和

满足共轭条件：

所以有：

将

代入，有：

由于之前Point 2中证明过的

，所以

由于各迭代点均为一维优化点，所以个点的梯度矢量均正交（证明略），即

故有：

移项得

如果上一次迭代时，初始点记为

，迭代点记为

， 两个迭代方向分别记为

和

故有：

;

所以：

所以：

所以：

Point 4 共轭梯度法的步长

当我们在给定的初始点，有了明确的迭代方向，那么又回到了在给定点和给定方向上找到最优步长的问题。跟工程优化设计与Matlab实现——一维搜索方法（总结）中提到的一维搜索解决多维问题，工程优化设计与Matlab实现——无约束问题的导数解法（一）中最速下降法求最优步长，工程优化设计与Matlab实现——无约束问题的导数解法（三）中阻尼牛顿法求最优步长一模一样。仍然是用一维搜索的方法。

当然根据Point3中的推导公式，也可以直接求出步长的表达式:

由Point 3中的

左右同时乘以

有：

因为

，所以：

即求得：

取每次的展开点

为

，则有

Point 5 共轭梯度法的特点

共轭梯度法实际上给出了相邻两次迭代方向之间的关系，由此对牛顿法的求Hessian矩阵（二次偏导）进行了规避，只需要求一次偏导，减小了计算量。

同时相对于最速下降法来说，共轭梯度法不仅考虑了当前点的下降速度，也考虑了相邻两次迭代方向之间的关系，使之有了更好的精度。

Point 6 栗子

利用共轭梯度法求目标函数

，在初始点为

时的极小值点和极小值

主程序：

clc

目标函数定义：

function

共轭梯度法函数定义（写的有些繁琐了，可以简化）：

%-----------共轭梯度法---主函数-----------%
function [GongErTiDuFa_x, GongErTiDuFa_xf, GongErTiDuFa_n] = GongErTiDuFa(e1, x0)
    syms t1 t2 t3 t4;
    x = [t1;t2];
    f = func(x);

    %求x0处jacobian矩阵、hessian矩阵的值
    yf = jacobian(f,x)';       %求jacobian矩阵的转置（梯度）
    hf = hessian(f,x);         %求hessian矩阵

    %第一次迭代
    n = 1;                     %迭代次数统计     
    yf_x0 = subs(yf, x, x0);     %将x0代入，得x0处的梯度
    hf_x0 = subs(hf, x, x0);     %将x0代入，得x0处的hessian矩阵
    d0 = -yf_x0;                                  %x0迭代方向
    yf_x0_mo = norm(yf_x0);                       %x0处梯度的模
    a0 = (yf_x0' *  yf_x0) / (d0' * hf_x0 * d0);  %x0迭代步长

    x1 = x0 + a0 * d0;                                 %求x1    
    f_x1 = subs(f, x, x1);            %x1的函数值
    yf_x1 = subs(yf, x, x1);         %x1的梯度值

    while norm(yf_x1) > e1
        n = n + 1;
        hf_x1 = subs(hf, x, x1);          %x1处的hessian矩阵
        yf_x1_mo = norm(yf_x1);            %x1处梯度的模

        b1 = (yf_x1' *  yf_x1) / (yf_x0' *  yf_x0);
        d1 = - yf_x1 + b1 * d0;          %x1迭代方向 
        a1 = (yf_x1' *  yf_x1) / (d1' * hf_x1 * d1);      %x1迭代步长

        x2 = x1 + a1 * d1;                                 %求x2
        f_x2 = subs(f, x, x2);            %x1的函数值
        yf_x2 = subs(yf, x, x2);         %x1的梯度值

        if norm(yf_x2) <= e1
            break
        else
            x1 = x2;
        end
    end

    if n ==1
        x = vpa(x1);
        xf = vpa(f_x1);
    else
        x = vpa(x2);
        xf = vpa(f_x2);
    end

    GongErTiDuFa_x = vpa(x);
    GongErTiDuFa_xf = vpa(xf);
    GongErTiDuFa_n = n;

使用一维搜索求步长，只需将求步长的公式替换成一维搜索的两个函数即可，不再赘述。

计算结果：

---

仔细品味，共轭梯度法相当于如下操作：

先像牛顿法那样在每次迭代目标函数的初始点处进行二次泰勒展开，把目标函数的变成二次问题（如果是二维二次问题其等值线就是一组椭圆）

然后再像最速下降法那样，按照根据梯度建立的共轭方向，去找极小值点作为下一次迭代点。

这可能就是它综合了二者的优点的原因吧。