最小二乘问题（十二）

当标准方程左侧的 $J^TJ$ 秩亏时，我们知道它不可逆。但我们不妨假设在数值计算中“阴差阳错”地使它可逆了，会发生什么呢？

导致这种现象通常的原因是浮点数的舍入误差。我们不去研究背后具体的机制了，假设 $\mathrm{det}(J^TJ) = \epsilon$，那么我们知道 $\|\Delta x\| = \|(J^TJ)^{-1}\|\|J^Tb\| = \|J^Tb\|/\epsilon$ 。这意味着在 $J^TJ$ 秩亏或者接近秩亏的时候，我们会得到一个巨大的迭代步！

当我们的问题存在很多局部最优时，这样的迭代步会导致我们的算法跳过当前局部最优的区域，进入到其它区域。很多时候这是我们不希望的。

回到 LM 算法，我们将 LM 的子优化问题放在概率框架下观察。其实这个优化的含义对应于：已知 $\Delta x$ 满足某个参数的高斯分布，求在这个先验条件下，原始线性化问题的最大后验概率解。我们知道 $\lambda$ 反比于这个高斯分布的方差，由此可以想到在原问题 $J^TJ$ 接近奇异时，我们希望 $\Delta x$ 被更加严格地约束，因此要提高 $\lambda$ 的值。反之，当 $J^TJ$ 性质良好时，我们愿意给 $\Delta x$ 更大的自由度，此时 $\lambda$ 可以适当减小。

这个高斯分布描述了我们对 $\Delta x$ 或者说对当前 $x$ 的置信区间，每一轮我们在当前的置信区间内寻找更优值，因此也说 LM 算法属于信赖域算法。

深藏功与名

之前到隔壁去取一些数据。这是一些建筑模型的数据，不过很奇怪的是，模型的 Bounding Box 经常会大出去一截，导致我需要重新计算。

一开始觉得是我的读取方式有误，于是问他们 BB 的解析算法。顺藤摸瓜发现他们是这么计算每个包围盒的坐标的：

float bb_min[3], bb_max[3];
for(int i = 0; i < 3; ++i) {
  bb_min[i] = std::numeric_limits<float>::max();
  bb_max[i] = std::numeric_limits<float>::min(); // min 给出的是最小的绝对值
}
for(int j = 0; j < model.vertex_num(); ++j) {
  for(int i=0;i<3;++i) {
    bb_min[i] = std::min(bb_min[i], model.vertex(j)[i]);
    bb_max[i] = std::max(bb_max[i], model.vertex(j)[i]); 
  }
}

看到我加注释的那行忽然就明白这个错误的原因了。而他们也一下子解决了一个据说两年没发现原因的 bug 。所以说，遇到不熟悉的接口，一定要查文档。

于是开心地拿到了正确的数据离开了，深藏功与名……

最小二乘问题（十一）

我们来看非线性最小二乘问题的线性化子问题：

假如我们不能保证 $J_r$ 列满秩，我们套用上一篇末尾的正则化方法，将这一子问题改写为：

这时，我们在每一步迭代时便会在“尽可能最小化当前线性化函数”和“尽可能保持步伐稳健”之间进行折衷。那么这会影响到我们最终收敛到局部极值么？注意看，如果我们已经在局部极值附近，此时 $\|\Delta x\|^2\to 0$，也就是说第二项自然就消失了。换句话说，在局部极值点附近，我们的问题与 Gauss-Newton 中的线性化子问题是相同的。

再来看看这个常数 $\lambda$ ，当 $\lambda \to 0$ 时，我们得到了 Gauss-Newton，此时我们通过二阶近似来逼近目标函数。当 $\lambda \to \infty$ 时，根据标准方程

这意味着，此时我们的优化步采用的是梯度下降步。

