題外話：
從開始學習Slam十四講第六章的時候就開始想寫一個文檔整理一下這些年遇到的優化算法，一周學一章，現在都學到第9章了，總算半整理半引用整理出來了...
如果學一個東西是不斷坑自己+自己去填坑的過程，下一次應該不會摔的那么疼了吧😀
對於一個最小二乘問題的求解，根據目標函數可分為線性最小二乘和非線性最小二乘；
對於非線性最小二乘問題，通常是進行泰勒展開將問題線性化，求解線性增量方程或是直接迭代找到最優值；
對於線性最小二乘問題，通常是直接進行展開、求導等於零，構造\(A\vec{x}=\vec{b}\)的解方程問題，使用直接分解法或是迭代法求解；
寫完后發現文檔較長，還是列一下有些什么東西吧：
- 梯度下降與其擴展算法（隨機梯度下降、mini-batch梯度下降以及批梯度下降）
- 牛頓法與其優化算法（擬牛頓法、BFGS、LBFGS、高斯牛頓法以及列文伯格-馬誇爾特法)
- 求解線性最小二乘問題的那些：1）直接分解（LU、LUP、Cholesky分解求解方陣線性方程組問題，QR分解解決欠定方程組問題以及超定方程組的最小二乘解）；2）迭代法（雅各比迭代、高斯賽德爾迭代、SOR以及超級好用的共軛梯度）
- 一些自己覺得不錯的博客介紹；
非線性最小二乘問題
對於非線性最小二乘問題，通常會將目標函數進行泰勒展開，並將問題轉換為一個線性求解問題：
設有一個最小二乘問題：
有\(\vec{x}\in {R^n}, f\)是非線性函數，求解這個問題的常規思路是：
- 給定某個初始值\(\vec{x}_0\)
- 對於第k次迭代，尋找一個增量\(\Delta\vec{x}_k\)，使得\(||f(\vec{x}_k+\Delta\vec{x}_k)||_2^2\)
- \(\Delta\vec{x}_k\)足夠小，則停止
- 否則，令\(\vec{x}_{k+1}=\vec{x}_k + \Delta\vec{x}_k\)，返回第2步
將非線性最小二乘問題求解的目標：從尋找最優值轉換為尋找最小的\(\Delta\vec{x}_k\)，當函數下降到\(\Delta\vec{x}_k\)很小的時候，則等價為已找到最優值。
二階泰勒展開：
一階收斂
——僅考慮一階泰勒展開項，\(\Delta\vec{x}_k = -J(\vec{x}_k)\)即可保證函數朝減小的方向改變，值得注意的是，日常的負梯度方向泛化表示就是\(-J(\vec{x}_k)\)，這種形式在SLAM中使用較多，在機器學習領域的話，更多的是使用梯度就好，下面還是簡單總結下從參考博客1里學習到的梯度下降的一些變種：梯度下降，隨機梯度下降，mini-batch梯度下降以及批梯度下降。
在機器學習領域，數據集都很大，通常僅考慮使用一階收斂；而在SLAM領域，數據量也不算大，為保證參數收斂到最優，通常還需要考慮二階收斂。
梯度下降法
優化思想：
用當前位置負梯度方向作為搜索方向，因為該方向為當前位置的最快下降方向，所以也被稱為是”最速下降法“。最速下降法越接近目標值，步長越小，前進越慢。
缺點：
（1）靠近極小值時收斂速度減慢；
（2）直線搜索時可能會產生一些問題；
（3）可能會“之字形”地下降。
批梯度下降、mini-batch梯度下降與隨機梯度下降
在機器學習中，考慮到訓練樣本數過於龐大的問題，衍生出了批梯度下降、mini-batch梯度下降以及隨機梯度下降三種方案來收斂loss。
批梯度下降
批梯度下降每迭代一步，會用到訓練集所有的數據，來得到一個全局最優解。
如果樣本數\(m\)很大，那么其迭代速度會相當的慢。設訓練樣本數為\(m\)，\(\vec{x}\)為\(n\)維向量，一次迭代需要把\(m\)個樣本全部帶入計算，迭代一次計算量為\(m*n^2\)。
隨機梯度下降
因此，出現了隨機梯度下降：
隨機梯度下降每次迭代只使用一個樣本，迭代一次計算量為\(n^2\)，當樣本個數\(m\)很大的時候，隨機梯度下降迭代一次的速度要遠高於批量梯度下降方法。
兩者的關系可以這樣理解：隨機梯度下降方法以損失很小的一部分精確度和增加一定數量的迭代次數為代價，換取了總體的優化效率的提升。增加的迭代次數遠遠小於樣本的數量。
mini-batch梯度下降
由於在樣本集比較大，比較復雜的時候，隨機梯度下降在訓練過程中出現loss波動較大、多次迭代后仍不能收斂的問題，因此出現了mini-batch梯度下降的方法，其實也就是將整個樣本集\(m\)分成多個mini-batch，每次迭代求取一個mini-batch的批梯度下降結果進行更新，遍歷完所有個mini-batch后完成一個epoch的更新。
總結
批量梯度下降---最小化所有訓練樣本的損失函數，使得最終求解的是全局的最優解，即求解的參數是使得風險函數最小，但是對於大規模樣本問題效率低下。
隨機梯度下降---最小化每條樣本的損失函數，雖然不是每次迭代得到的損失函數都向着全局最優方向， 但是大的整體的方向是向全局最優解的，最終的結果往往是在全局最優解附近，適用於大規模訓練樣本情況。
mini-batch梯度下—-單次迭代所考慮的樣本數為mini-batch大小（介於1和\(m\)之間，通常為\(2^n\)），一個epoch所經歷的迭代次數為\(\frac{m}{sizeof(minibatch)}\)（而非批梯度下降中的1或者隨機梯度下降中的\(m\)）
二階收斂
——通常需要考慮前面泰勒展開的二階項
牛頓法
思想很簡單，目標需要求最小的\(\Delta\vec{x}_k\)，直接對\(F(\vec{x}_k+\Delta\vec{x}_k)\approx{F(\vec{x}_k) + J(\vec{x}_k)^T\Delta\vec{x}_k+\frac{1}{2}\Delta\vec{x}_k^TH(\vec{x}_k)\Delta\vec{x}_k}\)等式右側對\(\Delta\vec{x}_k\)求導，等於0即可，最終得到下式：
具體思想就是，每迭代一步，需要求解一個線性方程組的解\(\Delta\vec{x}_k\)，求這個解，需要構造上面的線性方程組，也就是每迭代一次需要求一次Hessian矩陣\(H\),復雜度是很高的，所以出現了下面這些改進方法：擬牛頓法，高斯牛頓法，列文伯格—馬誇爾特法等。
擬牛頓法
——摘錄於引用3和4
在牛頓法中，需要求解巨大的\(H\)矩陣的逆來更新\(\Delta\vec{x}\)。因此，擬牛頓法改進思路是，如何通過已有的信息來從\(H_{n-1}^{-1}\)更新\(H_n^{-1}\)，而不是每次迭代都實打實的去求解一個\(H^{-1}\)，具體思路如下：
注意看第2和3部分的思路：
- Store the input and gradient deltas
- 計算\(\nabla{f(\vec{x}_{n+1})}\);
- \(\vec{s}_{n+1} = \Delta\vec{x}\);
- \(\vec{y}_{n+1} = \nabla{f(\vec{x}_{n+1})} - \nabla{f(\vec{x}_{n})}\)
- 將\(\vec{s}_{n+1}\),\(\vec{y}_{n+1}\)以及\(H_n^{-1}\)傳入QuasiUpdate計算新的\(H_{n+1}^{-1}\)
BFGS
——命名來源於Broyden–Fletcher–Goldfarb–Shanno (BFGS)四個人的名字縮寫
主要提出的是QuasiUpdate部分的計算方法：
算法偽代碼如下：
LBFGS
BFGS擬牛頓近似算法雖然免去了計算Hessian矩陣的煩惱，但是我們仍然需要保存每次迭代的\(\vec{s}_n\)和\(\vec{y}_n\)的歷史值。這依然沒有減輕內存負擔，要知道我們選用擬牛頓法的初衷就是減小內存占用。
L-BFGS是limited BFGS的縮寫，簡單地只使用最近的m個\(\vec{s}_n\)和\(\vec{y}_n\)記錄值。也就是只儲存\(\vec{s}_n, \vec{s}_{n-1}, ..., \vec{s}_{n-m-1}\)和\(\vec{y}_n, \vec{y}_{n-1}, ..., \vec{y}_{n-m-1}\),用它們去近似計算\(H_n^{-1}\vec{g}_n\)。初值\(H_0^{-1}\)依然可以選取任意對稱的正定矩陣。
Notions:
What form should QuasiUpdate take? Well, if we have QuasiUpdate always return the identity matrix (ignoring its inputs), then this corresponds to simple gradient descent, since the search direction is always ∇fn. While this actually yields a valid procedure which will converge to x∗ for convex f, intuitively this choice of QuasiUpdate isn’t attempting to capture second-order information about f.
Note that the only use we have of the hessian is via it’s product with the gradient direction.
高斯牛頓
——阻尼牛頓法
從式(3)可知，用牛頓法求解非線性最小二乘問題每迭代一次，需要求解一次Hessian矩陣，高斯牛頓嘗試用一階的Jacobian矩陣去逼近二階的解，思路如下：
對於\(\min_{\vec{x}}F(\vec{x})=\frac{1}{2}||f(\vec{x})||_2^2\)，選擇對非線性函數\(f(x)\)進行一階泰勒展開，得到\(f(\vec{x}_k+\Delta\vec{x}_k)\approx{f(\vec{x}_k) + J(\vec{x}_k)^T\Delta\vec{x}_k}\)，顯然此時的\(f(\vec{x}_k+\Delta\vec{x}_k)\) 對 \(\Delta\vec{x}_k\) 是線性函數；
因此，高斯牛頓的核心思想是同非線性最小二乘求解的本質思路一致：求解最小的 \(\Delta\vec{x}_k\)，但處理邏輯是將非線性函數\(f(\vec{x})\)泰勒展開為關於 \(\Delta\vec{x}_k\)的線性函數，將該問題轉換為求解關於 \(\Delta\vec{x}_k\)的線性最小二乘問題；
后面邏輯就很簡單了，像求線性最小二乘問題一樣，展開這個二次函數，得到以下式子：
對式(4)求導，等於0，有：
有：
高斯牛頓通過求解上面關於\(\Delta\vec{x}\)的線性增量方程組，即可在每一步迭代一個增量用於更新；
從另一個角度看來，高斯牛頓法用\(JJ^T\)來逼近Hessian矩陣H，規避每一次迭代都要求Hessian矩陣的問題。
然而，上面的\(JJ^T\)在問題比較病態的時候或是奇異矩陣的時候，\(JJ^T\)的逆就會有問題，因此，出現了列文伯格-馬誇爾特法，在展開點加個橢球限制。
列文伯格——馬誇爾特（簡稱LM法）
如果說，上面那些方法都是無約束的優化方法，那LM則是加了不等式約束的一種優化方法。
因為高斯牛頓中近似二階泰勒展開只在展開點附近小范圍內有較好的近似效果，LM對這個展開點增加了一個trust region，范圍用近似模型跟實際函數的差異來確定，定義如下公式來表征近似好壞程度：
而LM的算法偽代碼則如下：
用拉格朗日乘子把約束放到目標函數中則有：
並對式5對\(\Delta\vec{{x}}\)求導等於0，推導出線性增量方程有：
\(D\)通常為非負數對角陣，考慮\(D=I\)，則有\((H+\lambda I)\vec{x}_k=g\),則可作下述分析：
以下摘自SLAM十四講第2版，Page 131:
- 在LM中，當\(\lambda\)比較小的時候，H占主要作用，說明二次近似模型在該范圍內比較好，LM在此時更接近高斯牛頓法；當\(\lambda\)比較大的時候，\(\lambda I\)占主要作用，LM跟接近於梯度下降法，說明二次近似模型在該范圍內不夠好.
- 因此，LM可以起到規避在高斯牛頓中因為線性方程組非奇異或是病態問題導致的增量不穩定問題；
- 當問題較良好的時候，使用高斯牛頓進行優化，當問題較病態的時候，使用LM進行優化。
線性最小二乘問題
線性最小二乘問題通常都可以展開為以下形式：
然后通過對\(\vec{x}\)求導，等於0，即可把問題轉換為求解\(A\vec{x}=\vec{b}\)的解形式；
值得注意的是，在求解非線性最小二乘問題的時候，也會需要將問題轉化為求一個線性方程組的形式，所以下面的內容主要是整理一些求解線性方程組的解方法。
直接法
主要用到的就是LU分解、LUP分解、Cholesky分解以及QR分解。
LU, LUP, Cholesky分解
——\(A\)為方陣情況下的矩陣分解
對矩陣\(A\)的要求：
LU分解、LUP分解對矩陣\(A\)的要求是非奇異即可；
Cholesky分解對矩陣\(A\)的要求是對稱正定矩陣；
數值穩定性：
Cholesky分解是穩定的，LUP穩定性隨着n的增加會有誤差存在，LU穩定型較差；
Notion:
- LU分解就是將\(A\)分解為上三角矩陣\(U\)和下三角矩陣\(L\)
- LUP分解針對LU分解中遍歷矩陣行作消元的時候，對應行對角線元素為零的情況作了處理(加了行置換(pivoting))以增加矩陣分解的穩定性，所以叫做LUP分解；
- Cholesky分解是LU分解的特例，對於\(A\)矩陣為對稱正定矩陣進行分解的過程，分解結果滿足\(U=L^T\)；
下面博客有很詳細的從原理到流程到算法復雜度介紹：
QR分解
——\(A\vec{x}=\vec{b}\)為欠定時的矩陣分解
將\(A\in{R}^{m\times n}, m \geq n\)分解為一個正交矩陣\(Q \in R^{m \times m}\)和\(R\)矩陣，其中：
分解方法的幾點note：
- Gram–Schmidt Orthogonalization（常規矩陣論中施密特正交化）
- python matlab中使用的QR分解方法
- Givens Rotations
- 當A是稀疏矩陣，Givens Rotations大小為0的元素可以直接被忽略。另一個優勢是，Givens Rotations更容易並行化，因為Givens Rotations只對兩個元素進行操作，處理不同列的時候可以完全的獨立。
以上摘自同一個作者的QR分解：
Note:
對於超定方程組是存在最小二乘解的: \(\vec{b}=(A^TA)^{-1}A^T\vec{y}\)
詳細證明可見：
https://wenku.baidu.com/view/abb9dd06eefdc8d376ee32ca.html
上面作者在他的專欄里也有相關描述：
迭代法
主要在\(A\)矩陣維度比較大比較稀疏的情況常用，通過迭代的方式求得一個近似解。
一階線性定常迭代法
將矩陣\(A\)拆分成\(D-P\), \(D\)為可逆矩陣；
則有\(A\vec{x}=\vec{b}\) => \((D-P)\vec{x}=\vec{b}\) => \(\vec{x}=D^{-1}P\vec{x}+D^{-1}\vec{b}\)
由此則有一般迭代公式：
\(\vec{x}^{(k+1)}=G\vec{x}^{(k)}+d\)
此式收斂的充要條件是：
\(G\)的譜半徑小於1(譜半徑的定義式\(G\)的最大特征值)
Jacob迭代
——矩陣\(A\)為可逆矩陣
若\(D\)為\(A\)的對角元素構成的對角矩陣，且對角線元素全不為零，當\(A\)的最大特征值(也稱譜半徑)小於1時，雅各比迭代收斂
具體思路為：
將矩陣\(A\)拆分成\(L+D+U\), \(L\)為嚴格下三角矩陣，\(D\)為\(A\)的對角線元素構成的對角矩陣，\(U\)為嚴格下三角矩陣，顯然此時的迭代公式即可變為：
用向量分量表示即有：
其中，\((i=1,2,..,n; k=0,1,2,...)\)，\(x^{(0)} = (x_1^{(0)}, x_2^{(0)}, ..., x_n^{(0)})^T\)為初始向量
整個表達還是比較清楚直觀的。每次迭代需要計算一次矩陣與向量的乘法，需要存儲\(\vec{x}^{(k)}\)和\(\vec{x}^{(k+1)}\)
高斯賽德迭代
——Gauss–Seidel
主要內容：在求的\(x_i^{(k+1)}\)之前，\(x_{0}^{(k+1)}, ..., x_{i-1}^{(k+1)}\)就都已經得到，但根據雅各比的計算式\(x_i^{(k+1)}=\frac{1}{a_{ii}}[b_i-\sum_{j=1, i \ne j}^n{a_{ii}x_j^{(k)}}]\)，還在用\(\vec{x}^{(k)}\)來計算\(x_i^{(k+1)}\)；
高斯賽德的思想就是在求\(x_i^{(k+1)}\)的時候，將已有變量划分為已經更新的和沒有更新的變量，對於已經更新的就用求得的\(x_{0}^{(k+1)}, ..., x_{i-1}^{(k+1)}\)來計算，對於未更新的就還是用\(\vec{x}^{(k)}\)來計算。
這樣一來可以“在巨人的肩膀上前行”加快收斂速度，二來可以摒棄已用完的舊變量，從而優化雅各比迭代內存占用和收斂速度而提出，用矩陣進行表達就是：
一些迭代收斂判定可用的條件：
充分條件：
- 矩陣\(A\)的某種范數(e.g. \(||A||_{1}\), \(||A||_{\infty}\)等)小於1，因為可證明\(A\)的譜半徑總小於\(||A||_{r}\); 此條件可用於快速分析迭代的收斂性，因為在有些時候\(||A||_{r}\)比\(A\)的特征值分解簡單。
- \(A\)矩陣按行嚴格對角占優或按列對角嚴格占優，即滿足:\(|a_{ii}| \gt \sum_{j=1, i \ne j}^n{|a_{ij}|}\)或者\(|a_{ii}| \gt \sum_{j=1, i \ne j}^n{|a_{ji}|}\)
充要條件：
- 矩陣\(A\)的最大特征值（譜半徑）小於1；
滿足上述條件之一：高斯賽德或者雅各比迭代都會收斂
若\(A\)為對稱正定矩陣，高斯賽德迭代一定收斂；
逐次超松弛迭代(SOR)
其實就是在高斯賽德迭代的思想上加了松弛因子，看它的矩陣表示的更新式即可清楚了解SOR的思路：
當\(\omega\ge1\)的時候，即為超松弛迭代；
當\(\omega=1\)的時候，即為高斯賽德迭代；
當\(\omega\le1\)的時候，即為低松弛迭代；
涉及以下變量：
\(\vec{x}^{(0)}=(x_1^{(0)}, ..., x_n^{(0)})^T\)
\(x_i^{(k+1)}=x_i^{(k)} + \Delta x_i\)
\(\Delta x_i = \omega \frac{b_i - \sum_{j=1}^{i-1}{a_{ij}x_j^{(k+1)}} - \sum_{j=1}^n{a_{ij}x_j^{(k)}}}{a_{ii}}\)
\(i = 1, 2, ..., n; k = 0, 1, ...\)
\(\omega\)為松弛因子；
當\(\max_{1 \le i \le n}|\Delta x_i| = \max_{1 \le i \le n}|x_i^{(k+1)}-x_i^{(k)}| < \epsilon\)迭代終止。
總結：
摘自：
SOR和SSOR迭代
Jacobi迭代是將第i個方程中的第i個分量單獨挪到一邊，每個方程單獨迭代的過程。Gauss-Seidel迭代是在Jacobi迭代過程中，對於迭代變量，求解方程中若某分量產生了新值用新值而不用舊值的過程，SOR迭代是在新值和舊值之間取個加權平均的過程。
共軛梯度法
主要思想
對於\(n \times n\)維矩陣\(A\)，其思想就是找到\(n\)個兩兩共軛的共軛方向，每次沿着一個方向優化得到該方向上的極小值，后面再沿其它方向求極小值的時候，不會影響前面已經得到的沿哪些方向上的極小值，所以理論上對\(n\)個方向都求出極小值就得到了\(n\)維問題的極小值；
在此博客中有較為清晰的推導與介紹，直接參考即可：共軛梯度法詳細推導分析