梯度下降法的缺点

之前的内容已经完全学习了梯度下降算法，对于这个方法之所以在线性方程组求解领域没有更广泛地流行起来，最大的原因是每一步的迭代中近似值移动的方向是沿着负梯度方向，而相邻迭代步的残差向量(本质就是负梯度向量)$r$之间是相互正交的，这个现象会导致可能在同一个方向上搜索（移动）了多次，造成迭代效率过低，当然迭代效率好与坏还受到初始值选择的影响，但是这个现象仍然不能够被忽略。为了说明这个现象，也为了后续方法的引入，这里还是用一个可视化的方式来反应这一现象，这也有助于对于后续共轭方向的理解。

这里采用的这个$2 \times 2$系数矩阵的例子在前一部分的内容中已经出现过，之所以对这个二维解向量空间的方程组进行求解与分析，主要原因在于它的解向量是一个二维空间中的向量，对于可视化较为方便与明晰。方程组如下，由于这里不在牵扯到对应的二次型，因此对于未知量设定为$x_1, x_2$:

$$ \begin{bmatrix} 2 & 2\\ 2 & 5 \end{bmatrix} \begin{bmatrix} x_1\\ x_2 \end{bmatrix} = \begin{bmatrix} 6\\ 3 \end{bmatrix} \tag{1} $$

梯度下降法的实现函数如下，之前的内容已经给出。

import numpy as np

# 每一步的残差向量r0都会被以行向量的形式存放在r矩阵中
def gd_method(A, b, x0, tol, max_time):
    r0 = b - np.dot(A, x0)
    number = 0
    r = np.zeros([max_time, b.shape[0]])
    x = np.zeros([max_time, b.shape[0]])
    x[0:, :] = x0
    r[0, :] = r0
    while np.linalg.norm(r0, 2) > tol:
        r0 = b - np.dot(A, x0)
        r[number, :] = r0
        alpha = np.dot(r0.T, r0) / np.dot(np.dot(r0,A), r0)
        x0 = x0 + alpha * r0
        x[number + 1, :] = x0
        number += 1
        if number > max_time:
            print('warning:the result is not stable!')
            return
    return x, number, r

求解如下，这里初始值选择为$(0, 0)$，迭代了39次达到精度要求。这里函数返回的对象是一个元组，需要进行解包，第一个元素为x矩阵，存放着每一步的解向量，以行进行存储；第二个为迭代次数，第三为每一步中的残差向量，也是按照行进行存储。迭代39次后获得了一个可以被接受的近似解$(3.99999669, -0.99999779)$。为了比较好的理解线性方程组，我们依旧将其对应的二次型函数的等值线图画出来，为了直观地感受梯度，同时提供了这个函数的在每一个采样点处的斜率方向，即斜率场。

A1 = np.array([[2, 2], [2, 5]])
b1 = np.array([6, 3])
x01 = np.array([0, 0])
x_gd1 = gd_method(A1,b1, x01, 0.00001,500)
print(x_gd1[1])
print(x_gd1[0][0:x_gd1[1] + 1])

39
[[ 0.          0.        ]
 [ 1.42857143  0.71428571]
 [ 2.04081633 -0.51020408]
 [ 2.74052478 -0.16034985]
 [ 3.04039983 -0.76009996]
 [ 3.38311418 -0.58874279]
 [ 3.52999176 -0.88249794]
 [ 3.69785184 -0.7985679 ]
 [ 3.76979188 -0.94244797]
 [ 3.85200907 -0.90133938]
 [ 3.887245   -0.97181125]
 [ 3.92751464 -0.95167643]
 [ 3.94477306 -0.98619327]
 [ 3.96449697 -0.97633131]
 [ 3.97295007 -0.99323752]
 [ 3.98261076 -0.98840717]
 [ 3.98675106 -0.99668776]
 [ 3.99148282 -0.99432188]
 [ 3.99351072 -0.99837768]
 [ 3.99582832 -0.99721888]
 [ 3.99682158 -0.99920539]
 [ 3.99795673 -0.99863782]
 [ 3.99844322 -0.99961081]
 [ 3.99899921 -0.99933281]
 [ 3.9992375  -0.99980937]
 [ 3.99950982 -0.99967321]
 [ 3.99962653 -0.99990663]
 [ 3.99975991 -0.99983994]
 [ 3.99981708 -0.99995427]
 [ 3.99988241 -0.9999216 ]
 [ 3.9999104  -0.9999776 ]
 [ 3.9999424  -0.9999616 ]
 [ 3.99995612 -0.99998903]
 [ 3.99997179 -0.99998119]
 [ 3.99997851 -0.99999463]
 [ 3.99998618 -0.99999079]
 [ 3.99998947 -0.99999737]
 [ 3.99999323 -0.99999549]
 [ 3.99999484 -0.99999871]
 [ 3.99999669 -0.99999779]]

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0, 5, 100)
y = np.linspace(-2, 2, 100)
xx, yy = np.meshgrid(x, y)
f_xy = xx * (xx + yy) - 6 * xx + yy * (xx + 2.5 * yy) - 3 * yy
g_x, g_y = np.gradient(f_xy)
g_xy = np.hypot(g_x, g_y)

