當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

共轭梯度法的简单直观理解

發(fā)布時間：2023/12/20 编程问答 33 豆豆

生活随笔收集整理的這篇文章主要介紹了共轭梯度法的简单直观理解小編覺得挺不錯的,現(xiàn)在分享給大家,幫大家做個參考.

共軛梯度法的簡單直觀理解

參考資料
What is: 什么是共軛梯度法？(簡單直觀理解）
- 共軛向量
- 共軛方向
- 誤差與殘差
- 搜索方向的確定
- 步長，或者說系數(shù)alpha
How to: 怎么用共軛梯度法？(完整算法）
Python numpy代碼
Why: 為什么共軛梯度法能求解Ax=b?
- 二次型
- 將Ax=b問題轉(zhuǎn)化為最優(yōu)化問題
拓展：改進——預(yù)處理的共軛梯度法

參考資料

本文是參考以下內(nèi)容，結(jié)合自己的理解做的筆記。請盡量直接訪問下述網(wǎng)頁。

矩陣與數(shù)值計算（11）——共軛梯度法

共軛梯度法（一）：線性共軛梯度

無痛版共軛梯度法介紹(更新到第五章)

共軛梯度法通俗講義
第4個資料尤其清晰完備，很多都是參考它的。

What is: 什么是共軛梯度法？(簡單直觀理解）

共軛梯度法可以看作是梯度下降法（又稱最速下降法）的一個改進。

對梯度下降來說
$x?i+1=x?i?α?f\vec x_{i+1}=\vec x_i - \alpha\nabla f$
其中 $α\alpha$ 控制了一步要走多遠，因此被稱為步長，在機器學(xué)習里面又稱為學(xué)習率。

梯度下降法x移動的方向正是函數(shù)f的負梯度方向，這代表了局部上f減小最快的方向。

但是局部上減小最快的方向并不代表全局上指向最終解的方向。所以梯度下降法會出現(xiàn)像醉漢下山一樣走出zig-zag的路線。如下圖

圖1 梯度下降法在2維解空間（也就是解向量只有兩個維度）走出的路徑示意圖。
注：假如A是正定對稱陣，其2維解空間必定是橢圓的。

為什么會走出這一Z形線呢？因為梯度下降的方向恰好與f垂直，也就是說和等高線垂直。沿著垂直于等高線的方向，一定能讓函數(shù)減小，也就是最快地下了一個臺階。但是最快下臺階并不意味著最快到達目標位置（即最優(yōu)解），因為你最終的目標并不是直對著臺階的。

為了修正這一路線，采用另一個方向：即共軛向量的方向。

我們先暫且給出共軛梯度法最后的形式，方便字母的定義：
$x?i+1=x?i?αd?i\vec x_{i+1}=\vec x_i - \alpha \vec d_i$
對照梯度下降法，每次向下走的方向不是梯度了，而是專門的一個方向 $d?\vec d$ 。除此之外和梯度下降法幾乎一樣。

在推進下一步之前，我們來看看什么是向量共軛。

共軛向量

下面先簡要介紹共軛向量

所謂共軛向量，在數(shù)學(xué)上即：
$p_i^TAp_j=0$

其中A是一個對稱正定矩陣。
$p_i$ 和 $p_j$ 是一對共軛的向量。

可見，共軛是正交的推廣化，因為向量正交的定義為：
$p_i^Tp_j=0$
共軛比正交中間只多了個矩陣A，而矩陣的幾何意義正是對一個向量進行線性變換（可見Gilber Strang的線代公開課）。因此共軛向量的意思就是一個向量經(jīng)過線性變換（縮放剪切和旋轉(zhuǎn)）之后與另一個向量正交。

共軛方向

言歸正傳，如何找到一個更好的方向呢？我們首先可以看看最完美的方向是什么樣的。

下面這張圖展示的就是最完美的方向。圖中向量e代表的是誤差。向量d就是方向向量。

誤差e即當前迭代所得到的解與精確解的差值：
$e?i=x?i?x??\vec e_i=\vec x_i- \vec x^*$

可惜我們并不能找到誤差向量e，因為我們不知道精確解。

那么退而求其次，我們就找誤差向量的共軛向量。因為圖中可以看出，誤差向量是與方向向量垂直的，即正交。剛才說了，共軛就是正交的推廣。一個向量乘以一個矩陣之后與另一個方向正交，就是共軛。

即找到
$d?TAe?=0\vec d ^T A \vec e =0$

但是這個公式里面仍然含有e，我們必須想辦法去掉它，換成一個我們可以計算的量。

在推進下一步之前，我們先來看看誤差與殘差這兩個概念的區(qū)別。

誤差與殘差

前面寫道：

誤差error 即當前迭代所得到的解與精確解的差值：
$e?i=x?i?x??\vec e_i=\vec x_i- \vec x^*$

但是這種定義顯然是沒法直接用的，因為我們不知道精確解 $x^*$

那么退而求其次，我們想到，當誤差收斂為0的時候，必然滿足方程Ax=b，那么由此就可以定義出殘差residual：

$r?i=b??Ax?i\vec r_i=\vec b-A\vec x_i$

這個定義的精妙之處在于，它定義了Ax接近b的距離，當距離為0的時候，恰好就是精確解。但是又能避開精確解本身。

在實際的程序中，我們還常常定義相對殘差，即上一步迭代和這一步迭代的殘差的相對變化率，這里就不再贅述。

顯然，誤差和殘差之間就差了一個矩陣A，他們兩者的關(guān)系是這樣的：

$r?i=b??A(e?i+x??)=b??Ax???Ae?i=?Ae?i\vec r_i=\vec b - A(\vec e_i+\vec x^*)=\vec b - A \vec x^* -A\vec e_i = -A\vec e_i$

可見除了A，還多了個負號。

搜索方向的確定

言歸正傳，利用殘差，我們終于可以把誤差e給替換掉了：
于是前面的式子就變成了
$d?iTAe?i=?d?iTr?i=0\vec d_i ^T A \vec e_i =-\vec d_i ^T \vec r_i=0$

那么，這告訴我們：方向向量d，正是與殘差向量正交的方向！

接下來我們只需要構(gòu)建一個與殘差正交的向量就可以了。這部分內(nèi)容是由施密特正交化（更嚴謹一點的說法是共軛格萊姆-施密特過程）完成的。由于只是一個計算的方法，對概念的理解沒有幫助，所以我們跳過，直接給出結(jié)論。

每一步搜索方向的時候，這一方向與殘差以及前一步的方向有關(guān)
$d?i+1=r?i+1+βi+1d?i\vec d_{i+1} = \vec r_{i+1} +\beta_{i+1} \vec d_i$
其中系數(shù) $β\beta$ 可以這樣計算：
$βi+1=r?i+1Tr?i+1r?iTr?i\beta_{i+1} = \frac{ \vec r_{i+1}^T \vec r_{i+1} } {\vec r_{i}^T \vec r_{i} }$

這個系數(shù)beta其實很好記，因為分子就是殘差的內(nèi)積（下一步），分母也是殘差的內(nèi)積（這一步）。
或者說分子就是殘差長度的平方（下一步），分母也是殘差長度的平方（這一步）。（向量自己和自己的內(nèi)積就是它的長度）

從另一個角度額外補充一點理解：
每次走的方向恰好是與殘差正交的，這意味著：
每走一步恰好能消除殘差的一個方向！
所以，當消除了殘差所有投影方向上的值，那么就消除了整個殘差！

步長，或者說系數(shù)alpha

實際上，還有一點沒有解決，就是系數(shù) $α\alpha$ 怎么算？

這點的解釋我們以后有機會再說，直接給出結(jié)論。
$αi=r?i+1Tr?i+1d?iTAd?i\alpha_i = \frac{ \vec r_{i+1}^T \vec r_{i+1} } {\vec d_{i}^T A\vec d_{i} }$
這個alpha的分子和beta一樣，就是殘差的內(nèi)積。分母則是方向向量在乘以矩陣A之后的內(nèi)積。

How to: 怎么用共軛梯度法？(完整算法）

設(shè)定初值
$d?0=r?0=b??Ax?0\vec d_0=\vec r_0 = \vec b - A \vec x_0 \\$

計算系數(shù)alpha
$αi=r?i+1Tr?i+1d?iTAd?i\alpha_i = \frac{ \vec r_{i+1}^T \vec r_{i+1} } {\vec d_{i}^T A\vec d_{i} }$

迭代一步（向下走一步）
$x?i+1=x?i?αid?i\vec x_{i+1}=\vec x_i - \alpha_i \vec d_i$

計算殘差（此處已經(jīng)被修改，原文沒有被50整除那一個公式 2022-05-27）
如果迭代次數(shù)可以被50整除
$r?i+1=r?i?αiAx?\vec r_{i+1}=\vec r_i - \alpha_i A\vec x$
否則
$r?i+1=r?i?αiAd\vec r_{i+1}=\vec r_i - \alpha_i A d$

計算系數(shù)beta
$βi+1=r?i+1Tr?i+1r?iTr?i\beta_{i+1} = \frac{ \vec r_{i+1}^T \vec r_{i+1} } {\vec r_{i}^T \vec r_{i} }$

計算搜索方向 $d?\vec d$
$d?i+1=r?i+1+βi+1d?i\vec d_{i+1} = \vec r_{i+1} +\beta_{i+1} \vec d_i$

重復(fù)2~6，直到殘差足夠小

Python numpy代碼

import numpy as np import scipy.linalg as sl import matplotlib.pyplot as pltnn=4 #矩陣的規(guī)模 FIXME: 當規(guī)模>5的時候會出現(xiàn)震蕩，為什么？ accuracy = 1e-6 #使用共軛梯度法， A矩陣有兩個條件：1. 正定（特征值全為正數(shù)） 2.對稱 A = sl.pascal(nn, exact=False) # A是對稱正定矩陣, 10階帕斯卡矩陣, exact=False示用float元素而非默認的uint b = np.arange(1., nn+1., 1.) x0 = np.array([2.0]*nn) #x0def Conjugate_Gradient_Method(A,b,x0):#1. Set initial valuex = x0r = b - matrixVecProd(A, x)d = rrr = vecVecProd(r, r) #在計算beta的時候可以復(fù)用iter = 0relativeResidual = 0.1while( relativeResidual > accuracy or iter <1): #2. Compute alphaAd = matrixVecProd(A, d) # 在計算r_new和alpha的時候可以復(fù)用alpha = rr / vecVecProd(d, Ad)#3. step forwardx = x + alpha * d#4. compute residualif iter % 50 == 0 :r_new = b - matrixVecProd(A, x)else :r_new = r - alpha * Ad#5. compute betarr_new = vecVecProd(r_new, r_new) beta = rr_new / rrrr = rr_new#6. compute search directiond = r_new + beta * diter += 1 relativeResidual = np.linalg.norm(r_new) / np.linalg.norm(r)r = r_newprint("iter",iter, "relativeResidual",relativeResidual)return xdef matrixVecProd(A, vec):res = np.dot(A, vec)return resdef vecVecProd(vec1, vec2):res = np.dot(vec1, vec2)return res# ----------------TEST------------- def TEST_A(A,b):print(A)#A矩陣有兩個條件：1. 正定（特征值全為正數(shù)） 2.對稱eig = np.linalg.eig(A)print("eig=",eig[0]) #eig[0]取的是特征值print("b",b) # ----------------ENDTEST-------------def main():res = Conjugate_Gradient_Method(A, b, x0)print("-------------numerical result-------------------")print(res)print("-------------accurate result-------------------")accRes=np.dot(np.linalg.inv(A), b)print(accRes)if __name__ == "__main__":# TEST_A(A,b) # 可以先打印出來看看main()

這個代碼仍然是有問題的，主要是矩陣規(guī)模大的時候就會震蕩，我也不清楚為什么。
這里照抄一下劉天添課上的算法

Why: 為什么共軛梯度法能求解Ax=b?

說了這么多，其實有一個關(guān)鍵問題沒有講，那就是：為什么共軛梯度法能求解Ax=b?

按理說，共軛梯度法是函數(shù)最優(yōu)化的方法，怎么就扯上了求解Ax=b了呢？

實際上使用共軛梯度法的兩個條件

A是對稱的

A是正定的

也和這個原理有關(guān)。

數(shù)學(xué)家求解問題的思路是：把不會的問題轉(zhuǎn)化成會的問題，再套用會的問題的思路求解問題。

為了說明這一點，我們要從線性代數(shù)的二次型入手。我們可以先復(fù)習一下二次型，了解一下它是什么。

二次型

二次型就是關(guān)于向量的二次函數(shù)。

我們高中學(xué)過的二次函數(shù)通用表達式為
$f(x) = a x^2 +bx+c$

如果把其中的x替換為向量x，并且把a x^2 替換為
x^T A x 就得到了

$f(x) = x^T A x +bx+c$

這就是二次型。

二次型求導(dǎo)得到
$\frac{1}{2}( A x + A^T x)+b$

將Ax=b問題轉(zhuǎn)化為最優(yōu)化問題

我們本來求解的是
$Ax=bA\mathbf x=\bf b$

這個問題被轉(zhuǎn)化為了求某個函數(shù)的導(dǎo)數(shù)等于0的問題，即駐值問題。

我們設(shè)這個函數(shù)為g(x)。我們的問題即：

$g'(x)=0\\ x^*=argmin_x g(x)$
argmin_x的意思就是我們求取最小值的時候的x，而不是最小值本身。

這個 $x^*$ 就是最終解。

那么怎么聯(lián)系到Ax=b呢？

我們只要改造這個函數(shù)g，讓它的導(dǎo)數(shù)恰好就是Ax-b=0就好了！！
而這個函數(shù)，恰好就是二次型函數(shù)！

即
$g^{'} (x) = A x ? b$

于是求最小值得問題就能夠被轉(zhuǎn)化為求Ax=b的問題！

這里有個小小的瑕疵：
實際上，二次型g(x)的導(dǎo)數(shù)是
$g'(x)=1/2 (A^T+A)x-b$

所以我們就要限定 $A^T=A$ ，即限定A是對稱的。這就是第一個條件的由來！

to be continued
2022-05-20

拓展：改進——預(yù)處理的共軛梯度法

總結(jié)

以上是生活随笔為你收集整理的共轭梯度法的简单直观理解的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

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

上一篇：全国计算机vb考试经典程序设计,全国计算
下一篇：计算机网络基础专业找工作,2021计算机