3. 最速下降法

3.1 最速下降法的基本思想

最速下降法是典型的线搜索方法. 设 $f$ 是 ${\mathbb{R}}^n$ 上可微函数，在第 $k$ 次迭代时，若$\nabla f(x_k)=0$,则认为$x_k$ 已经是 $f$ 的极小点，否则令$$\begin{array}{r}{d}_k=-\nabla f(x_k).\end{array}$$则 $\nabla f(x_k{)}^T{d}_k=-\Vert \nabla f(x_k){\Vert }^2<0$.以 $-\nabla f(x_k)$ 为线搜索方向的方法称为最速下降法.

梯度方向是函数值增长最快的方向，那么梯度的负方向就是函数值下降最快的方向。也就是说梯度下降法选取了函数下降最快的方向来作为搜索方向.

需要说明的是，函数下降最快不代表函数值的下降幅度最大.下面证明这点.

设目标函数 $f$ 连续可微，将 $f$ 在 $x_k$ 处Taylor展开：
$$f(x)=f(x_k)+\nabla f(x_k{)}^T(x-x_k)+o(||x-x_k||)$$令 $x=x_{k+1}$ ,结合 $x_{k+1}=x_k+{\alpha }_k{d}_k$ 迭代法可得：
$$f(x_{k+1})=f(x_k)+{\alpha }_k\nabla f(x_k{)}^T{d}_k+o(||x_{k+1}-x_k||)$$如果迭代的点列 { x k } \{x_{k}\} {xk​} 是收敛的，那么很显然有 $x_{k+1}-x_k\to 0(k\to \infty )$ 。因此，我们有：$f(x_{k+1})=f(x_k)+{\alpha }_k\nabla f(x_k{)}^T{d}_k$ 。

为了使得本次迭代有意义，即 $f(x_{k+1})<f(x_k)$ ,需要确保 $f(x_{k+1})$ 随着 ${\alpha }_k$ 的增大而减小 ( $x_t$
的小领域内),即 $\frac{\partial f}{\partial \alpha }=\nabla f(x_k{)}^T{d}_k<0$ ,由于希望 $f(x_{t+1})$ 下降得尽可能快，也就是想要的满足：
$$d_k^{\ast }=\arg \underset{d_k}{\min }\frac{\partial f}{\partial \alpha }=\arg \underset{d_k}{\max }-\frac{\partial f}{\partial \alpha }$$由Cauchy不等式有：$$-\frac{\partial f}{\partial \alpha }=-\nabla f(x_k{)}^T{d}_k\le ||\nabla f(x_k)||\cdot ||d_k||$$根据Cauchy不等式成立的条件，当且仅当 $d_k$ 与 $-\nabla f(x_t)$ 同向时，上述不等式的等号成立。所以下降最快的方向 $d_k^{\ast }$ 与 $-\nabla f(x_k)$ 同向.也就是说梯度下降法只考虑了当前的函数关于 $alpha$ 下降最快的方向，但不一定是函数值下降幅度最大的.

3.2 基于精确搜索的最速下降法

当 $f(x_k+\alpha d_k)$ 关于 $\alpha$ 可微时，精确搜索算法中由于 ${\alpha }_k$ 是最优步长， 那么必然有 $\frac{d}{dt}f(x_k+t{d}_k){|}_{t={\alpha }_k}=0$ ，从而精确线搜索的迭代点列 $x_{k+1}=x_k+{\alpha }_k{d}_k$ 满足$\nabla f(x_{k+1}{)}^T{d}_k=0$, 即当 $\nabla f(x_k)\ne 0$ 时，有
$$\begin{array}{rl}\nabla f(x_{k+1}{)}^T{d}_k& =0,\text{其中}{d}_k:=-\nabla f(x_k).\end{array}$$考虑二次函数
$$f(x)=\frac{1}{2}{x}^{T}Ax+b^{T}x+c,$$其中 $A\in {\mathbb{S}}_{++}^n,b\in {\mathbb{R}}^n,c\in \mathbb{R}$.由于 $\nabla f(x_{k+1})=A{x}_{k+1}+b$. 代入(2)便得$$(A{x}_{k+1}+b{)}^T{d}_k=0.$$将 $x_{k+1}=x_k+{\alpha }_k{d}_k$ 代入上式，并利用 $\nabla f(x_k)=A{x}_k+b$ 整理便得
$$\begin{array}{r}{\alpha }_k:=\underset{\alpha >0}{\mathrm{argmin}}f(x_k+\alpha d_k)=-\frac{(A{x}_k+b{)}^T{d}_k}{d_k^{T}A{d}_k}=-\frac{\nabla f(x_k{)}^T{d}_k}{d_k^{T}A{d}_k}.\end{array}$$