fig, ax = plt.subplots()
cb = ax.pcolormesh(xx, yy, f_xy, cmap='Spectral_r')
ax.contour(xx, yy, f_xy, colors='gray')
ax.quiver(xx[::10, ::10], yy[::10, ::10], 
          g_x[::10, ::10] / g_xy[::10, ::10],
          g_y[::10, ::10] / g_xy[::10, ::10])
ax.set_xlabel('$x_{1}$')
ax.set_ylabel('$x_{2}$')
ax.set_title('the function of Ax=b')
fig.colorbar(cb, ax=ax)
fig.savefig('1.jpg', dpi=600)
plt.show()

梯度下降法中，相邻残差向量是相互正交的，从图中可以看出，每一步的向量都与前后步的向量正交，在二维平面上正交就意味着垂直，那么相隔的残差向量必定是平行的，即共线。梯度下降法的迭代或者说是搜索效率低，因为在某一个方向进行了多次的重复搜索，使得一个很简单的二元方程组的计算都需要数十次迭代，虽然这个迭代次数会受到初值的影响，但是迭代效率低是一个不争的事实。根据这个图中搜索步的显示，实际上很法线不论搜索了多少步，会发现这些向量只有两个方向，因为相隔的都是共线的，即本质上只存在两个方向，只不过没一个方向上进行了很多次的搜索。试想一下，如果对这些向量移动下，让它们首尾相接，那么似乎好像可以只进行两步就可以达到最终的点，这样岂不是效率更高。这就是共轭方向法的构想。

fig, ax = plt.subplots()
cs = ax.contour(xx, yy, f_xy, colors='k')
ax.plot(x_gd1[0][0:x_gd1[1]][:, 0], 
        x_gd1[0][0:x_gd1[1]][:, 1], 
        color='tab:red', marker='o')
ax.clabel(cs, inline=True, colors='k', fontsize=14)
ax.set_xlabel('$x_{1}$')
ax.set_ylabel('$x_{2}$')
fig.colorbar(cs, ax=ax)
fig.savefig('2.jpg', dpi=600)
plt.show()

一个改进的方向--共轭方向

依据之前的阐述，我们已经明确了一个思路就是，如果是一个$n \times n$的方程组，那么其解是一个$n$维向量，之前的例子是一个二维的，这里推广到更一般的情况。一个好的想法是找到n个“方向向量”，每一次的移动搜索都能完全的搞定一个未知数$x_{i}$，这样效率应该是最高的，那么只要$n$次后，所有的$x_{i}$都将被确定！这绝对是一个好的想法，当然按照这个思路来看，点移动的方向似乎不再是负梯度方向，而是方向向量$d_{i}$了，那么自然得到，下一个解可以表示为：

$$ x_{i + 1} = x_{i} + \alpha_{i} d_{i} \tag{2} $$

对于这个$2 \times 2$的方程组而言，初始值$x_0$向量沿着方向向量$d_0$移动，当然还需要乘以一个步长才行，获得第一步迭代结果$x_1$向量，依据前述的原理，解向量$x$中的第一个元素也就是二维坐标轴的第一个轴的值得到确定，此时的向量$x_1$与精确解$x$之间的距离，我们称之为误差向量，$e_1 = x_1 - x$，然后从这个点按照$d_1$移动，直接得到精确解，这样一共两步就可以完成计算。

按照上述的过程，可以确定一个关系，方向向量$d_i$与误差向量$e_{i + 1}$保持正交！

$$ d_{i}^{\mathrm{T}} \cdot e_{i + 1} = 0\\ e_{i + 1} = x_{i + 1} - x\\ d_{i}^{\mathrm{T}} \cdot (e_{i} + \alpha_{i} d_{i}) = 0\\ \alpha_{i} = - \frac{d_{i}^{\mathrm{T}} e_{i}}{d_{i}^{\mathrm{T}} d_{i}} \tag{3} $$

由于精确解本身是未知的，因此，误差向量其实也是未知的，因此，只是依据这个式子，还无法确定出步长,这里采用另一种正交方式。