也就是说，通过调节 $\lambda$ ，我们的优化算法在 Gauss-Newton 步和梯度下降步之间平衡。此外，采用了这种方法后，在每次子问题求解后，更新变量时无需再采用步长控制，可以直接使用 $x\gets x + \Delta x$ 。 这是因为 $\lambda$ 自带了步长控制的效果。

这个方法被称为 Levenberg-Marquardt 算法，是一种广泛使用的非线性最小二乘问题的数值优化算法。

那么，吹了一篇 LM 算法，究竟什么时候要 GN 步，什么时候要梯度下降步呢？Jacobian 不满秩时 GN 有多糟糕呢？ 为啥 LM 就好了呢？ 我们下一篇介绍。

小软件：Hourglass

Hourglass - The simple countdown timer for Windows

Hourglass is the most advanced simple countdown timer for Windows. Just enter a time in just about any format, and hit Enter.

一个漂亮并且精简的倒计时小工具~

之后也会陆续记录各种遇到的实用小软件，免得每次想用时候记不得叫啥~

最小二乘问题（十）

在我们这一系列文章的最一开始，我们对标准形式的 $\min \|Ax-b\|^2$ 有着如下的要求：

这其中对 $A$ 的满秩要求是非常重要的，因为我们优化这一问题的方式是求解标准方程，得到：

结合前面的条件，我们知道方程有解当且仅当 $A^TA$ 可逆。而如果 $\mathrm{rank}(A)<n$ ，$\mathrm{rank}(A^TA)\leq\mathrm{rank}(A)<n$ 一定不可逆。

在秩亏的情形下，我们的解 $x$ 有无穷多，这是因为对于任意的 $z\in\mathrm{ker}(A^TA)$ ，$A^TA(x+z)=A^TAx$ 。此时，我们的问题是欠约束的。

为了解决秩亏问题，我们可以引入 $A$ 的伪逆 $A^\dagger$ 。将问题的解写作：

利用 SVD 分解，若矩阵 $A=U\Sigma V^T$，则 $A^\dagger = V\Sigma^\dagger U^T$ ，这里 $\Sigma^\dagger = \mathrm{diag}(\mathrm{inv}(\sigma_1), \mathrm{inv}(\sigma_2), \dots, \mathrm{inv}(\sigma_n))$. 相比一般的矩阵求逆，由于秩亏矩阵存在为 0 的奇异值，我们特别引入了专用的倒数函数：

我们回到原始的最小化问题，$Ax-b$ 几何上对应于求 $b$ 在 $A$ 的列空间上的投影，利用上面的 SVD 分解结果，我们知道 $U^T b$ 中对应非零奇异值的行便对应了这个投影，剩余的部分是垂直于 $A$ 列空间的余项（$Ax$ 永远无法表达的部分）。完整的 $A^\dagger b$ 恰好给出了不包含余项的部分，同时它在 $A$ 的右零空间投影也是零，所以它是对 $Ax-b=0$ 的最好逼近，同时也保证了 $\|x\|$ 最小。

有关上面内容的详细证明，可以参考其它资料，或者尝试自行推导。从结论上讲，$x=A^\dagger b$ 求解了下面的问题：

这种处理方法虽然可以解决任意的线性最小二乘，但它的缺点也是明显的：SVD 分解通常需要很大的运算量。通过观察问题的解，我们“希望”最小化原始问题的同时最小化 $x$ 的模长，这就诱使我们转而采用下面的带正则化最小二乘：

它的标准方程为 $(\lambda I+A^TA)x = A^Tb$ ，当 $\lambda>0$ 时，$x$ 左侧的系数矩阵一定可逆。为了能够更好近似原始最小问题的结果，我们需要 $\lambda$ 尽可能小。否则我们就会在两项之间进行一个折衷，最终的结果并不能真正最小化 $\|Ax-b\|^2$ 。

但在非线性最小二乘中，我们会看到这一正则化会带来一种新的最小化算法。它可以保证我们收敛到局部最优，同时它比标准的 Gauss-Newton 具有更好的稳定性。