题外话：
从开始学习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\)维问题的极小值；
在此博客中有较为清晰的推导与介绍，直接参考即可：共轭梯度法详细推导分析