共軛梯度法(1)--的基本原理

之前已經搞明白了，梯度下降法的基本原理，當然解釋的調度是從求函數極值的角度出發的，事實上從這個角度來理解，個人感覺是一個最為直接的理解角度，其完完全全是建立在多變量函數的微分系統中的。事實上這個方法(思想)在實際中是被用于求解線性方程的，當然單純的梯度下降形式并沒有被直接采用，被廣泛用于求解對稱、正系數矩陣方程組的方法是，基于梯度下降原理實現的共軛梯度法，這一方法在大型稀疏線性方程組的求解中，有極高的效率。但是，之前對于梯度下降思想的描述，好像認為這個方法是用于求函數極值的方法，似乎與線性方程組求解問題無關。因此，對于對稱、正定系數矩陣的二次型的理解是很必要的，它能讓我們明白矩陣系統中線性方程組與函數極值的關系，能夠解釋為什么一個求函數極值的思想可以被用于求解線性方程組系統。

$Ax = b$對應的二次型

這里就不嚴格去給出嚴格的數學表示了，因為對于非數學專業的我們來說，好像沒啥必要，知曉核心的思想更為重要。注意本次針對的系數矩陣$A$都是對稱、正定的，線性方程組系統$Ax = b$對應的是二次型對于自變量導數為0時的表達式，不嚴謹地可將其看作是函數吧，即二次型為

$$

f(x) = \frac{1}{2} x^{T} A x - b^{T} x + c \tag{1}

$$

對于$Ax = b$來說，$c$是一個0向量。我們知道求$f(x)$極小值的方法最直接的一種就是利用其導數為0，獲得極值點從而獲得極小值，因此，求導得到

$$

f'(x) = \frac{1}{2}A^T x + \frac{1}{2}A x - b = 0 \tag{2}

$$

因此，求解線性方程組$Ax = b$就等價于求解$f'(x) = 0$這個表達式，那么一階導數為0這個表達換一種說法，不就是去求二次型$f(x)$的函數的極值嘛！這也就解釋了為什么一個求解函數極值的方法卻可以應用到線性方程組的求解當中來。

為了更好的理解線性方程組$Ax = b$的求解與函數極值之間的關系，這里舉一個簡單的，可以被可視化的例子，若方程組為：

$$

\begin{bmatrix}

2 & 2\\

2 & 5

\end{bmatrix}

\begin{bmatrix}

x\\

y

\end{bmatrix}

=

\begin{bmatrix}

6\\

3

\end{bmatrix} \tag{3}

$$

這樣一個對稱正定的方程組，為了待求的未知數為了更符合函數表達，這里用了$x, y$來表示，依據二次型公式，得到

$$

f(x, y) = \frac{1}{2}

\begin{bmatrix}

x\\

y

\end{bmatrix}

\begin{bmatrix}

2 & 2\\

2 & 5

\end{bmatrix}

\begin{bmatrix}

x & y

\end{bmatrix} -

\begin{bmatrix}

6 & 3

\end{bmatrix}

\begin{bmatrix}

x\\

y

\end{bmatrix}

$$

這樣一個$2 \times 2$的線性方程組對應的二次型就是一個二元函數，為了得到表達式，這里利用sympy庫進行符號運算。import numpy as np

from sympy import *

# 先定義符號變量x, y

x, y = symbols('x y')

A = np.array([[2, 2], [2, 5]])

b = np.array([6, 3])

# 自變量向量這里用x1表示

x1 = np.array([x, y])

f_xy = 1/2 * x1.T @ A @ x1 - b.T @ x1

init_printing(use_unicode = True) # 更美觀的顯示

pprint(f_xy.expand())2 2

1.0?x + 2.0?x?y - 6?x + 2.5?y - 3?y

利用函數圖像來直觀地可視化函數極值情況。import matplotlib.pyplot as plt

import numpy as np

from mpl_toolkits.mplot3d import Axes3D

x = np.linspace(-10, 10, 100)

y = np.linspace(-10, 10, 100)

xx, yy = np.meshgrid(x, y)

f_xy = xx * (xx + yy) - 6 * xx + yy * (xx + 2.5 * yy) - 3 * yy

# 這里用等值線云圖與偽色彩圖來可視化，也順便對比下兩者的區別

fig, ax = plt.subplots(1, 2)

cf1 = ax[0].pcolormesh(xx, yy, f_xy, cmap='RdBu_r')

cs1 = ax[0].contour(x, y, f_xy, colors='gray', levels=15)

ax[0].clabel(cs1, inline=True, colors='k')