$$ \alpha_{i} = \frac{r_{i}^{\mathrm{T}} r_{i}}{d_{i}^{\mathrm{T}} A d_{i}} \tag{4} $$

这个式子的推导利用了，矩阵与向量之间的正交！并且，依据误差向量与残差向量的相互关系$r_{i} = - A e_{i}$得出。

共轭梯度法

因此现在只要能找到这样的一组相互正交的方向向量$d_{i}$，那么实际上这个问题似乎就可以得到解决！而这一过程通过Gram-Schmidt过程来完成，而梯度下降法中残差向量$r_{i}$中相邻的向量互为正交，那么可以考虑使用残差作为这一组基，在此基础上进行共轭。用残差向量产生共轭后方向向量，在$n$步之后，即可得到结果。

$$ r_{i + 1} = - A e_{i + 1}\\ = -A (e_{i} + \alpha_{i} d_{i})\\ = r_{i} - \alpha_{i} A d_{i} \tag{5} $$

从残差向量构造方向向量，即

$$ d_{i} = r_{i} + \sum_{k=0} ^{i - 1} \beta_{ik}d_{k} \quad (i > k) \tag{6} $$

两侧构造矩阵向量正交,可以得到$\beta_{ij}$

$$ \beta_{ij} = -\frac{r_{i}^{\mathrm{T}} A d_{j}}{d_{j}^{\mathrm{T}} A d_{j}} \tag{7} $$

因为式(5)成立，可以继续进行化简，依据矩阵向量的正交性，$j$实际上只能取固定值，就可以得到关于$\beta_{ij}$的一个公式，这样的完全可以只写为$\beta_i$.

$$ \beta_{i} = \frac{r_{i + 1}^{\mathrm{T}} r_{i + 1}}{r_{i}^{\mathrm{T}} r_{i}} \tag{8} $$

这样，共轭梯度法的涉及到的计算过程及公式如下：

$$ d_{0} = r_{0} = b - Ax_{0}\\ \alpha_{i} = \frac{r_{i}^{\mathrm{T}} r_{i}}{d_{i}^{\mathrm{T}} A d_{i}}\\ x_{i + 1} = x_{i} + \alpha_{i} d_{i}\\ r_{i + 1} = r_{i} - \alpha_{i} A d_{i}\\ \beta_{i} = \frac{r_{i + 1}^{\mathrm{T}} r_{i + 1}}{r_{i}^{\mathrm{T}} r_{i}}\\ d_{i + 1} = r_{i + 1} + \beta_{i + 1} d_{i} \tag{9} $$

以上就是共轭梯度的全部过程和原理，以及实现的思想，如果你对这个过程完全了解，那么一定对这个方法很惊叹！因为，从迭代法的角度来看，它似乎有着超高的效率，有几个未知数就只需要迭代几次，极大地改进了梯度下降法的缺陷，对于一个迭代法而言，传统的迭代法都是获得的近似值，但是CG法确实通过有限步的迭代获得是精确解，最为重要的是它的迭代次数是完全确定的，这与传统的迭代法有着很大的区别，从这一个角度来看，CG法应该是一种直接法。从直接法的角度来看，它又不是完全意义上的直接法，因为该方法中需要给定的一个初始值才能求解，这一点又与迭代法有点相似。对于大型稀疏的线性方程系统来说，CG法是最好的方法！

def cg_method(A, b, x0):
    r0 = b - np.dot(A,x0)
    d0 = r0
    x = np.zeros([b.shape[0] + 1, b.shape[0]])
    x[0,:] = x0
    for i in np.arange(b.shape[0]):
        r0 = b - np.dot(A,x0)
        alpha = np.dot(r0.T, r0) / np.dot(np.dot(d0.T, A), d0)
        x0 = x0 + alpha * d0
        x[i + 1,:] = x0
        r1 = r0 - alpha * np.dot(A, d0)
        beta = np.dot(r1, r1) / np.dot(r0, r0)
        d0 = r1 + beta * d0
    return x

xl=cg_method(A1,b1, x01)

fig, ax = plt.subplots()
cs = ax.contour(xx, yy, f_xy, colors='k')
ax.plot(xl[:, 0], xl[:, 1], marker='*')
ax.clabel(cs, inline=True, colors='k', fontsize=14)
ax.set_xlabel('$x_{1}$')
ax.set_ylabel('$x_{2}$')
fig.colorbar(cs, ax=ax)
fig.savefig('3.jpg', dpi=600)
plt.show()

本站文章如非特别说明，均为原创。未经本人同意，请勿转载！转载请务必注明出处！