下面看看用最速下降法求二次函数之收敛速度的估计. 记
$$\begin{array}{r}\Vert x{\Vert }_A^2:=x^{T}Ax,\forall x\in {\mathbb{R}}^n.\end{array}$$实际上，最速下降法求二次函数的收敛速度是可求的，具体地如下面的命题所示.

命题 3.2 (精确搜索最速下降法之收敛速度) 设目标函数为 $f(x)=\frac{1}{2}{x}^{T}Ax+b^{T}x+c$, 其中 $A$是一个$n$ 阶正定矩阵，$b\in {\mathbb{R}}^n,c\in \mathbb{R}$. 那么，任取初值 $x_0\in {\mathbb{R}}^n$, 使用精确搜索的最速下降法所得的点列 { x k } \{x_k\} {xk​} 当 $x_k\ne x^{\ast }$ 时，有$$\begin{array}{r}\frac{\Vert x_{k+1}-x^{\ast }{\Vert }_A}{\Vert x_k-x^{\ast }{\Vert }_A}\le \frac{\mathbf{cond}(A{)}_2-1}{\mathbf{cond}(A{)}_2+1},\end{array}$$其中，$x^{\ast }$为 $f$ 的极小值点，cond$(A{)}_2:={\lambda }_{\max }/{\lambda }_{\min }$ 是矩阵 $A$ 的条件数，${\lambda }_{\max },{\lambda }_{\min }$ 分别为 $A$ 的最大和最小特征值.

证.设 $x^{\ast }:={\mathrm{argmin}}_{x\in {\mathbb{R}}^n}f(x)$.则 $\nabla f(x^{\ast })=0$.由于 ${\nabla }^{2}f(x_k)=A$ ，且上文中我们规定了$\begin{array}{r}\Vert x{\Vert }_A^2:=x^{T}Ax,\forall x\in {\mathbb{R}}^n.\end{array}$ ，于是可以得到
$$\begin{array}{rl}f(x_k)& =f(x^{\ast })+\frac{1}{2}(x_k-x^{\ast }{)}^{T}A(x_k-x^{\ast })=f(x^{\ast })+\frac{1}{2}\Vert x_k-x^{\ast }{\Vert }_A^2.\end{array}$$将 $k$ 换成 $k+1$, 有$$\begin{array}{r}f(x_{k+1})=f(x^{\ast })+\frac{1}{2}\Vert x_{k+1}-x^{\ast }{\Vert }_A^2.\end{array}$$两式相减，得$$f(x_{k+1})-f(x_k)=\frac{1}{2}(\Vert x_{k+1}-x^{\ast }{\Vert }_A^2-\Vert x_k-x^{\ast }{\Vert }_A^2).$$所以$$\begin{array}{r}\frac{\Vert x_{k+1}-x^{\ast }{\Vert }_A^2}{\Vert x_k-x^{\ast }{\Vert }_A^2}=1+2\frac{f(x_{k+1})-f(x_k)}{\Vert x_k-x^{\ast }{\Vert }_A^2}.\end{array}$$记 $g_k:=\nabla f(x_k)$，可得$$\begin{array}{r}\begin{array}{rl}f(x_{k+1})-f(x_k)& =\nabla f(x_k{)}^T(x_{k+1}-x_k)+\frac{1}{2}(x_{k+1}-x_k{)}^{T}A(x_{k+1}-x_k)\\ & ={\alpha }_k{g}_k^T{d}_k+\frac{1}{2}{\alpha }_k^2{d}_k^{T}A{d}_k.\end{array}\end{array}$$