fig.colorbar(cf1, orientation='vertical')

ax[0].set_title('function by pcolormesh')

cf2 = ax[1].contourf(x, y, f_xy, cmap='RdBu_r')

cs2 = ax[1].contour(x, y, f_xy, colors='gray', levels=15)

ax[1].clabel(cs2, inline=True, colors='k')

ax[1].set_title('function by contourf')

fig.savefig('function2.jpg', dpi=600)

plt.show()

# 三維的函數圖像

fig1 = plt.figure()

ax = Axes3D(fig1)

cd = ax.plot_surface(xx, yy, f_xy, cmap='RdBu_r')

ax.set_title('the function of 3D image')

fig1.colorbar(cd)

fig1.savefig('image3d.jpg', dpi=600)

plt.show()

通過函數圖像的可視化，可以看到線性方程組本身是對稱正定的，那么對應的二次型的函數圖像則必定都位于在$z= 0 $平面以上，函數均是大于0的，且從函數的等值線及三維圖像中可以看到是存在極小值的。明白了這個關系之后，對于梯度下降法或者共軛梯度法應用于方程組(系數矩陣為對稱正定)的求解就不足為奇了。

Gradient descent method求解方程組

共軛梯度法的實現全來自于或者說是基于梯度下降法，因此，理解梯度下降法求解線性方程系統很有必要。在之前的梯度下降法的原理中，已經明白了這個方法對于函數極值的求解，是通過沿著梯度方向負向進行的，每一次的計算實質上獲得是自變量的增量向量，而現在面對方程組時，由于之前的二次型的內容已經得出了線性方程組與多元函數極值之間的關系，即$f'(x_{i}) = Ax_{i} - b = 0$。那么梯度下降法可以被這樣描述與構建：每一次產生的近似值$x_{i}$向量，近似解自然會有一個殘差向量(residual vector)為$r_{i}$， 下表$i$表示第幾步的迭代。

$$

r_{i} = b - A x_{i}\\

r_{i} = - f'(x_{i}) \tag{4}

$$

式(4)中的$r_{i}$不就是梯度向量嘛，表示了點下一步搜索或者移動的方向，當然要想確定自變量的增量大小，還需要知道步長，這一點在梯度下降原理中已經給出了解釋；在這里步長設定為$\alpha_{i}$，所以，產生的新的點$x_{i + 1}$為：

$$

x_{i + 1} = x_{i} + \alpha_{i} r_{i} \tag{5}

$$

到這里好像似乎與之前的原理中闡述并無大的區別，實際上現在問題在于如何確定步長$\alpha_i$，而在之前的介紹中步長是一個固定的，對于迭代計算來說，固定的步長的計算效率顯然是不能夠令人滿意的。對于函數的極小值的確定來說，在這里每一步獲得的新的函數值為$f(x_{i+1})$，對于整個求解過來說，只要這個更新的函數值達到最小值那么任務就算是完成了，依據式(5)，函數值可以看做是關于步長的函數，為了獲得極小值，通過導數為0來完成，即

$$

\frac{\mathrmozvdkddzhkzdf'(x_{i + 1})}{\mathrmozvdkddzhkzd \alpha_{i}} = 0 \tag{6}

$$

結合式(5)來看，這是一個復合函數，利用鏈式求導法則，得到

$$

式(6) = f'(x_{i+1}) ^{\mathrm{T}} \cdot \frac{\mathrmozvdkddzhkzdx_{i + 1}}{\mathrmozvdkddzhkzd\alpha_{i}} = 0 \tag{7}

$$

這個式子中，$f'(x)$不就是梯度向量嘛，依據式(4)給出的定義，在這里它也稱作殘差向量$-r_{i + 1}$，$x_{i + 1}$求導的結果可以根據式(5)為$r_{i}$，所以，要保證函數取得極值，就要使得前后兩步的梯度向量(也就是殘差向量$r$)保證正交！這里寫為內積形式：

$$

r_{i + 1} ^ {\mathrm{T}} \cdot r_{i} = 0 \tag{8}

$$

那么依據這個關系，帶入殘差向量的表達，得到：

$$

(b- A x_{i + 1}) ^ {\mathrm{T}} \cdot r_{i} = 0\\

(b - A (x + \alpha_{i} r_{i}))^{\mathrm{T}} \cdot r_{i} = 0

\tag{9}

$$

通過這個式子就可以得到步長$\alpha_{i}$的計算公式：

$$

\alpha_{i} = \frac{r_{i} ^{\mathrm{T}} r_{i}}{r_{i} ^{\mathrm{T}} A r_{i}} \tag{10}

$$

到這里，所有量的表達都得到了解決，那么對于線性方程組系統來說，梯度下降法的迭代過程為：

$$

repeat:\\

r_{i} = b - A x_{i}\\

\alpha_{i} = \frac{r_{i}^{\mathrm{T}} r_{i}}{r_{i} ^{\mathrm{T}} A r_{i}}\\

x_{i + 1} = x_{i} + \alpha_{i} r_{i}\\

until: x_{i}的增量幅度滿足設定精度\\

end

\tag{11}

$$

迭代計算的代碼如下：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]])

