Epipolar a personal journal

最小二乘问题(十四)

与固定变量相对的,我们还可以“忽略”某个变量。这里的忽略当然不是从我们的问题中去掉这个变量,而是允许它取任何可能的值,在它全体可能取值范围内求解最优的剩余变量来最小化我们的目标函数。

我们把这种“忽略”叫做边缘化(Marginalization)。在一个一般的非线性最小二乘中直接边缘化变量的情形比较复杂,因此这里我们只考虑在每一步线性子问题中边缘化掉某个梯度方向。

我们依旧套用上一篇中的标准方程。完整的系统中的标准方程包含了全体梯度方向。如果我们想要边缘化某个方向,不妨是 $\Delta x_1$ 。意味着我们想得到一个新的标准方程,在这个标准方程中不包含 $\Delta x_1$ 作为变量。

“可是上一篇中我们不是已经有了一个新的并且没有 $\Delta x_1$ 的标准方程了吗?”是这样没错,但回忆上一篇,我们得到的那个问题中,$x_1$ 变量被固定了,也就是说并不是 $\Delta x_1$ 被忽略了,而是它被固定为了 $\Delta x_1 = 0$ 。而这里我们所希望的是让它可以随心所欲,因此我们需要用另外一个技巧:高斯消元。

我们将原始的标准方程的矩阵重新分块,然后将标准方程写为下面这样:

我们用 $\bar{1}$ 表示除去 1 之外的部分,$\Delta x_{\bar{1}}$ 便是我们想要保留的全部变量,其它部分也同样类比。我们对方程进行块状行变换,将方程左右的第一行($x_1$ 对应的行)左侧乘上$-J_{\bar{1}}^TJ_1(J_1^TJ_1)^{-1}$ 后加到第二行上,于是第二行便成为了:

这个形式看起来复杂,但是左侧的系数矩阵其实正是原先完整的系数矩阵关于 $J_1^TJ_1$ 部分的 Schur 补。这一操作我们在前面的文章中介绍过。在这个新的方程中,我们的变量不再包含 $\Delta x_1$ ,但可以发现无论是左侧的系数还是右侧的余项中,都包含了来自 $\Delta x_1$ 的梯度信息 $J_1$ 。正是因为我们加入了这部分信息,才允许我们可以“忽略”掉对应的变量直接求解。

如果我们求解了这第二行的部分,接下来是可以将 $\Delta x_{\bar{1}}$ 代入回第一行继续解出对应当前最优的 $\Delta x_1$ 的。不过如果我们强行让 $\Delta x_{\bar{1}} = 0$ 并代入,得到了什么?没错,这时我们就得到了上一篇文章中描述的固定变量情形了。

进一步的,我们可以给 $\Delta x_{\bar{1}}$ 赋予任意需要的值,然后利用第一行求解对应的 $\Delta x_1$ ,这实际上对应了给定前者作为条件时后者的最大似然估计。这一系列的概率背景就不在这里展开介绍了,有机会在另外的文字里描述一下。

径向畸变的效果

Effect of Radial Distortion

Explanation: Black dots should form a rectangular grid. However, if the camera lens has radial distortion, image will be distorted. This animation simulates the photo under different radial distortion parameters. More specifically, the distorted image point is given by $p \mapsto L(r)\times p$, where $L(r) = 1+k_1\times r$, and $r = \|p-c\|$.

最小二乘问题(十三)

我们考虑一种情况:优化过程中,如果我们要固定某一个变量的值,问题需要如何改变。

为了方便研究,对于一般的最小二乘问题,我们把它所有项的残差向量 $r_i$ 纵向连接成一个残差向量 $r$ ,得到下面的最小二乘问题:

我们用 $J_i$ 表示 $r$ 关于 $x_i$ 的 Jacobian 矩阵。那么这个问题的标准方程就是:

假如我们希望将某个变量固定,不失一般性地,设这个变量为 $x_1$ 。由于新的优化问题中 $x_1$ 被视为常量,$r$ 的 Jacobian 矩阵将不包含 $J_1$ 项。新的优化问题的标准方程就是:

这个结论可以自然地向固定任意多个变量推广,当我们要固定某些变量时,只需要求解未固定变量对应的区块的子线性系统。

让我们再来看这个线性的子问题,上述的固定实际上可以看成“令 $\Delta x_1 = 0$ ”,这么一来,我们的线性最小化问题中,与它有关的系数就全都不见了。

一般地,如果我们有线性最小化问题 $\min_{x,y} \|(A_1 \; A_2){x \choose y}-b\|^2$。如果我们选择固定 $x=x_0$ ,这个问题就变成了 $\min_y \|A_2y+(A_1x_0-b)\|^2$ ,它的标准方程也会对应地改变。我们可以换一种说法,把这个新的问题叫做“以 $x=x_0$ 为条件时的最小化”,这种操作也因此叫做条件化。

最小二乘问题(十二)

当标准方程左侧的 $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 。所以说,遇到不熟悉的接口,一定要查文档。

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