利用 $d_k=-g_k$ 和(3)有
$$\begin{array}{c}{\alpha }_k=-\frac{g_k^T{d}_k}{d_k^{T}A{d}_k}=\frac{g_k^T{g}_k}{g_k^{T}A{g}_k},{\alpha }_k{g}_k^T{d}_k=-\frac{(g_k^T{g}_k{)}^2}{g_k^{T}A{g}_k},\frac{1}{2}{\alpha }_k^2{d}_k^{T}A{d}_k=\frac{1}{2}\frac{(g_k^T{g}_k{)}^2}{g_k^{T}A{g}_k}.\end{array}$$代入(6)便得$$\begin{array}{r}f(x_{k+1})-f(x_k)=-\frac{1}{2}\frac{(g_k^T{g}_k{)}^2}{g_k^{T}A{g}_k}.\end{array}$$

因 $A(x^{\ast }+b)=\nabla f(x^{\ast })=0$, 我们有 $g_k=A{x}_k+b=A(x_k-x^{\ast })$.即 $x_k-x^{\ast }=A^{-1}{g}_k.$所以有
$$\begin{array}{rlr}\Vert x_k-x^{\ast }{\Vert }_A^2& =g_k^T(A^{-1}{)}^{T}A{A}^{-1}{g}_k& =g_k^T{A}^{-1}{g}_k.\end{array}$$
代入(5)便得
$$\begin{array}{c}\begin{array}{rl}\frac{\Vert x_{k+1}-x^{\ast }{\Vert }_A^2}{\Vert x_k-x^{\ast }{\Vert }_A^2}& =1-\frac{(g_k^T{g}_k{)}^2}{(g_k^{T}A{g}_k)(g_k^T{A}^{-1}{g}_k)}.\end{array}\end{array}$$利用 Kantorovich 不等式
( x T x ) 2 ( x T A x ) ( x T A − 1 x ) ≥ 4 λ max ⁡ λ min ⁡ ( λ max ⁡ + λ min ⁡ ) 2 , ∀ x ∈ R n ∖ { 0 } , \begin{aligned}\frac{(x^Tx)^2}{(x^TAx)(x^TA^{-1}x)}\geq\frac{4\lambda_{\max}\lambda_{\min}}{(\lambda_{\max}+\lambda_{\min})^2},\quad\forall x\in\mathbb{R}^n\setminus\{0\},\end{aligned} (xTAx)(xTA−1x)(xTx)2​≥(λmax​+λmin​)24λmax​λmin​​,∀x∈Rn∖{0},​

其中 ${\lambda }_{max},{\lambda }_{min}$分别为 $A$ 的最大和最小特征值，化简即得(4).

根据上面的证明可以发现利用最速下降法求二次函数的极小点的收敛阶至少是线性的.并且由于矩阵 $A$ 是对称正定矩阵，那么根据解析几何相关的知识可以知道 $f(x)$ 的等值面是一个椭球，而椭球的轴长就是 $A$ 的各特征值，并且其最长轴和最短轴之比恰为 $A$ 的谱条件数 $\mathbf{cond}(A{)}_2$，并且可以知道条件数的值越大，则 $f(x)$ 的等高面的长轴与短轴相差越大(也就是说椭球会越扁)， 迭代点列的收敛就越慢.当谱条件数 $\mathbf{cond}(A{)}_2$接近 1 时，等高线接近一个球，此时收敛很快.特别地，当谱条件数 $\mathbf{cond}(A{)}_2$ = 1 时 ，$\Vert x_1-x\ast \Vert A=0$， 也就是说只需要一次迭代就可以达到精确解.

3.3 基于精确搜索的最速下降法的程序实现

最优步长的推导

在此之前我们还需要讨论一下基于精确搜索的二次目标函数的最优步长应该如何选取.设讨论的二次型函数如下：$$f(x)=\frac{1}{2}{x}^{T}Qx+b^{T}x+c$$其中​ $x$ 是​ $n$ 维列向量， $A$ 是实对称正定矩阵， ​$b$ 是与​ $x$ 维度相同的列向量， $c\in \mathbb{R}$.

