优化算法：最速下降法、梯度下降法、共轭梯度法---理论分析与实践

优化算法：最速下降法、梯度下降法、共轭梯度法—理论分析与实践

最近被TRPO算法搞得头大，对于其中用到的共轭梯度法等优化算法也不了解，在此详细分析一下最速下降法、梯度下降法、共轭梯度法三者的异同。

参考：
数值分析6(3共轭梯度法)
共轭梯度法简介

0、优化目标

考虑最小化二次函数：
$$?(x)=\frac{1}{2}{x}^{T}Ax?x^{T}b$$
其中 $b,x\in {\mathbb{R}}^n,A\in {\mathbb{R}}^{n\times n}$且假设矩阵$A$是对称正定的（SPD）。该函数的最小值$x^{?}$可以根据一阶最优条件得到，即导数为零
$$??(x^{?})=A{x}^{?}?b=0$$
这也意味着最小化$?(x)$等价于求解线性方程$Ax=b$。此外，在TRPO算法中，由于二次函数的Hessian矩阵是半正定的，该解具有唯一性。

1、最速下降法

最速下降法的策略是：在任何给定点$x$中，函数$?(x)$的负梯度给出的搜索方向是最速下降的方向。换句话说，负梯度方向是局部最优的搜索方向，如下图A点所标注的就是负梯度方向。注意对于二次函数而言它的梯度为$Ax?b$，我们也将它称为系统的残差$r$

有了方向，那么该如何知道要沿着这个方向走多远呢，最速下降法给出的结论是走到函数值不再下降为止，也就是如果继续走下去函数值会增加。上面的概念容易使人混淆，给出下面的一个优化案例进行解释。

下面是使用最速下降法得到的优化过程图像，其中心位置为最小值点，大致的搜索方向是由右上到中心，每一次转折点就代表已经走到函数值不再下降为止，需要重新进行方向搜索。至于为什么已经走到函数值不再下降为止呢，图中有直观的解释，可以看到每个转折点的原来的曲线都是和等高线圆相切，并且重新进行了方向搜索之后的方向都是和等高线切线垂直（也就是当前局部下降最快的方向）。

由此可见，最速下降法是个蛮牛一样的角色，每次迭代是不撞南墙不回头的类型（走到函数值不再下降为止），这种策略有助于减少总的迭代次数，但是走出的路线可能比较曲折。

回到正题，如何确定每次要走多远呢？也就是要走多远才能达到南墙（走到函数值不再下降为止）。回顾所有的已知条件：

1. 迭代式子为：$x_{k+1}=x_k?{\alpha }_k??(x_k)$，其中${\alpha }_k$为步长。
2. 下降方向为：$?f(x_k)=A_{k}x?b=r_k$，这个式子是对优化函数$?(x)=\frac{1}{2}{x}^{T}Ax?x^{T}b$求导得到的。

我们的目标就是计算出走到函数值不再下降为止的${\alpha }_k$，也就是找到函数$?(x)$关于变量${\alpha }_k$的极值点。首先对于${\alpha }_k$，我们有：
$$\begin{array}{l}{\alpha }_k=\arg {\min }_{{\alpha }_k}?(x_{k+1})\\ =\arg {\min }_{{\alpha }_k}?(x_k?{\alpha }_k??(x_k))\\ =\arg {\min }_{{\alpha }_k}?(x_k?{\alpha }_k{r}_k)\\ =\arg {\min }_{{\alpha }_k}\left[\frac{1}{2}{x}_{k+1}^{T}A{x}_{k+1}?x_{k+1}^{T}b\right]\end{array}$$

对$?(x_{k+1})$关于变量${\alpha }_k$求导并使得其导数为0（此处使用复合函数求导公式，因为$x_{k+1}$中含有${\alpha }_k$，得：

$$\begin{array}{l}\frac{??(x_{k+1})}{?{\alpha }_k}=\frac{?\left[\frac{1}{2}{x}_{k+1}^{T}A{x}_{k+1}?x_{k+1}^{T}b\right]}{?{\alpha }_k}\\ =?r_k^{T}A{x}_{k+1}+r_k^{T}b\\ =?r_k^{T}A\left(x_k?{\alpha }_k{r}_k\right)+r_k^{T}b\\ =?r_k^{T}A{x}_k+r_k^{T}A{\alpha }_k{r}_k+r_k^{T}b\\ =?r_k^T\left(A{x}_k?b\right)+r_k^{T}A{\alpha }_k{r}_k\\ =?r_k^T{r}_k+r_k^{T}A{\alpha }_k{r}_k\\ =0\end{array}$$
由此可得最佳步长的计算结果为：
$${\alpha }_k=\frac{r_k^T{r}_k}{r_k^{T}A{r}_k}$$