while np.linalg.norm(r0, np.inf) > 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, :] = x0

number += 1

if number > max_time:

print('warning:the result is not stable!')

return

return x, number, r

這里對一個$3\times3$大小的對稱正定矩陣進行驗算，梯度下降算法迭代了67次得到結果，$[0.99983945,0.99976565, 1.99978575]$，而這個方程組的精確解為$[1,1,2]$。可以迭代結果還是比較可信與可靠的。A = np.array([[4, -2, -1], [-2, 4, -2], [-1, -2, 3]])

b = np.array([0, -2, 3])

x0 = np.array([1, 1, 1])

x_gd = gd_method(A,b, x0, 0.0001,500)

print(x_gd[0][x_gd[1]-1,:])

print(x_gd[1])

print(x_gd[2])[0.99983945 0.99976565 1.99978575]

67

[[-1. -2. 3. ]

[-0.39130435 0.43478261 0.15942029]

[ 0.09201888 0.02436289 0.15942029]

...

[ 0. 0. 0. ]

[ 0. 0. 0. ]

[ 0. 0. 0. ]]x = np.linalg.solve(A,b)

print(x)[1. 1. 2.]

依據上述的分析以及這個算例的計算，實際上對于線性方程組迭代的迭代計算會有一個新的認識。首先，從最為原始的迭代的角度來看，通過每一迭代步中都會產生一個殘差向量，似乎在方程組的求解角度來看，殘差向量僅僅也就表示一個誤差吧，好像也并沒有很直觀的一個認識。但是從梯度下降的算法角度來看，線性方程組迭代步中的殘差向量實際上就是二次型(函數)的負梯度向量，那么也就是說，這個向量$r$不斷地在修正點移動的方向，使得移動方向不斷的向函數的極小值點處靠近！由于gd_method函數中會記錄每一步的參差向量，因為這個是一個3個未知數的方程組，因此殘差向量$r$本身就是一個$R^{3}$空間中的向量，同時返回得到的$x$矩陣中存放著每一步的近似解，這樣進行可視化，可以看出梯度下降中前后的參差向量確實是正交的。因為三元函數無法進行可視化，不是很容易看出梯度下降的搜尋過程。可能二元函數的化對于理解更方便吧。xs = x_gd[0][0:x_gd[1]-1,:]

plt.style.use('ggplot')

fig = plt.figure()

ax = Axes3D(fig)

ax.plot(xs[:,0], xs[:,1], xs[:,2], marker='o')

ax.set_title('point in every step')

ax.set_xlabel('x')

ax.set_ylabel('y')

ax.set_zlabel('z')

fig.savefig('gdimage.jpg', dpi=600)

plt.show()

在上述的分析中，我們會發現實際好像似乎梯度下降這個迭代算法的效率并不是很高，很簡單的方程組都要迭代67次，顯然這樣的效率去求解大型的稀疏方程組是不能令人滿意的。這個方法之所以沒有被廣泛地應用，主要原因是因為每一步的殘差向量也就是梯度向量，都必須與相鄰步的殘差向量保持正交，如果為了便于理解，以二維空間為例，就是每一次的搜索(移動)方向都是相互垂直的，那么間隔的殘差向量必定平行，也就說從初始值向精確解移動是呈折線的方式，一個搜索方向會被重復了好多次，這樣使得迭代效率并不高。但是梯度下降法實現的思想確是很有意義的，為函數極值的確定與線性方程組的求解之間建立了聯系。真正有意義的方法是在這個方法基礎上實現的共軛梯度法，它對搜索方向進行了共軛處理,在理解了梯度下降的原理后，再對共軛梯度法進行學習可能更好一點。

本站文章如非特別說明，均為原創。未經本人同意，請勿轉載！轉載請務必注明出處！

總結

以上是生活随笔為你收集整理的python共轭梯度法_Numerical Analysis: 共轭梯度法(1)--基本原理的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。