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

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

$Ax = b$對應(yīng)的二次型

這里就不嚴(yán)格去給出嚴(yán)格的數(shù)學(xué)表示了，因?yàn)閷τ诜菙?shù)學(xué)專業(yè)的我們來說，好像沒啥必要，知曉核心的思想更為重要。注意本次針對的系數(shù)矩陣$A$都是對稱、正定的，線性方程組系統(tǒng)$Ax = b$對應(yīng)的是二次型對于自變量導(dǎo)數(shù)為0時(shí)的表達(dá)式，不嚴(yán)謹(jǐn)?shù)乜蓪⑵淇醋魇呛瘮?shù)吧，即二次型為

$$

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

$$

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

$$

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

$$

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

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

$$

\begin{bmatrix}

2 & 2\\

2 & 5

\end{bmatrix}

\begin{bmatrix}

x\\

y

\end{bmatrix}

=

\begin{bmatrix}

6\\

3

\end{bmatrix} \tag{3}

$$

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

$$

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}

$$

這樣一個(gè)$2 \times 2$的線性方程組對應(yīng)的二次型就是一個(gè)二元函數(shù)，為了得到表達(dá)式，這里利用sympy庫進(jìn)行符號運(yùn)算。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

利用函數(shù)圖像來直觀地可視化函數(shù)極值情況。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

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

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()

# 三維的函數(shù)圖像

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()

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

Gradient descent method求解方程組

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

$$

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

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

$$

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

$$

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

$$

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

$$

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

$$

結(jié)合式(5)來看，這是一個(gè)復(fù)合函數(shù)，利用鏈?zhǔn)角髮?dǎo)法則，得到

$$

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

$$

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

$$

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

$$

那么依據(jù)這個(gè)關(guān)系，帶入殘差向量的表達(dá)，得到：

$$

(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}

$$

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

$$

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

$$

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

$$

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}的增量幅度滿足設(shè)定精度\\

end

\tag{11}

$$

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

這里對一個(gè)$3\times3$大小的對稱正定矩陣進(jìn)行驗(yàn)算，梯度下降算法迭代了67次得到結(jié)果，$[0.99983945,0.99976565, 1.99978575]$，而這個(gè)方程組的精確解為$[1,1,2]$。可以迭代結(jié)果還是比較可信與可靠的。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.]

依據(jù)上述的分析以及這個(gè)算例的計(jì)算，實(shí)際上對于線性方程組迭代的迭代計(jì)算會有一個(gè)新的認(rèn)識。首先，從最為原始的迭代的角度來看，通過每一迭代步中都會產(chǎn)生一個(gè)殘差向量，似乎在方程組的求解角度來看，殘差向量僅僅也就表示一個(gè)誤差吧，好像也并沒有很直觀的一個(gè)認(rèn)識。但是從梯度下降的算法角度來看，線性方程組迭代步中的殘差向量實(shí)際上就是二次型(函數(shù))的負(fù)梯度向量，那么也就是說，這個(gè)向量$r$不斷地在修正點(diǎn)移動(dòng)的方向，使得移動(dòng)方向不斷的向函數(shù)的極小值點(diǎn)處靠近！由于gd_method函數(shù)中會記錄每一步的參差向量，因?yàn)檫@個(gè)是一個(gè)3個(gè)未知數(shù)的方程組，因此殘差向量$r$本身就是一個(gè)$R^{3}$空間中的向量，同時(shí)返回得到的$x$矩陣中存放著每一步的近似解，這樣進(jìn)行可視化，可以看出梯度下降中前后的參差向量確實(shí)是正交的。因?yàn)槿瘮?shù)無法進(jìn)行可視化，不是很容易看出梯度下降的搜尋過程。可能二元函數(shù)的化對于理解更方便吧。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()

在上述的分析中，我們會發(fā)現(xiàn)實(shí)際好像似乎梯度下降這個(gè)迭代算法的效率并不是很高，很簡單的方程組都要迭代67次，顯然這樣的效率去求解大型的稀疏方程組是不能令人滿意的。這個(gè)方法之所以沒有被廣泛地應(yīng)用，主要原因是因?yàn)槊恳徊降臍埐钕蛄恳簿褪翘荻认蛄?#xff0c;都必須與相鄰步的殘差向量保持正交，如果為了便于理解，以二維空間為例，就是每一次的搜索(移動(dòng))方向都是相互垂直的，那么間隔的殘差向量必定平行，也就說從初始值向精確解移動(dòng)是呈折線的方式，一個(gè)搜索方向會被重復(fù)了好多次，這樣使得迭代效率并不高。但是梯度下降法實(shí)現(xiàn)的思想確是很有意義的，為函數(shù)極值的確定與線性方程組的求解之間建立了聯(lián)系。真正有意義的方法是在這個(gè)方法基礎(chǔ)上實(shí)現(xiàn)的共軛梯度法，它對搜索方向進(jìn)行了共軛處理,在理解了梯度下降的原理后，再對共軛梯度法進(jìn)行學(xué)習(xí)可能更好一點(diǎn)。

本站文章如非特別說明，均為原創(chuàng)。未經(jīng)本人同意，請勿轉(zhuǎn)載！轉(zhuǎn)載請務(wù)必注明出處！

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。