《神经网络与机器学习》第4讲最优化基础
神经网络与机器学习
第4章 性能曲面与最优化基础
§4.1 性能曲面函数的泰勒展开与梯度定义
网络性能，比如均方误差，交叉熵，准确率，接收机性能曲线等等，这些有些称为损失函数或代价函数，还有时候统一称为目标函数。从本质上讲它们之间并无区别。将机器学习问题建模成一个最优化问题已经是现在最主流的建模手段，这些都只是最优化问题中优化目标的在各种背景下的别名罢了。这些后面会详细讲解。
网络性能是网络权系数和输入的函数，因此有必要回顾一下网络性能函数曲面和数值优化算法等基础知识。
泰勒级数经常用来在最优解附近的一种近似，为何要近似？因为原函数可能非常复杂，近似后虽然函数值有些不同，但是它们对应的最优解$x^*$非常接近，最优解$x^*$附近的函数曲面也相近，这对神经网络性能研究带来非常大的便利。一元函数的泰勒展开是
\[f(x)\approx f(x^*)+\sum_{n=1}^{\infty}\frac{1}{n!}\frac{d^n}{dx^n}f(x)\Bigl|_{x=x^*}(x-x^*)^n\]
比如函数
\[f(x)=cos(x)\]
在$x^*=0$附近泰勒展开是
\[f(x)\approx 1-\frac{1}{2}x+\frac{1}{24}x^4-\cdots\]
图4.1 余弦函数零点附近的各阶泰勒展开
图4.1可以看出，零点的值各阶展开都非常精确地近似，但是高阶近似函数更加逼近原函数。 实际中网络性能是所有网络参数的多元函数，不是标量函数，设n元函数
\[f(x)=f(x_1,x_2,\cdots,x_n)\]
在一点附近的泰勒级数展开式为
\[f(x)=f(x^*)+\frac{\partial}{\partial x_1}f(x)\Bigl|_{x=x^*}(x_1-x_1^*)+\cdots+\frac{\partial}{\partial x_n}f(x)\Bigl|_{x=x^*}(x_n-x_n^*)+\frac{1}{2}\frac{\partial^2}{\partial x_1^2}f(x)\Bigl|_{x=x^*}(x_1-x_1^*)^2+\cdots+\frac{1}{2}\frac{\partial^2}{\partial x_n^2}f(x)\Bigl|_{x=x^*}(x_n-x_n^*)^2+\frac{1}{2}\frac{\partial^2}{\partial x_1x_2}f(x)\Bigl|_{x=x^*}(x_1-x_1^*)(x_2-x_2^*)+\cdots+\frac{1}{2}\frac{\partial^2}{\partial x_1x_n}f(x)\Bigl|_{x=x^*}(x_1-x_1^*)(x_n-x_n^*)+\cdots\]
这样太繁琐，我们可以定义梯度向量(gradient)
\[\triangledown f(x)=\begin{bmatrix} \frac{\partial f(\bm{x})}{\partial x_1} & \frac{\partial f(\bm{x})}{\partial x_2} & \cdots & \frac{\partial f(\bm{x})}{\partial x_n} \end{bmatrix}^{\mathrm{T}}\]
对应二阶导数的合起来称为Hessian矩阵
\[\triangledown^2f(\bm{x})=\begin{bmatrix} \frac{\partial^2f(\bm{x})}{\partial x_1^2} & \cdots & \frac{\partial^2f(\bm{x})}{\partial x_1x_n}\\ \vdots & \ddots & \vdots\\ \frac{\partial^2f(\bm{x})}{\partial x_nx_1} & \cdots & \frac{\partial^2f(\bm{x})}{\partial x_n^2} \end{bmatrix}\]
那么上面的多元函数泰勒展开可以写成
\[f(x)=f(x^*)+\triangledown f(x)^{\mathrm{T}}\Bigl|_{x=x^*}(x-x^*)+\frac{1}{2}(x-x^*)^{\mathrm{T}}\triangledown ^2f(x)(x-x^*)+\cdots\]
图4.2 函数的切线，方向导数和梯度
梯度是一个矢量，其方向上的方向导数最大，其大小正好是此最大方向导数。比如如图一个二元二次函数$f(x_1,x_2)$，设单位向量$e=cos\theta\bm{i}+sin\theta\bm{j}$ ，在点$(x_1,x_2)$的极限
\[\lim_{t\rightarrow 0}\frac{f(x_1+tcos\theta,x_2+tsin\theta)-f(x_1,x_2)}{t}=\frac{\partial f(\bm{x})}{\partial x_1}cos\theta+\frac{\partial f(\bm{x})}{\partial x_2}sin\theta\]
称为方向导数
\[D_ef(\bm{x})=\bm{e}^{\mathrm{T}}\triangledown f(\bm{x})=\triangledown f(\bm{x})\cdot \bm{e}=|\triangledown f(\bm{x})||\bm{e}|cos\alpha\]
$\alpha$是梯度$\triangledown f(\bm{x})$和单位向量$\bm{e}$的夹角，可以看出，当$\alpha =0$，方向导数和梯度相同方向，函数增长最快，因此梯度其方向上的方向导数最大，其大小正好是最大方向导数！而$\alpha =\pi$，方向导数和梯度方向相反，函数减小最快，而其他角度在二者之间，其中如果垂直的话，函数不增不减（等高线）。
定义：任意方向$p$上的一阶、二阶方向导数为
\[\frac{p^{\mathrm{T}}\triangledown f(x)}{\left \| p \right \|},\frac{p^{\mathrm{T}}\triangledown ^2f(x)p}{\left \| p \right \|^2}\]
$\|p\|$，$\|p\|^2$起到归一化的作用。
极小点：给定函数$f(x)$和点$x^*$，如果存在一个标量$\delta$使得点$x^*$附近的点$x+\Delta x$，且$|\delta x|<\delta$，满足
\[f(x^*)<f(x^*+\Delta x)\]
称为强极小点，而
\[f(x^*)\leq f(x^*+\Delta x)\]
叫弱极小点。当对于所有的点（$|\Delta x|$无穷大），$x=x^*+\Delta x$在整个定义域内，满足$f(x^*)<f(x^*+\Delta x)$，称为全局极小点。
图4.3二元函数$f(x)=(x_1-x_2)^4+8x_1x_2-x_1+x_2+3$的三维等高线，一个局部最小点（-0.42，0.42），一个全局最小点（0.55，-0.55）
§4.2 性能曲面函数优化的充分必要条件
一阶必要条件
$f(x)=f(x^*)+\triangledown f(x)^{\mathrm{T}}\Bigl|_{x=x^*}(x-x^*)+\frac{1}{2}(x-x^*)^{\mathrm{T}}\triangledown ^2f(x)(x-x^*)+\cdots$
再次看上面的泰勒展开，当$|x-x^*|=\Delta x$很小，也就是说$x=x^*+\Delta x$或，$x=x^*-\Delta x$在附近的点$x^*$，用一阶泰勒展开来近似，即
$f(x)=f(x^*)+\triangledown f(x)^{\mathrm{T}}\Bigl|_{x=x^*}(x-x^*)=f(x^*)\pm \triangledown f(x)^{\mathrm{T}}\Bigl|_{x=x^*}\Delta x$
如果$x^*$是对应了极小值最优解，那么必然有
$f(x^*)\pm \triangledown f(x)^{\mathrm{T}}\Bigl|_{x=x^*}\Delta x>f(x^*)$
$\triangledown f(x)^{\mathrm{T}}\Bigl|_{x=x^*}\Delta x=0$
而$\Delta x$是任意的，所以和一元函数求驻点一样，对于网络性能函数最优化的一阶必要条件是
$\triangledown f(x)^{\mathrm{T}}\Bigl|_{x=x^*}=0$
任何满足此条件的点都叫驻点，所以是一阶极小值的必要条件，因为鞍点也满足。
二阶充分条件
如果满足了上面一阶必要条件，我们再看二阶充分条件，那么泰勒展开到二阶
$f(x)=f(x^*)+\triangledown f(x)^{\mathrm{T}}\Bigl|_{x=x^*}(x-x^*)+\frac{1}{2}(x-x^*)^{\mathrm{T}}\triangledown ^2f(x)(x-x^*)$
$=f(x^*)+\frac{1}{2}\Delta x^{\mathrm{T}}\triangledown ^2f(x)\Delta x>f(x^*)$
得到
$\Delta x^{\mathrm{T}}\triangledown ^2f(x)\Delta x>0$
而$\Delta x$是任意的，所以二阶条件是Hessian矩阵正定
$\triangledown ^2f(x)$正定矩阵
其特征值都是正数。
图4.4函数$f(x)=x^2$以及$x=\pm 1$处的切线
我们可以用一个简单的例子来理解上面的这些必要和充分条件，比如$f(x)=x^2$，它的一阶导数"梯度"是$\triangledown f(x)=2x$，得到$x=0$是个驻点，而二阶导数"Hessian矩阵"是$1\times 1$的"矩阵"$\triangledown ^2f(x)=2>0$，仔细看函数图像可以知道一阶导数"梯度"是$\triangledown f(x)=2x$在$x<0$, 也小于0，在$x>0$, 也大于0，并且整个区间逐渐增加的，其变化率就是二阶导数"Hessian矩阵"是$1\times 1$的"矩阵"$\triangledown ^2f(x)>0$。
为了更好地理解这两个最优条件，我们来看一个经典的二次函数(也是很多目标函数的二阶泰勒展开式)
$f(x)=\frac{1}{2}x^{\mathrm{T}}Ax+d^{\mathrm{T}}x+c$
这里矩阵$A$为对称矩阵，可以不对称，那么
$f(x)+(f(x))^{\mathrm{T}}=2f(x)=\frac{1}{2}x^{\mathrm{T}}(A+A^{\mathrm{T}})x+2d^{\mathrm{T}}x+2c$
则$A+A^{\mathrm{T}}$肯定是对称矩阵，因此即使$A$不对称，也能转化过来。那么其梯度为
$\triangledown f(x)=Ax+d$
二阶导数"Hessian矩阵"是
$\triangledown ^2f(x)=A$
二次函数的另外表示形式
$f(x)=\frac{1}{2}x^{\mathrm{T}}Ax+d^{\mathrm{T}}x+c$
其中驻点（这里是对应的全局最小）
$\triangledown f(x^*)=Ax^*+d=0\Rightarrow x^*=-A^{-1}d$
代入得到
$f_{min}=-\frac{1}{2}d^{\mathrm{T}}A^{-1}d+c$
那么和一元二次函数一样进行配方
$f(x)=\frac{1}{2}x^{\mathrm{T}}Ax+d^{\mathrm{T}}x+c=f_{min}+\frac{1}{2}(x-x^*)^{\mathrm{T}}A(x-x^*)$
$f_{min}$起到的作用是函数值的大小改变，而$x-x^*$决定了驻点（取最小值$f_{min}$）位置为$x^*$。因此我们只需要研究一个标准型
$f(x)=\frac{1}{2}x^{\mathrm{T}}Ax$
那么二阶导数"Hessian矩阵"$\triangledown ^2f(x)=A$的特征值和特征向量决定了函数曲面的形状。比如如下两个例子
$f(x)=\frac{1}{2}\begin{bmatrix} x_1 & x_2 \end{bmatrix}\begin{bmatrix} 2 & 0\\ 0 & 2 \end{bmatrix}\begin{bmatrix} x_1\\ x_2 \end{bmatrix}$
其Hessian矩阵$\triangledown ^2f(x)=\begin{bmatrix} 2 & 0\\ 0 & 2 \end{bmatrix}$，特征值和特征向量分别为
$\lambda_1=2$，$v_1=\begin{bmatrix} 1\\ 0 \end{bmatrix}$和$\lambda_2=2$，$v_2=\begin{bmatrix} 0\\ 1 \end{bmatrix}$
|
x = -2:0.1:2; y = x; [X,Y] = meshgrid(x); Z = X.^2+Y.^2; surfc(X,Y,Z) |
图4.5圆形二次曲面
$f(x)=\frac{1}{2}\begin{bmatrix} x_1 & x_2 \end{bmatrix}\begin{bmatrix} 2 & 1\\ 1 & 2 \end{bmatrix}\begin{bmatrix} x_1\\ x_2 \end{bmatrix}$
其Hessian矩阵$\triangledown ^2f(x)=\begin{bmatrix} 2 & 1\\ 1 & 2 \end{bmatrix}$，特征值和特征向量分别为
$\lambda_1=1$，$v_1=\begin{bmatrix} 1\\ -1 \end{bmatrix}$和$\lambda_2=3$，$v_2=\begin{bmatrix} 1\\ 1 \end{bmatrix}$
|
x = -2:0.1:2; y = x; [X,Y] = meshgrid(x); Z = X.^2+Y.^2; surf(X,Y,Z) hold on contour(X,Y,Z,20) |
图4.6椭圆二次曲面
$v_2$方向的等高线更加密，代表了最大曲率。
作业：鞍点和驻点
$f(x)=\frac{1}{2}\begin{bmatrix} x_1 & x_2 \end{bmatrix}\begin{bmatrix} -0.5 & -1.5\\ -1.5 & -0.5 \end{bmatrix}\begin{bmatrix} x_1\\ x_2 \end{bmatrix}$
$f(x)=\frac{1}{2}\begin{bmatrix} x_1 & x_2 \end{bmatrix}\begin{bmatrix} 1 & -1\\ -1 & 1 \end{bmatrix}\begin{bmatrix} x_1\\ x_2 \end{bmatrix}$
***Hessian矩阵的特征系统
将矩阵的特征向量作为列向量并标准化，组成特征矩阵
$\bm{V}=[v_1,v_2,\cdots,v_n]$
由于对称矩阵的特征向量总是相互正交的，所以我们得到
$V^{\mathrm{T}}AV=\begin{bmatrix} \lambda_1 & 0 & 0\\ 0 & \ddots & 0\\ 0 & 0 & \lambda_n \end{bmatrix}=\wedge$
比如$A=\begin{bmatrix} 2 & 1\\ 1 & 2 \end{bmatrix}$特征值和标准化的特征向量分别为
$\lambda_1=1$，$v_1=\begin{bmatrix} \frac{1}{\sqrt 2}\\ - \frac{1}{\sqrt 2} \end{bmatrix}$和$\lambda_2=3$，$v_2=\begin{bmatrix} \frac{1}{\sqrt 2}\\ \frac{1}{\sqrt 2} \end{bmatrix}$
$V^{\mathrm{T}}AV=\begin{bmatrix} \lambda_1 & 0 \\ 0 & 2 \end{bmatrix}$
由于$V$是正交矩阵，即$V^{\mathrm{T}}V=I$，那么还有
$A=V\wedge V^{\mathrm{T}}$
代入标准型
$f(x)=\frac{1}{2}x^{\mathrm{T}}Ax=\frac{1}{2}x^{\mathrm{T}}V\wedge V^{\mathrm{T}}x=\frac{1}{2}(V^{\mathrm{T}}x)^{\mathrm{T}}\wedge V^{\mathrm{T}}x=\frac{1}{2}\widetilde{x}^{\mathrm{T}}\wedge \widetilde{x}$
这里旋转变换
$\widetilde{x}=V^{\mathrm{T}}x$
因为这样变换之后，所有的交叉项都消失了，只剩下平方项。同时，原来的单位向量$[1,0]^{\mathrm{T}}$和$[0,1]^{\mathrm{T}}$也旋转为特征向量的所指示的主向量。
回顾定义：任意方向$p$上的一阶、二阶方向导数为
\[\frac{p^{\mathrm{T}}\triangledown f(x)}{\left \| p \right \|},\frac{p^{\mathrm{T}}\triangledown ^2f(x)p}{\left \| p \right \|^2}\]
$\|p\|$，$\|p\|^2$起到归一化的作用。设
$p=Vc,c=V^{\mathrm{T}}p$
$c$是向量$p$在特征向量$V$上的一个表达，那么二阶方向导数
$\frac{p^{\mathrm{T}}\triangledown ^2f(x)p}{\left \| p \right \|^2}=\frac{p^{\mathrm{T}}Ap}{\left \| p \right \|^2}=\frac{c^{\mathrm{T}}V^{\mathrm{T}}AVc}{\left \| c^{\mathrm{T}}V^{\mathrm{T}}Vc \right \|^2}=\frac{c^{\mathrm{T}}\wedge c}{\left \| c^{\mathrm{T}}c \right \|^2}=\frac{\sum_{i=1}^{n}\lambda_ic_i^2}{\sum_{i=1}^{n}c_i^2}$
因此我们发现
$\lambda_{min}\leq \frac{p^{\mathrm{T}}\triangledown ^2f(x)p}{\left \| p \right \|^2}\leq \lambda_{max}$
二阶方向导数取最大就是对应的特征向量方向$p=v_{max}$。
图4.7椭圆二次曲面等高线和两个特征向量方向
§4.3 性能曲面函数优化算法
这里我们介绍3种优化算法：最速下降法、牛顿下山法和共轭梯度法，这也是训练神经网络的基本算法。虽然这些算法产生于数百年前，但是由于计算机的实现，使得优化理论成了神经网络训练的基础。
我们从一个猜测值$x_0$出发，按照下面的规则进行更新
\[x_{k+1}=x_k+\Delta x_k=x_k+\alpha_{k}p_k\]
标量$\alpha_{k}$是学习率，向量$p_k$是表示搜索方向。最速下降法、牛顿下山法和共轭梯度法的不同在于搜索方向$p_k$不同，当然标量$\alpha_{k}$学习率也有不同选择方式，我们一一在下面进行介绍。
最速下降法：当我们不断更新$x_{k+1}=x_k+\Delta x_k$时，我们希望函数值
$f(x_{k+1})<f(x_k)$
我们再来看在上一次猜测值处的一阶泰勒展开
\[f(x_{k+1})=f(x_k)+\triangledown f(x)^{\mathrm{T}}\Bigl|_{x=x_k}\Delta x_k=f(x_k)+g_k^{\mathrm{T}}\Delta x_k\]
为了表达方便，我们把猜测值$x_k$处的梯度设为
\[g_k=\triangledown f(x)^{\mathrm{T}}\Bigl|_{x=x_k}\]
我们希望函数值$f(x_{k+1})<f(x_k)$越迭代越小，因此$(\alpha_{k}>0)$
\[g_k^{\mathrm{T}}\Delta x_k=\alpha_{k}g_k^{\mathrm{T}}p_k<0\Rightarrow g_k^{\mathrm{T}}p_k<0\]
满足这样的向量$p_k$很多，假设这些方向都是相同的长度（模），$ g_k^{\mathrm{T}}p_k$取最小就是
\[p_k=-g_k\]
在这个方向上下降的函数值最大，因此称为最速下降法（负梯度法）
\[x_{k+1}=x_k+\Delta x_k=x_k-\alpha_{k}g_k\]
例子：最速下降法
$f(x)=\frac{1}{2}\begin{bmatrix} x_1 & x_2 \end{bmatrix}\begin{bmatrix} 2 & 0\\ 0 & 50 \end{bmatrix}\begin{bmatrix} x_1\\ x_2 \end{bmatrix}=x_1^2+25x_2^2$
图4.8学习率$\alpha_{k}=0.01,0.035$的最速下降轨迹
%输入: fun, gfun分别是目标函数及其梯度, x0是初始点,epsilon为容许误差
%输出: k是迭代次数, x, val分别是近似最优点和最优值
maxk=50; %最大迭代次数
x0=[0.5;0.5];
alphak=0.01;
x=zeros(2,maxk+1);
x(:,1)=x0;
for k=1:maxk
gk=gfun(x0); %计算梯度
x0=x0-alphak*gk; %最速下降
x(:,k+1)=x0;
end
x0
val=fun(x0)
dx = -1:0.05:1;
dy = dx;
[dX,dY] = meshgrid(dx);
dZ=dX.^2+25*dY.^2;
contour(dX,dY,dZ,50) %等高线图
hold on
plot(x(1,:),x(2,:),'r-o')
%目标函数
function f=fun(x)
f=x(1).^2+25*x(2).^2;
end
% 梯度
function gf=gfun(x)
gf=[2*x(1); 50*x(2)];
end
--------------------------------------------------------
图3.8给出了两种学习率的最速下降法的轨迹，学习率大了反而出现振荡，因为沿着当前点的负梯度方向走的太过了，那么多大的学习率是合适的呢？这里因为很多函数都用二次函数来近似，那么通过二次函数的Hessian矩阵我们可以给出学习率的一个上界。
对于
\[f(x)=\frac{1}{2}x^{\mathrm{T}}Ax\]
$g_k=\triangledown f(x)^{\mathrm{T}}\Bigl|_{x=x_k}=Ax_k$
那么最速下降法更新公式
\[x_{k+1}=x_k-\alpha_{k}g_k=[I-\alpha_{k}A]x_k\]
前面的旋转变换，
\[\widetilde{x}=V^{\mathrm{T}}x,V\widetilde{x}_k=x_k\]
那么得到
\[V\widetilde{x}_{k+1}=[I-\alpha_{k}A]V\widetilde{x}_{k}\]
对于恒定的学习率$\alpha_k=\alpha$
\[\widetilde{x}_{k+1}=V^{\mathrm{T}}[I-\alpha A]V\widetilde{x}_{k}=[I-\alpha V^{\mathrm{T}}AV]\widetilde{x}_{k}=[I-\alpha \wedge]\widetilde{x}_k=[I-\alpha \wedge]^k\widetilde{x}_0\]
旋转变换不改变最优点的值，仅仅改变的是方向，所以$f(x)=\frac{1}{2}x^{\mathrm{T}}Ax$的最优点（对应最小值）的向量$\widetilde{x}=V^{\mathrm{T}}x=0$，$\widetilde{x}_0=V^{\mathrm{T}}x_0$表示了初始向量的正交变换，是给定的，那么我们得到
\[\lim_{k\rightarrow \infty}[I-\alpha\wedge ]^k=0\]
$I-\alpha \wedge$是对角阵，它的$k$次幂也是对角阵，所以上面为n维差分动力学方程，每个对角线上写成
\[\widetilde{x}_{k,i}=(1-\alpha \lambda_n)^k\widetilde{x}_{0,i}\]
$i=1,2,\cdots,n$。我们希望每迭代一步$\widetilde{x}_{k,i}$都减小，因此要求每一个乘子$(1-\alpha \lambda_n)$满足
\[-1<\left ( 1-\alpha \lambda_n \right )<1\]
才有可能收敛到0。对于正定矩阵$A$，特征值都是正数，那么等价于
\[0<\alpha<\frac{2}{\lambda_{max}}\]
这是最速下降法收敛的充要条件。
*最大学习率与二次函数的最大曲率$\lambda_{max}$，成反比，梯度变化太快，那么学习率则需要变小，防止步长太大跳过最小值点。大家可以把图3.8中的学习率取大于$\alpha<\frac{2}{\lambda_{max}}=0.04$，数值仿真发现不会收敛。
上面学习率$\alpha_k=\alpha$选择的恒定，如果每一步都调整，调整原则是使得下一步的函数值最小，即
\[\min_{\alpha_k}f(x_k-\alpha_kg_k)\]
这里对于函数的泰勒展开
\[f(x_k-\alpha_kg_k)\approx f(x_k)-\triangledown f(x)^{\mathrm{T}}\Bigl|_{x=x_k}\alpha_kg_k+\frac{1}{2}\alpha_kg_k^{\mathrm{T}}\triangledown ^2f(x)\Bigl|_{x=x_k}\alpha_kg_k=f(x_k)-\alpha_kg_k^{\mathrm{T}}g_k+\frac{1}{2}\alpha_k^2g_k^{\mathrm{T}}A_kg_k\]
这里$\triangledown ^2f(x)\Bigl|_{x=x_k}=A_k$。学习率的选择非常简单，就是直接对其进行求导
\[\frac{df(x_k-\alpha_kg_k)}{d\alpha_k}=-g_k^{\mathrm{T}}g_k+\alpha_kg_k^{\mathrm{T}}A_kg_k=0\]
得到
\[\alpha_k=\frac{g_k^{\mathrm{T}}g_k}{g_k^{\mathrm{T}}A_kg_k}\]
所以每一步的学习率都是当前数据$x_k$的梯度$g_k$和Hessian矩阵$A_k$确定的。上面的例子我们重新求解如下，可以看出同样的迭代次数，最优点收敛更快。
图4.9学习率$\alpha_k$最小化的最速下降轨迹
%输入: fun, gfun分别是目标函数及其梯度, x0是初始点,epsilon为容许误差
%输出: k是迭代次数, x, val分别是近似最优点和最优值
maxk=50; %最大迭代次数
x0=[0.5;0.5];
alphak=0.035;
x=zeros(2,maxk+1);
x(:,1)=x0;
Ak=[2 0;0 50];%Hessian均值是常数矩阵
for k=1:maxk
gk=gfun(x0); %计算梯度
alphak=gk'*gk./(gk'*Ak*gk);%每次计算学习率直线最小化
x0=x0-alphak*gk;
x(:,k+1)=x0;
end
x0
val=fun(x0)
dx = -1:0.05:1;
dy = dx;
[dX,dY] = meshgrid(dx);
dZ=dX.^2+25*dY.^2;
contour(dX,dY,dZ,50)
hold on
plot(x(1,:),x(2,:),'r-o')
%目标函数
function f=fun(x)
f=x(1).^2+25*x(2).^2;
end
% 梯度
function gf=gfun(x)
gf=[2*x(1); 50*x(2)];
end
*补充性质：每次迭代方向正交。这可以从导数的链式法则可以看出，因为每一步都是求导数为0
\[\frac{df(x_{k+1})}{d\alpha_k}=0=\triangledown f(x)^{\mathrm{T}}\Bigl|_{x=x_{k+1}}\frac{dx_{k+1}}{d\alpha_k}=g_{k+1}^{\mathrm{T}}(-g_{k})\]
因此我们可以看出相邻两步的梯度是正交的
\[g_{k+1}^{\mathrm{T}}g_k=0\]
牛顿法：最速下降法是利用一阶泰勒展开，牛顿法是利用二阶泰勒展开得出
\[f(x_{k+1})\approx f(x_k)+g_{k+1}^{\mathrm{T}}\Delta x_k+\frac{1}{2}\Delta x_k^{\mathrm{T}}A_k\Delta x_k\]
那么该变量如何使得$f(\bm{x}_{k+1})$最小？对上面展开求导数
\[g_k+A_k\Delta x_k=0\]
我们直接得到
\[\Delta x_k=-A_k^{-1}g_k\]
因此牛顿法的更新规则定义为
\[\Delta x_k=-A_k^{-1}g_kx_{k+1}=x_k-A_k^{-1}g_k\]
但是这种方法局限性非常大，因为它不能分辨极小、极大点和鞍点的区别。
例子：
\[f(x)=\left ( x_1-x_2 \right )^4+8x_1x_2-x_1+x_2+3\]
这个非线性函数有3个驻点
$x=\begin{bmatrix} -0.41878\\ 0.41878 \end{bmatrix}$（局部极小），$\begin{bmatrix} -0.134797\\ 0.134797 \end{bmatrix}$（鞍点），$\begin{bmatrix} 0.55358\\ -0.55358 \end{bmatrix}$（全局极小）
图4.10 非线性函数曲面的三个驻点
我们选择初始点$x_0=\begin{bmatrix} -1.5\\ 0 \end{bmatrix}$，$x_0=\begin{bmatrix} 0.75\\ 0.75 \end{bmatrix}$,$x_0=\begin{bmatrix} 1.15\\ 0.75\end{bmatrix}$在此点求求梯度和Hessian矩阵，发现都不能到达全局最低点，必须对于此方法进行变形，而且它需要求矩阵的逆阵，这导致了计算复杂性。
图4.11 牛顿法的三次初值进行更新的轨迹
dx = -2:0.01:2;
dy = dx;
[dX,dY] = meshgrid(dx);
dZ=(dX-dY).^4+8*dX.*dY-dX+dY+3;
>> contour(dX,dY,dZ,200)
>> hold on
x=-1.5;y=0;
gx = -4*(y-x)^3 + 8*y - 1;
gy = 4*(y-x)^3 + 8*x + 1;
grad = [gx; gy];
% CREATE HESSIAN
temp = 12*(y-x)^2;
hess = [temp 8-temp;8-temp temp];
x1=[x;y]-inv(hess)*grad;
plot([x x1(1)],[y x1(2)],'ro-')
x=0.75;y=0.75;
gx = -4*(y-x)^3 + 8*y - 1;
gy = 4*(y-x)^3 + 8*x + 1;
grad = [gx; gy];
% CREATE HESSIAN
temp = 12*(y-x)^2;
hess = [temp 8-temp;8-temp temp];
x1=[x;y]-inv(hess)*grad;
plot([x x1(1)],[y x1(2)],'bo-')
x=1.15;y=0.75;
gx = -4*(y-x)^3 + 8*y - 1;
gy = 4*(y-x)^3 + 8*x + 1;
grad = [gx; gy];
% CREATE HESSIAN
temp = 12*(y-x)^2;
hess = [temp 8-temp;8-temp temp];
x1=[x;y]-inv(hess)*grad;
plot([x x1(1)],[y x1(2)],'mo-')
注意后面第5章讲反向（BP）传播算法时候，我们会进一步发展牛顿法，一种牛顿法的变形Levenberg-Marquardt算法。
共轭梯度法（**选讲）：最速下降法和牛顿法之间折衷的产物，不需要计算Hessian逆矩阵。适合那些大量参数需要优化，Hessian矩阵逆阵计算复杂的问题。而且上面牛顿法的搜索方向并不是最好的。
共轭向量集合：设对称正定矩阵$A\in R^{n\times n}$，存在n维的向量组${p_i}$满足
\[p_iAp_j=0,i\neq j\]
\[p_iAp_i>0\]
那么向量组${p_i}$是$A$共轭的，且线性无关。比如$A$的特征值和特征向量就是关于$A$共轭的。
我们就是避免计算Hessian矩阵，所以特征值和特征向量这里并没有什么作用，我们试着从初始点的梯度开始构造关于$A$共轭的向量组${p_i}$。考虑正定二次函数
\[f(x)=\frac{1}{2}x^{\mathrm{T}}Ax+d^{\mathrm{T}}x+c\]
的极小点线性共轭方法，
\[x_{k+1}=x_{k+1}+\alpha_kp_k,\alpha_k>0\]
再来回忆$f(x_{k+1})$在上一次猜测值$x_k$处的一阶泰勒展开
\[f(x_{k+1})=f(x_k)+\triangledown f(x)^{\mathrm{T}}\Bigl|_{x=x_k}\alpha_kp_k=f(x_k)+\alpha_kg_k^{\mathrm{T}}p_k\]
为了表达方便，我们把猜测值$x_k$处的梯度设为
\[g_k=\triangledown f(x)^{\mathrm{T}}\Bigl|_{x=x_k}=Ax_k+d^{\mathrm{T}}\]
我们希望函数值$f(x_{k+1})<f(x_k)$越迭代越小，因此要求方向$p_k$满足
\[\alpha_kg_k^{\mathrm{T}}p_k<0\Rightarrow g_k^{\mathrm{T}}p_k<0\]
那么$\alpha_k$按照前面所要求的线性优化，我们更新$x_{k+1}=x_k+\alpha_kp_k$，并计算此点的梯度
\[g_{k+1}=\triangledown f(x)^{\mathrm{T}}\Bigl|_{x=x_{k+1}}=Ax_{k+1}+d^{\mathrm{T}}\]
选取下一步的迭代方向$p_{k+1}$满足下降形状
\[g_{k+1}^{\mathrm{T}}p_{k+1}<0\]
还要满足共轭性
\[p_{k+1}Ap_{j}=0,j=0,1,\cdots,k\]
这样形成了线性共轭的下降方向向量组${p_0,p_1,\cdots,p_{k+1}}$。
但是这种想法还是不能实现，这里引入一个定理：每一步的迭代点$x_{k+1}$都是由初始点和共轭方向向量组所张成的线性流形（就是正交线性组合）
\[S=\left \{ x|x=x_0+\sum_{i=0}^k\alpha_ip_i \right \}\]
中的极小点。我们利用归纳法证明，已构造${p_0,p_1,\cdots,p_k}$线性无关，且共轭$p_kAp_j=0,j=0,1,\cdots,k$。只需要证明点
\[x_{k+1}=x_k+\alpha_kp_k=x_0+\sum _{i=0}^k\alpha_ip_i\in S\]
是$S={x|x=x_0+\sum_{i=0}^k\alpha_ip_i}$中的极小点即可。$S$中任意的一点
\[x=x_0+\sum _{i=0}^k\gamma_ip_i\]
那么与$x_{k+1}$之间的差为
\[h_{k+1}=x-x_{k+1}=\sum_{i=1}^k\left ( \gamma_i-\alpha_i \right )p_i\]
函数在$x=x_{k+1}+h_{k+1}$点的泰勒展开为
\[f(x)=f(x_{k+1})+g_{k+1}^{\mathrm{T}}h_{k+1}+\cdots>f(x_{k+1})+g_{k+1}^{\mathrm{T}}h_{k+1}=f(x_{k+1})+\sum_{i+0}^k\left ( \gamma_i-\alpha_i \right )g_{k+1}^{\mathrm{T}}p_i\]
如果$\sum_{i=0}^k\left ( \gamma_i-\alpha_i \right ) g_{k+1} ^{\mathrm{T}}p_i =0$，那么上面就证明了$x_{k+1}$是\[S=\left \{ x|x=x_0+\sum_{i=0}^k \alpha_ip_i\right \}\]中的极小点。那么我们可以找到一种方式使得
\[\sum_{i=0}^k(\gamma_i-\alpha_i)g_{k+1}^{\mathrm{T}}p_i=0,\Rightarrow g_{k+1}^{\mathrm{T}}p_i=0,i=0,1,\cdots,k\]
对于二次函数
\[g_{k+1}=\triangledown f(x)^{\mathrm{T}}\Bigl|_{x=x_{k+1}}=Ax_{k+1}+d^{\mathrm{T}}\]
\[g_{k+1}-g_k=A(x_{k+1}-x_k)=\alpha_kAp_k\]
那么可以将$i=0,1,\cdots,k$
\[g_1^{\mathrm{T}}p_0=g_1^{\mathrm{T}}p_0\]
\[g_2^{\mathrm{T}}p_0=g_1^{\mathrm{T}}p_0+\left ( g_2-g_1 \right )^{\mathrm{T}}p_0\]
\[g_3^{\mathrm{T}}p_0=g_1^{\mathrm{T}}p_0+\left ( g_2-g_1 \right )^{\mathrm{T}}p_0+\left ( g_3-g_2 \right )^{\mathrm{T}}p_0\]
$\vdots$
\[g_{k+1}^{\mathrm{T}}p_i=g_{i+1}^{\mathrm{T}}p_i+\sum_{j=i+1}^k\left ( g_{j+1}-g_j \right )^{\mathrm{T}}p_i=g_{i+1}^{\mathrm{T}}p_i+\sum_{j=i+1}^k\alpha_jp_j^{\mathrm{T}}Ap_i\]
由于共轭性$\sum_{j=i+1}^k\alpha_jp_j^{\mathrm{T}}Ap_i=0$,只要使得$g_{i+1}^{\mathrm{T}}p_i=0$即可，这可以用类似Gram-Schmidt正交方法实现。
我们一般开始选取初值$x_0$，计算当前梯度$g_0$，选择下降方向为最速下降方向$p_0=-g_0$，那么更新
\[x_1=x_0+\alpha_0p_0=x_0-\alpha_0g_0\]
这里
\[\min_{\alpha_k}f(x_0-\alpha_0g_0)\]
选择最优步长$\alpha_0$，这里对于函数的泰勒展开
\[f(x_0-\alpha_0g_0)\approx f(x_0)-\alpha_0g_0^{\mathrm{T}}g_0+\frac{1}{2}\alpha_0^2g_0^{\mathrm{T}}Ag_0\]
这里二次函数Hessian矩阵是恒定的，$\triangledown^2f(x)\Bigl|_{x=x_k}=A$。或者直接代入二次函数推导一样，注意$g_0=Ax_0+d$即可。学习率的选择非常简单，就是直接对其进行求导
\[\frac{df(x_0-\alpha_0g_0)}{d\alpha_0}=-g_0^{\mathrm{T}}g_0+\alpha_0g_0^{\mathrm{T}}Ag_0=0\]
得到
\[\alpha_0=\frac{g_0^{\mathrm{T}}g_0}{g_0^{\mathrm{T}}Ag_0}\]
为满足上面定理的$g_{k+1}^{\mathrm{T}}p_i=0$，现在只需要$g_1^{\mathrm{T}}p_0=0$，由于
\[g_1=Ax_1+d=A\left ( x_0-\frac{g_0^{\mathrm{T}}g_0}{g_0^{\mathrm{T}}Ag_0}g_0 \right )+d\]
\[g_1^{\mathrm{T}}p_0=p_0^{\mathrm{T}}g_1=-g_0^{\mathrm{T}}A\left ( x_0-\frac{g_0^{\mathrm{T}}g_0}{g_0^{\mathrm{T}}Ag_0}g_0 \right )-g_0^{\mathrm{T}}d\]
\[=-g_0^{\mathrm{T}}\left ( Ax_0+d \right )+g_0^{\mathrm{T}}A\frac{g_0^{\mathrm{T}}g_0}{g_0^{\mathrm{T}}Ag_0}g_0\]
\[=-g_0^{\mathrm{T}}g_0+g_0^{\mathrm{T}}g_0=0\]
这样构造是满足的。那么如何进行选择新的下降方向$p_1$呢？选择
\[p_1=-g_1+\beta_0p_0\]
就是$g_0$，$g_1$的线性组合，且满足共轭条件
\[p_1^{\mathrm{T}}Ap_0=0\]
注意$g_{k+1}-g_k=A(x_{k+1}-x_k)=\alpha_kAp_k$，我们得到
\[-g_1^{\mathrm{T}}Ap_0+\beta_0p_0^{\mathrm{T}}Ap_0=0\]
\[\beta_0=\frac{g_1^{\mathrm{T}}Ap_0}{p_0^{\mathrm{T}}Ap_0}=\frac{g_1^{\mathrm{T}}\left ( g_1-g_0 \right )}{p_0^{\mathrm{T}}\left ( g_1-g_0 \right )}=\frac{g_1^{\mathrm{T}}g_1}{g_0^{\mathrm{T}}g_0}\left ( p_0=-g_0,g_1^{\mathrm{T}}p_0=0 \right )\]
因此我们得到$p_1$，以此类推，我们得到共轭梯度算法为
1 给定初始点$x_0$，梯度，$g_0=Ax_0+d$,下降方向$p_0=-g_0$
2 计算
\[\alpha_k=\frac{g_k^{\mathrm{T}}g_k}{g_k^{\mathrm{T}}Ag_k}\]
\[x_{k+1}=x_k+\alpha_kp_k\]
\[g_{k+1}=Ax_{k+1}+d\]
\[\beta_k=\frac{g_{k+1}^{\mathrm{T}}g_{k+1}}{g_k^{\mathrm{T}}g_k}\]
\[p_{k+1}=-g_{k+1}+\beta_kp_k\]
3 当$|g_k|<\epsilon$达到允许误差停止。
例子：比较沿直线最小化的最速下降法和共轭梯度法求解
$f(x)=\frac{1}{2}\begin{bmatrix} x_1 & x_2 \end{bmatrix}\begin{bmatrix} 2 & 0\\ 0 & 50 \end{bmatrix}\begin{bmatrix} x_1\\ x_2 \end{bmatrix}=x_1^2+25x_2^2$
图4.12 最速下降法（○）与共轭梯度法（△）的轨迹
最速下降法需要5步达到规定的精度，而共轭梯度法仅仅需要2步达到
\[|g_k|<\epsilon=10^{-5}\]
共轭梯度法程序：
clc; clear all;
maxk=50;%最大迭代次数
epsilon=1e-5;
x0=[0.5;0.5];
x(:,1)=x0;
A=[2 0;0 50];%f(x)=0.5*x*A*x+d'x+c
d=[0;0];
gk=A*x0+d;
pk=-gk;
for k=1:maxk
grad_old=gk;
alphak=-gk'*pk./(pk'*A*pk);
x0=x0+alphak*pk;
gk=A*x0+d;
betak=gk'*gk/(grad_old'*grad_old);
pk=-gk+betak*pk;
x(:,k+1)=x0;
if (norm(gk)<epsilon),
break,
end
end
x0
val=0.5*x0'*A*x0
dx = -0.1:0.01:0.55;
dy = dx;
[dX,dY] = meshgrid(dx);
dZ=dX.^2+25*dY.^2;
contour(dX,dY,dZ,100)
hold on
plot(x(1,:),x(2,:),'b-o')
思考：共轭梯度法中依然用到了二次函数的Hessian矩阵$A$，对于下面的非线性方程求解
$f(x)=2x_1^2+x_2^2-2x_1x_2-2x_2$
依然不合适，因为每一步都要求最优步长
\[\alpha_k=\frac{g_k^{\mathrm{T}}g_k}{g_k^{\mathrm{T}}Ag_k}\]
Hessian矩阵$A$对于二次函数是固定，但是一般非线性函数不是，是变化的，依然没有节省存储计算Hessian矩阵$A$，那么我们利用Armijo线搜索方法求解最优步长，可以解决这个问题，得到改进的非线性共轭算法。Armijo线搜索方法是，选择一个常数$\sigma \in (0,0.5)$,$\theta \in(0,1)$，步长为$\alpha_k=\theta^{m_k}$，整数$m_k$是满足下面不等式的最小正整数
\[f(x_k+\theta^{m_k}p_k)\leq f(x_k)+\sigma \theta^{m_k}g_k^{\mathrm{T}}p_k\]
对于下面的非线性方程
$f(x)=2x_1^2+x_2^2-2x_1x_2-2x_2$
精确解为$[1\quad 2]^{\mathrm{T}}$，最小值为-2。程序如下：
function [k,x,val]=frcg(fun,gfun,x0,epsilon,N)
% 功能: FR非线性共轭梯度法求解无约束问题: min f(x)
% 输入: fun, gfun分别是目标函数及其梯度, x0是初始点,
% epsilon是容许误差, N是最大迭代次数
if nargin<5, N=1000; end
if nargin<4, epsilon=1.e-5; end
beta=0.6; sigma=0.4;
n=length(x0); k=0;
while(k<N)
gk=feval(gfun,x0); %计算梯度
itern=k-(n+1)*floor(k/(n+1));
itern=itern+1; %计算搜索方向
if(itern==1)
dk=-gk;
else
betak=(gk'*gk)/(g0'*g0);
dk=-gk+betak*d0; gd=gk'*dk;
if(gd>=0.0), dk=-gk; end
end
if(norm(gk)<epsilon), break; end %检验终止条件
m=0; mk=0;
while (m<20) %用Armijo搜索求步长
if (feval(fun,x0+beta^m*dk)<=feval(fun,x0)+sigma*beta^m*gk'*dk)
mk=m; break;
end
m=m+1;
end
x=x0+beta^mk*dk;
g0=gk; d0=dk;
x0=x; k=k+1;
end
val=feval(fun,x);
function f=fun(x)
f=2*x(1)^2+x(2)^2-2*x(1)*x(2)-2*x(2);
end
function gf=gfun(x)
gf=[4*x(1)-2*x(2);2*x(2)-2*x(1)-2];
end
给定初值x0和误差epsilon，以及最大迭代次数N即可。
浙公网安备 33010602011771号