我们记 $g_k=\nabla f(x_k)$ 。由矩阵求导的性质以及 $A$ 是对称矩阵，那么有 $\nabla f(x_k)=A{x}_k+b$ 时，就可以开始在$d_k=-\nabla f(x_k)$ 的方向上实施精确线搜索，从而得到最优步长 ${\alpha }^{\ast }$ 了：$$\begin{array}{rl}{\alpha }^{\ast }& =\arg \underset{\alpha }{\min }f(x_k+\alpha d_k)\\ & =\arg \underset{\alpha }{\min }f(x_k-\alpha g_k)\\ & =\arg \underset{\alpha }{\min }\left(\frac{1}{2}{g}_k^{T}A{g}_k{\alpha }^2-g_k^T{g}_k\alpha +f(x_k)\right)\end{array}$$注意到上述函数是关于 $\alpha$ 的二次函数，故最优步长 ${\alpha }^{\ast }$ 为：
$${\alpha }^{\ast }=\frac{g_k^T{g}_k}{g_k^{T}A{g}_k}$$因此在对正定二次型做梯度下降时，通过上述上述上述的精确线搜索得到最优步长后，更新迭代式为$$x_{k+1}=x_k-\frac{g_k^T{g}_k}{g_k^{T}A{g}_k}(A{x}_k+b)$$下面做一个正定二次型的最速下降法实现的总结，已知迭代的终止条件参数$ϵ$。初始迭代点为$x_0$,记$g_k=\nabla f(x_k)$,重复以下操作 :

1.计算$d_k=-\nabla f(x_k)$,若$||d_t||<ϵ$,迭代终止，否则执行下一步

2.计算精确线搜索的最优步长：${\alpha }^{\ast }=\frac{g_k^T{g}_k}{g_k^{T}A{g}_k}$

3.通过迭代式 $x_{k+1}=x_k-\frac{g_k^T{g}_k}{g_k^{T}A{g}_k}(A{x}_k+b)$ 更新得到下一个迭代点，然后回退到第一步

程序实现

实际上，可以证明利用最速下降法来求二次函数的最优化问题时，如果存在局部最优点，那么该点是一定全局最优点，也就是说该优化问题满足全局收敛性.

下面基于此事实，利用Python来看看基于精确所搜的最速下降法的程序实现.

设需要优化的目标二次函数如下：$$f(x)=\frac{1}{2}{x}^{T}Qx+b^{T}x,Q=\left(\begin{array}{cc}10& -9\\ -9& 10\end{array}\right),b=(4,-15{)}^T$$为了检验迭代过程的可靠性，我们可以解析地求出该优化问题的最优点：$$x^{\ast }=-Q^{-1}b=(5,6{)}^T$$此外我们还需要定义目标函数中的 $Q,b$ 还有原本的函数 $f$ 和函数梯度 $\nabla f$ 的映射func和gradient.准备工作部分的代码如下：

# 准备工作
import numpy as np
import matplotlib.pyplot as pltQ = np.array([[10, -9], [-9, 10]], dtype="float32") # 矩阵Q
b = np.array([4, -15], dtype="float32").reshape([-1, 1]) # 向量b
# function and its gradient defined in the question
func = lambda x: 0.5 * np.dot(x.T, np.dot(Q, x)).squeeze() + np.dot(b.T, x).squeeze() # 利用lambda表达式定义目标函数：1/2 x^T*Q*x + b^Tx
gradient = lambda x: np.dot(Q, x) + b # 定义梯度 = f' = Qx + b
x_0 = np.array([5, 6]).reshape([-1, 1]) # 定义最优点

准备工作完成后，编写基于精确搜索的最速下降法的函数实现：

# GD algorithm
def gradient_descent(start_point, func, gradient, epsilon=0.01):""":param start_point: start point of GD 初始值:param func: map of plain function 目标函数:param gradient: gradient map of plain function 梯度函数:param epsilon: threshold to stop the iteration 迭代终止的阈值:return: converge point, # iterations"""assert isinstance(start_point, np.ndarray)  # assert that input start point is ndarrayglobal Q, b, x_0     # claim the global variencex_k_1, iter_num, loss = start_point, 0, []xs = [x_k_1] # xs是迭代点的序列，最后一个就是最后一次求的点while True:g_k = gradient(x_k_1).reshape([-1, 1])if np.sqrt(np.sum(g_k ** 2)) < epsilon:breakalpha_k = np.dot(g_k.T, g_k).squeeze() / (np.dot(g_k.T, np.dot(Q, g_k))).squeeze()x_k_2 = x_k_1 - alpha_k * g_kiter_num += 1xs.append(x_k_2)loss.append(float(np.fabs(func(x_k_2) - func(x_0))))if np.fabs(func(x_k_2) - func(x_k_1)) < epsilon:breakx_k_1 = x_k_2return xs, iter_num, loss

