Jinyu Li a personal journal

Cholesky 分解(续)

Cholesky 分解有一个变体,它无需计算平方根。它就是 LDL 分解。

其实思路也很简单,我们把一个对称正定矩阵 $A$ 分解为 $A=LDL^T$ ,这里 $L$ 是下三角矩阵且对角线上全为1,$D$ 是对角阵且对角线上全大于零。

不妨按照前面的思路对矩阵分块,得到

\[\begin{pmatrix} 1 & 0 \\ l_{21} & L_{22} \end{pmatrix} \begin{pmatrix} d_1 & 0 \\ 0 & D_2 \end{pmatrix} \begin{pmatrix} 1 & l_{21}^T \\ 0 & L_{22}^T \end{pmatrix} = \begin{pmatrix} a_{11} & a_{21}^T \\ a_{21} & A_{22} \end{pmatrix}.\]

计算乘积并且对应项系数相等,就能得到结果:

\[\begin{aligned} d_1 &= a_{11} \\ l_{21} &= \frac1{d_1}a_{21} \\ L_{22}D_2L_{22}^T &= A_{22} - d_1l_{21}l_{21}^T = A_{22}-l_{21}a_{21}^T. \end{aligned}\]

从上面的结果进一步观察就可以知道 LDL 分解与 Cholesky 分解有相同的 fill-in 。从运算量和数值性能上,LDL 要优于原始的 Cholesky 分解。

最小二乘问题(十四)

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

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

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

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

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

\[\begin{pmatrix} J_1^TJ_1 & J_1^TJ_{\bar{1}} \\ J_{\bar{1}}^TJ_1 & J_{\bar{1}}^TJ_{\bar{1}} \end{pmatrix}\begin{pmatrix}\Delta x_1 \\ \Delta x_{\bar{1}}\end{pmatrix} = -\begin{pmatrix}J_1^T \\ J_{\bar{1}}^T\end{pmatrix}r.\]

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

\[(J_{\bar{1}}^TJ_{\bar{1}}- J_{\bar{1}}^TJ_1 (J_1^TJ_1)^{-1} J_1^TJ_{\bar{1}}) \Delta x_{\bar{1}} = -(J_{\bar{1}}^T - J_{\bar{1}}^TJ_1 (J_1^TJ_1)^{-1} J_1^T)r.\]

这个形式看起来复杂,但是左侧的系数矩阵其实正是原先完整的系数矩阵关于 $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$ ,得到下面的最小二乘问题:

\[\min \|r(x_1, x_2, \dots, x_n)\|^2\]

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

\[\begin{pmatrix} J_1^TJ_1 & J_1^TJ_2 & \cdots & J_1^TJ_n \\ J_2^TJ_1 & J_2^TJ_2 & \cdots & J_2^TJ_n \\ \vdots & \vdots & \ddots & \vdots \\ J_n^TJ_1 & J_n^TJ_2 & \cdots & J_n^TJ_n \end{pmatrix}\begin{pmatrix} \Delta x_1 \\ \Delta x_2 \\ \vdots \\ \Delta x_n \end{pmatrix} = -\begin{pmatrix}J_1^T \\ J_2^T \\ \vdots \\ J_n^T\end{pmatrix} r.\]

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

\[\begin{pmatrix} J_2^TJ_2 & \cdots & J_2^TJ_n \\ \vdots & \ddots & \vdots \\ J_n^TJ_2 & \cdots & J_n^TJ_n \end{pmatrix}\begin{pmatrix} \Delta x_2 \\ \vdots \\ \Delta x_n \end{pmatrix} = -\begin{pmatrix}J_2^T \\ \vdots \\ J_n^T\end{pmatrix} r.\]

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

让我们再来看这个线性的子问题,上述的固定实际上可以看成“令 $\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$ 为条件时的最小化”,这种操作也因此叫做条件化。