終於要正式進入求解拉!
在開始之前，我必須要介紹一篇非常適合入門的教學文章
An Introduction to the Conjugate Gradient Method Without the Agonizing Pain
對於iterative解\(Ax=b\)，希望讀者已經了解目標式的定義與normal equation
本文就其中Steepest Descent Method開始(其實也只是翻譯)
Introduction to Gradient Descent
\[f(x)=\left \|Ax-b\right \|^2=x^TA^TAx-2x^TA^Tb+b^Tb\]
目標是求\(x\)使得\(f(x)\)最小，一次偏微分求導數為0
\[f^{‘}{(x)}=Ax-b=0\]
residual \(r_{(i)}\)表第\(i\)次iteration的觀察值\(b\)與\(Ax_{(i)}\)的差，同時也表在\(x_{(i)}\)位置的梯度反方向，指向\(f(x)\)變小的位置
\[r_{(i)}=b-Ax_{(i)}=-f^{‘}{(x_{(i)})}\]
將\(x\)往梯度數值變小的方向\(r_{(i)}\)移動，步伐大小為\(\alpha\)。注意\(r_{(i)}\)非單位向量，\(\alpha\)是搭配\(r_{(i)}\)本身長度使用
\[x_{(i+1)}=x_{(i)}+{\alpha}r_{(i)}\]
上式就是Steepest Descent Method的迭代核心，不斷的更新\(x\)直到收斂為止；你可以隨便帶一個\(\alpha>0\)，就是一個標準的Gradient Descent解法，但我們其實可以更科學!
Optimized \(\alpha\) using Directional Derivative
現在出現了一個子問題，求最佳的\(\alpha\)，嘗試在\(r_{(i)}\)方向上進行line search找到\(f(x_{(i+1)})\)最小值；這是一個directional derivative問題，也是Steepest Descent Method中最困難的部分
\[\color{red}{\bigtriangledown_{r_{(i)}}f(x_{(i+1)})=0}\\f^{‘}(x_{(i+1)})\cdot{r_{(i)}}=0\\-r_{(i+1)}\cdot{r_{(i)}}=0\\r_{(i+1)}^{T}{r_{(i)}}=0\]
接下來想辦法把\(r_{(i+1)}\)中有關未知的\(x_{(i+1)}\)項消除，展開\(r_{(i+1)}\)並整理出\(\alpha\)
\[(b-Ax_{(i+1)})^{T}r_{(i)}=0 \\ (b-A(x_{(i)}+{\alpha}r_{(i)})^{T} r_{(i)}=0 \\ (b-(Ax_{(i)}+\alpha{Ar_{(i)}}))^{T}r_{(i)}=0 \\ ((b-Ax_{(i)})-\alpha{Ar_{(i)}})^{T}r_{(i)}=0 \\ (r_{(i)}-\alpha{Ar_{(i)}})^{T}r_{(i)}=0 \\ r^{T}_{(i)}r_{(i)}-\alpha{r^{T}_{(i)}}{A^T}r_{(i)}=0\]
\[\alpha=\frac{r^{T}_{(i)}r_{(i)}}{r^{T}_{(i)}A^{T}r_{(i)}}\]
離我們常見的\(\alpha\)還差一點，\(\alpha\)是常數所以轉置後相等
\[\alpha=\alpha^T=\frac{r^{T}_{(i)}r_{(i)}}{r^{T}_{(i)}Ar_{(i)}}\]
Optimize Computation Performance
整個迭代過程最大的overhead是在\(Ar_{(i)}\)上，其次就是兩個向量的內積部分，在large sparse \(A\)問題中通常會將這個矩陣乘法用其他方式取代，如FFT、spatial convolution。
再者就是每次要重新帶入\(x_{(i+1)}\)才能得到\(r_{(i+1)}\)也是個負擔，我們將residual公式修改，同乘上\(-A\)，再加上\(b\)
\[x_{(i+1)}=x_{(i)}+{\alpha}r_{(i)}\\-Ax_{(i+1)}=-Ax_{(i)}-{\alpha}Ar_{(i)}\\b-Ax_{(i+1)}=b-Ax_{(i)}-{\alpha}Ar_{(i)}\\r_{(i+1)}=r_{(i)}-{\alpha}Ar_{(i)}\]
我們得到\(r_{(i+1)}\)不用經過任何多餘的矩陣乘法，最複雜的\(Ar_{(i)}\)在計算\(\alpha\)時就會用到了。但在電腦上可能會存在浮點數的round-off error，隨著迭代次數的增加，這個誤差可能會不斷累加；或許可以考慮隔個幾次iteration再校正一次residual。
Steepest Descent Method
\(r_{(0)}=b-Ax_{(0)}\)
for \(i=0\) to convergence
\(\alpha=\frac{r^{T}_{(i)}r_{(i)}}{r^{T}_{(i)}Ar_{(i)}}\)
\(x_{(i+1)}=x_{(i)}+{\alpha}r_{(i)}\)
\(r_{(i+1)}=r_{(i)}-{\alpha}Ar_{(i)}\)
\(i=i+1\)
end for
Properties of Steepest Descent Method
這裡注意一個特性：\(r_{(i+1)}\)和\({r_{(i)}}\)互為orthogonal(如附圖綠色線\(r_{(i+1)}^{T}{r_{(i)}}=0\))，前後一次梯度方向是垂直的；我們發現了一個zigzag的低效移動方式，不斷的重複兩個方向，而每次更改方向就代表著額外的迭代計算，為什麼不要將同方向的移動量合併，一次走足呢? 於是乎conjugate gradient method(附圖紅色線)就出來啦~! 理論上在n維度空間中最多進行n次迭代就可以找到最佳解(更多時候遠小於n!)