假设此处的迭代初值为​ $x_0=(4,4{)}^T$ ，进行梯度下降法，并可视化迭代过程：

x0 = np.array([4,4], dtype="float32").reshape([-1, 1])
xs, iter_num, loss = gradient_descent(start_point=x0, func=func,gradient=gradient,epsilon=1e-6)
print(xs[-1])	# last point of the sequence
plt.style.use("seaborn") # 样式
plt.figure(figsize=[12, 6]) # 大小
plt.plot(loss) 
plt.xlabel("# iteration", fontsize=12)
plt.ylabel("Loss: $|f(x_k) - f(x^*)|$", fontsize=12)
plt.yscale("log") # 对y轴取log
plt.show()

运行结果：

[[4.999424][5.998852]]

可以看到最终迭代点到达了最优点 $(5,6{)}^T$。当纵坐标取对数 $log$ 时，收敛速率近似是线性的。也就是说梯度下降法用于优化是可行的，收敛速率也不错。但是结论还不能够下得太早，我们在更换几组初始点，并把相关的参数画在一个图像内看看效果.

我们再取四个迭代点：​$(0,0),(0.4,0),(10,0),(11,0)$

# create the list of all starting point x_0
starting_points = [np.array([num, 0]).astype(np.float).reshape([-1, 1]) for num in [0.0, 0.4, 10.0, 11.0]]plt.figure(dpi=150)xss = []
# implement GD
for idx, start_point in enumerate(starting_points):xs, iter_num, losses = gradient_descent(start_point, func, gradient, epsilon=1e-6)target_point = xs[-1]xss.append(xs)# plot the losses of $|f(x_k) - f(x^*)|$plt.plot(np.arange(len(losses)), np.array(losses), label=f"start point: ({start_point[0][0]}, {start_point[1][0]})")loss = np.fabs(func(target_point) - func(x_0))print(f"{idx + 1}: start point:{np.round(start_point, 5).tolist()}, "f"point after GD:{np.round(target_point, 5).tolist()}, "f"loss:{np.round(loss, 16)}, # iterations: {iter_num}")print("-" * 60)plt.grid(True)
plt.legend()
plt.xlabel("# iteration", fontsize=12)
plt.ylabel("Loss: $|f(x_k) - f(x^*)|$", fontsize=12)
plt.yscale("log")
plt.title("Loss-iteration given to different starting points of GD", fontsize=18)
plt.show()

输出结果：

1: start point:[[0.0], [0.0]], point after GD:[[4.99855], [5.99826]], loss:2.9518850226e-06, # iterations: 60------------------------------------------------------------2: start point:[[0.4], [0.0]], point after GD:[[4.99902], [5.99873]], loss:1.6873629107e-06, # iterations: 42------------------------------------------------------------3: start point:[[10.0], [0.0]], point after GD:[[5.0], [6.0]], loss:1.33369e-11, # iterations: 4------------------------------------------------------------4: start point:[[11.0], [0.0]], point after GD:[[5.0], [6.0]], loss:0.0, # iterations: 1------------------------------------------------------------

由上图可以看到可以看到不同的初值点，收敛速率不相同。也就是说不同的初始点到达最优点的所需的迭代次数不一样.为了更好滴感受不同初始点选取的差异，可以把上述过程在二维平面内可视化出来：

plt.figure(dpi=150)
X = np.linspace(-2, 12, 200)
Y = np.linspace(-2, 12, 200)
XX, YY = np.meshgrid(X, Y)
Z = [func(np.array([XX[i, j], YY[i, j]], dtype="float32").reshape([-1, 1])).tolist() for i in range(200) for j in range(200)]
Z = np.array(Z).reshape([200, 200])
plt.contourf(XX, YY, Z, cmap=plt.cm.BuGn)plt.annotate(f"$(5.0, 6.0)$",xy=(5, 6),xytext=(5 - 2, 6 + 2),arrowprops={"color" : "black","shrink" : 0.1,"width" : 0.6})# plot the scatter
for idx, start_point in enumerate(starting_points):xx = [xss[idx][i][0] for i, _ in enumerate(xss[idx])]yy = [xss[idx][i][1] for i, _ in enumerate(xss[idx])]plt.plot(xx, yy, "o--", label=f"start point: ({start_point[0][0]}, {start_point[1][0]})")# add some tips for start pointplt.annotate(f"$({start_point[0][0]}, {start_point[1][0]})$",xy=(start_point[0][0], start_point[1][0]),xytext=(start_point[0][0] - 1.5, start_point[1][0] + idx + 2),arrowprops={"color" : "black","shrink" : 0.1,"width" : 0.6})plt.grid(True)
plt.title("Line Search For Two-dimensional Diagrams", fontsize=18)
plt.xlabel("$x_1$", fontsize=12)
plt.ylabel("$x_2$", fontsize=12)
plt.legend()
plt.show()