方便理解，给出了Matlab编程代码的实现：

%% linear equation Ax=b; Ax-b=0
% min(phi(x)=0.5*x'Ax-x'*b)
clear
clc
A = [3,-2;-2,4];
b = [0;-2];

%% 最速下降法
x0 = [10;10];
x_buffer(1,:) = x0;

iter_max = 1000;
phi = zeros(1,iter_max+1);
phi(1) = 0.5*x0'*A*x0 - x0'*b;

r = [0;0];
for i = 1:iter_max
    r_old = r;
    % 计算残差，梯度下降方向
    r = A*x0 - b;
    % 计算最佳步长
    alpha = (r'*r)/(r'*A*r);
    % alpha = 0.1;
    % 执行最速下降
    x = x0 - alpha*r;
    % disp(r'*r_old)
    if norm(x-x0)<=10^(-8)
        break
    end
    x0 = x;
    x_buffer(i+1,:) = x0;
    phi(i+1) = 0.5*x'*A*x - x'*b;
end
phi = phi(1:i);
iter = i-1;

close all



figure
t=-10:.02:10;
[x,y]=meshgrid(t,t);%形成格点矩阵

for i= 1:1:length(t)
    for j= 1:1:length(t)
        x_tmp=[x(i,j);y(i,j)];
        phi_mesh(i, j) = 0.5*x_tmp'*A*x_tmp-x_tmp'*b;
    end
end
mesh(x,y,phi_mesh);
hold on
plot3(x_buffer(:,1),x_buffer(:,2),phi,'LineWidth',2,'Color','r');
% axis([-0.5 0.5 -0.5 0.5 -2 2]);
title('mesh')
colormap summer%cool是一种配色方案，还有其他方案如winter，summer····见help colormap
colorbar

figure
contour(x,y,phi_mesh,phi,'ShowText','on');
hold on
plot(x_buffer(:,1),x_buffer(:,2),'LineWidth',2,'MarkerSize', 32,'Color','r');
disp(['iter: ',num2str(iter),'       error:',num2str(sum(r))])

运行结果如下所示，迭代了31次得到了小于10^(-8)的误差：

2、梯度下降法

之前也提到了，最速下降法是个蛮牛一样的角色，每次迭代是不撞南墙不回头的类型（采用走到函数值不再下降为止的步长${\alpha }_k$），这种策略有助于减少总的迭代次数，但是走出的路线可能比较曲折。

在机器学习算法中，我们经常使用梯度下降法进行训练，此时，我们不再依赖每部迭代去获得步长${\alpha }_k$，而是往往直接指定一个超参数的学习率Lr，修改上述的Matlab代码，不再使用最速下降法的步长计算式（也就是强制alpha = 0.1）：

for i = 1:iter_max
    r_old = r;
    % 计算残差，梯度下降方向
    r = A*x0 - b;
    % 计算最佳步长
    % alpha = (r'*r)/(r'*A*r);
    alpha = 0.1;
    % 执行最速下降
    x = x0 - alpha*r;
    % disp(r'*r_old)
    if norm(x-x0)<=10^(-8)
        break
    end
    x0 = x;
    x_buffer(i+1,:) = x0;
    phi(i+1) = 0.5*x'*A*x - x'*b;
end

得到的结果如下所示（可以看到迭代次数显著变多，从31次变为124次，但是曲线变得平滑）：

3、共轭梯度法

3.1、最速下降法劣势

共轭梯度法主要是用于解决最速下降法的收敛慢的问题，那么首先，这个问题从何而来呢？事实上，最速下降法的迭代方向存在限制，相邻两次迭代的方向必然垂直（正交），由此其方向的选择收到约束，必须走之字形，也就是如下的形状：

相邻两次迭代的方向必然垂直的证明也非常简单，过程如下：

设相邻两次的迭代方向为：
$$\begin{array}{l}{r}_k=A{x}_k?b\\ r_{k+1}=A{x}_{k+1}?b\end{array}$$
根据更新关系可得：
$$\begin{array}{l}{r}_k=A{x}_k?b\\ r_{k+1}=A{x}_{k+1}?b\\ x_{k+1}=x_k?{\alpha }_k{r}_k=x_k?{\alpha }_k\left(A{x}_k?b\right)\end{array}$$

由此：
$$\begin{array}{l}{r}_{k+1}=A{x}_{k+1}?b=A\left[x_k?{\alpha }_k{r}_k\right]?b\\ =A{x}_k?b?A{\alpha }_k{r}_k=r_k?A{\alpha }_k{r}_k\end{array}$$
结合${\alpha }_k=\frac{r_k^T{r}_k}{r_k^{T}A{r}_k}$计算两次相邻梯度的积：
$$r_k^T{r}_{k+1}=r_k^T{r}_k?{\alpha }_k{r}_k^{T}A{r}_k=r_k^T{r}_k?\frac{r_k^T{r}_k}{r_k^{T}A{r}_k}{r}_k^{T}A{r}_k=0$$
积为0，则代表相邻两次迭代的方向必然垂直，这样搜索必然是歪歪扭扭的形状。

3.2、共轭的含义

要理解共轭梯度法，首先要知道共轭的含义，和负数中那种虚部相反的概念不同，此处的共轭的定义更加偏向于A正交，其定义如下：

回到优化目标上，函数$?(x)$的最小值$x^{?}$可以根据一阶最优条件得到。
即导数为零$??(x^{?})=A{x}^{?}?b=0$。我们求解的本质就是$A{x}^{?}?b=0$这个线性方程组，也就是如下的形式：
$$\left[\begin{array}{cccc}|& |& ?& |\\ a_1& a_2& ?& a_n\\ |& |& ?& |\end{array}\right]\left[\begin{array}{c}{x}_1\\ ?\\ x_n\end{array}\right]=\left[\begin{array}{c}|\\ b\\ |\end{array}\right]$$

3.3、A共轭向量构成变量x

假设给定空间$R^n$中一组彼此A共轭的向量组$P_1,P_2,?,P_n$，那我们可以用这个向量组来表示整个$R^n$空间的数值吗？只需要证明$P_1,P_2,?,P_n$线性无关即可。给定一组非零向量：$$p_1,p_2,...,p_k,p_j\in R^n,\&k\le n$$
若：
$$(p_i,p_j{)}_A=0,i\ne j$$
线性无关等价于仅当${\alpha }_l=0$时${\sum }_{l=1}^k{\alpha }_l{p}_l=0$成立。那么下面我们看看能否从${\sum }_{l=1}^k{\alpha }_l{p}_l=0$推导出${\alpha }_l=0$，推导过程非常简单，在两边乘以$p_{n}A$，再使用A共轭的性质进行化简约去求和的无关项即可（$(p_i,p_j{)}_A=0,i\ne j$）

$$\begin{array}{l}{p}_{1}A\sum _{l=1}^k{\alpha }_l{p}_l=0?{\alpha }_1{p}_{1}A{p}_1=0?{\alpha }_1=0\\ ...\\ p_{n}A\sum _{l=1}^k{\alpha }_l{p}_l=0?{\alpha }_n{p}_{n}A{p}_n=0?{\alpha }_n=0\end{array}$$

由此可得空间$R^n$中一组彼此A共轭的向量组$P_1,P_2,?,P_n$是线性无关的，那么$P_1,P_2,?,P_n$可以用于表示空间中的所有内容，使用$P_1,P_2,?,P_n$来表示变量要求解的变量$x$，可得：$$x=\sum _{j=1}^n{\alpha }_j{P}_j?Ax=b=\sum _{j=1}^n{\alpha }_{j}A{P}_j$$
同样的，利用性质$(p_i,p_j{)}_A=0,i\ne j$，我们在方差左右两边乘以$P_k^T$，即可约去求和的无关项，得到：
$$P_k^{T}b={\alpha }_k{P}_k^{T}A{P}_k$$
由此可计算：
$${\alpha }_k=\frac{(P_k,b)}{{\left(P_k,P_k\right)}_A}=\frac{P_k^{T}b}{P_k^{T}A{P}_k}$$
那么，如果每个${\alpha }_k$可以通过计算得出，那么变量$x$的求解岂不是也是易如反掌了嘛。

由此可见，如果能够给定空间$R^n$中一组彼此A共轭的向量组$P_1,P_2,?,P_n$，对于方程的求解可以无需迭代，直接进行。

3.4、A共轭向量$P_1,P_2,?,P_n$的构建

之前也说到了，A共轭也叫A正交，是一种特殊的正交形式。要获得维度为n的正交向量，我们可以使用Smith方法，当然此处是A正交，因此其Smith正交化公式有些许不同：
$$?\begin{array}{l}{P}_0=y_0\\ P_j=y_j?{\sum }_{l=0}^{j?1}\frac{(A{P}_l,y_j)}{(A{P}_l,P_l)}{P}_l,j\ge 1\end{array}$$
通过将原本的残差$r$转化为A正交的向量组$P_n$，原有的方程可以被轻松解决。

但是之前也提到了，如果能够给定空间$R^n$中一组彼此A共轭的向量组$P_1,P_2,?,P_n$，对于方程的求解可以无需迭代，直接进行。那么为什么在实际运用中共轭梯度法是一种迭代求解的算法呢？

简单来说，每次迭代我们都会获得一个残差$r_t$，对这个残差$r_t$进行Smith方法的A正交化，由此才可以得到一个所对应的A共轭的向量$P_t$。也就是说，如果迭代了n步，我们一定能找到空间$R^n$中完整的A共轭的向量组$P_1,P_2,?,P_n$，由此可以得到$x$的精确解。

但是，如果迭代少于n步时，倘若我们的误差已经满足了要求，我们可以不继续进行迭代。因此，共轭梯度法可以使用不超过n的迭代次数对方程进行精确的求解。共轭梯度法也是个蛮牛一样的角色，每次迭代是不撞南墙不回头的类型（采用走到函数值不再下降为止的步长${\alpha }_k$）。算法流程如下所示：

这实际上是一个子空间的问题，对于要求解的目标$x$，其处于空间$R^n$中，我们需要A共轭的向量组$P_1,P_2,?,P_n$对其进行精确的求解。但是，我们也可以使用少于n个向量构建的子空间对其最终结果进行逼近。

打个比方，我要造一个房子，需要钢筋、水泥、石灰、地板、地砖、洗衣机、电饭煲、床，但是一开始我两手空空，需要去采购一番（相当于迭代），我采购了5个材料并进行了施工（相当于迭代了五次），分别是钢筋、水泥、石灰、地板、地砖。虽然没有其他更好的东西，但是我的房子已经基本符合我的预期了，那我为了省钱可能就不去买洗衣机、电饭煲、床等等了。当然，我也可以再买一个洗衣机，这样房子会更加符合我的预期。

3.5、最速下降法实现

下面是结果，争对n=2的空间进行求解，迭代2次实现了10^-15的误差，实际上是精确求解的范畴，符合共轭梯度法可以使用不超过n的迭代次数对方程进行精确的求解

%% linear equation Ax=b; Ax-b=0
% min(phi(x)=0.5*x'Ax-x'*b)
clear
clc
A = [3,-2;-2,4];
b = [0;-2];


%% 共轭梯度法
x0 = [10;10];
x_buffer(1,:) = x0;
iter_max = 1000;
phi = zeros(1,iter_max);
phi(1) = 0.5*x0'*A*x0 - x0'*b;

r0 = A*x0 - b;
p0 = -r0;
for i = 1:iter_max
    alpha = (r0'*r0)/(p0'*A*p0);
    x = x0 + alpha*p0;
    r = r0 + alpha*A*p0;
    beta = (r'*r)/(r0'*r0);
    p = -r + beta*p0;
    if norm(x-x0)<=10^(-8)
        break
    end
    x0 = x;
    r0 = r;
    p0 = p;

    x_buffer(i+1,:) = x0;
    phi(i+1) = 0.5*x'*A*x - x'*b;
end
iter = i-1;

phi = phi(1:i);

close all
figure
t=-10:.1:10;
[x,y]=meshgrid(t,t);%形成格点矩阵

for i= 1:1:length(t)
    for j= 1:1:length(t)
        x_tmp=[x(i,j);y(i,j)];
        phi_mesh(i, j) = 0.5*x_tmp'*A*x_tmp-x_tmp'*b;

    end
end
mesh(x,y,phi_mesh);
hold on
plot3(x_buffer(:,1),x_buffer(:,2),phi,'LineWidth',2,'Color','r');
% axis([-0.5 0.5 -0.5 0.5 -2 2]);
title('mesh')
colormap summer%cool是一种配色方案，还有其他方案如winter，summer····见help colormap
colorbar

figure
contour(x,y,phi_mesh,phi);
hold on
plot(x_buffer(:,1),x_buffer(:,2),'LineWidth',2,'Color','r');

disp(['iter: ',num2str(iter),'       error:',num2str(sum(r0))])