运行结果：

这张图可以直观地感受四个迭代点的迭代步数的差距。究其原因，我们在命题 3.2中提到，等高线代表的椭圆的长轴和短轴比的含义，不同初始点的选取所位于的等高面是不一样的，如果位于狭长的区域则收敛的速度会相应地快些，这是因为梯度方向的值要大些，从上图中可以直观的感受.

3.4 基于精确搜索的最速下降法的缺点

首先拿简单的二次型来说，随着条件数 $cond(A{)}_2$ 的增大，等高线的椭圆族会越来越扁，那么靠近椭圆长轴的迭代初值迭代到最优点会更加困难；或者说，精确线搜的梯度下降法的迭代次数会对迭代初值的选取越来越敏感。由于实际计算中初始值的选取是随机的，我们并不希望算法对初始值的选取过于敏感而导致算法的不稳定性.

其次最速下降法的下降方向选取存在一定的问题：虽然最速下降法选取了函数下降最快的方向，但是这并不一定就是函数值幅度减少最大的方向.简要地说这是因为，最快的下降方向只是对于一个很小的领域而言，但是在这个领域之外的函数值改变的幅度并没有考虑进去 。此外加上精确线搜的要求，使得整个算法变得很不灵活，因为程序会每步都找到 的最小值，但是由于每次迭代的函数值下降的幅度并不是最大的，这样就造成了一种计算资源的浪费.

关于下降最快不代表下降幅度最大可以如下证明：

设目标函数 $f$ 连续可微，将 $f$ 在 $x_k$ 处Taylor展开：
$$f(x)=f(x_k)+\nabla f(x_k{)}^T(x-x_k)+o(||x-x_k||)$$令 $x=x_{k+1}$ ,结合 $x_{k+1}=x_k+{\alpha }_k{d}_k$ 迭代法可得：
$$f(x_{k+1})=f(x_k)+{\alpha }_t\nabla f(x_k{)}^T{d}_t+o(||x_{k+1}-x_k||)$$如果迭代的序列点 { x k } \{x_{k}\} {xk​} 是收敛的，那么很显然有 $x_{k+1}-x_k\to 0(t\to \infty )$ 。因此，我们有：$f(x_{k+1})=f(x_k)+{\alpha }_k\nabla f(x_k{)}^T{d}_k$ 。

要是每次迭代都是有意义的，则要求 $f(x_{k+1})<f(x_k)$ ,需要确保 $f(x_{k+1})$ 随着 ${\alpha }_k$ 的增大而减小( $x_k$
的小领域内),即 $\frac{\partial f}{\partial \alpha }=\nabla f(x_k{)}^T{d}_k<0$ ，同时我们希望 $f(x_{k+1})$ 下降得尽可能快，也就是满足：$$d_t^{\ast }=\arg \underset{d_t}{\min }\frac{\partial f}{\partial \alpha }=\arg \underset{d_t}{\max }-\frac{\partial f}{\partial \alpha }$$由Cauchy不等式：
$$-\frac{\partial f}{\partial \alpha }=-\nabla f(x_k{)}^T{d}_k\le ||\nabla f(x_k)||\cdot ||d_k||$$根据Cauchy不等式成立的条件，当且仅当 $d_k$ 与 $-\nabla f(x_k)$ 同向时，上述不等式中的等号成立.也就是说下降最快的方向$d_k^{\ast }$与$-\nabla f(x_k)$ 同向。

Reference

本笔记内容的程序实现不分参考如下了文章：
最优化方法复习笔记（一）梯度下降法、精确线搜索与非精确线搜索（推导+程序）