Jinyu Li a personal journal

最小二乘问题(七)

在各种最小二乘问题中,还有一类比较重要的问题。我们以 Wahba’s Problem 举个例子:

已知三维空间中的点云 $P = { p_1, p_2, \dots, p_n }$ 和 $Q = { q_1, q_2, \dots, q_n }$ ,求旋转矩阵 $R \in \mathbb{SO}(3)$ 将 $P$ 旋转后与 $Q$ 吻合,即求解

\[\min \sum_i \|Rp_i-q_i\|^2\ \text{s.t.}\ R\in\mathbb{SO}(3).\]

注意到哪里不同了么?对的,我们的解空间有了约束。对于这种带约束的最小二乘问题,是无法直接套用我们前面的方法的。

可能这里还有人在关心 Wahba’s Problem 怎么解,因为它有着明确且重要的几何意义,可以被用于各种类型的姿态配准。幸运的是,对于这个特殊的问题,我们可以用 SVD 分解直接求出它的解,无需经过迭代方法。这个算法叫做 Kabsch Algorithm 。

那么,其它带有类似约束的问题该怎么办呢?

首先,我们要说说“类似”指的是什么。我们前面问题里的约束是一个旋转群,它可以等价的转换成这样的两个约束:

\[R^TR = RR^T = I,\ \mathrm{det}\left(R\right) = 1.\]

注意到,这是两个等式约束,它定义了在 $\mathbb{R}^{3\times3}$ 空间里的一个曲面。我们的解正是被约束在了这个曲面上。

数学上我们管这个曲面叫做流形(Manifold)。它除了有形状,在每个点处还有微分性质。我们接下来要关心的正是这类约束在微分流形上最小二乘问题。不过这里我们先开个头,正文我们下回继续。

MDZZ

换了一种新程序语言之后连键盘都不会敲了……

然而当发现它其实就是披了马甲的 Java 之后智商马上又回来了……

这是什么毛病?

最小二乘问题(六)

前面提到的最小二乘有一个共同的假设,就是残差向量遵循高斯分布。

那么,如果不满足高斯分布假设呢?

通常我们模型的噪音都是符合高斯分布的,对于模型中蕴含的非高斯分布情形,我们不在这里讨论了。我们介绍一种在高斯假设下依旧会引入非高斯分布残差的情形,并介绍在最小二乘中处理它的一种手段。

通常我们在建模时认为数据点都是由模型产生并带有一定的噪音。但在真实中,数据可能掺杂着不属于当前模型的点,叫做外点(outlier)。这些外点会污染数据,表现上便是出现了不符合高斯分布的噪音。由于我们采用平方项,当与模型偏差很远的点存在时,它会对模型的估计带来严重的不良影响。

为了减小这种影响,我们可以采用加权的最小二乘,即

\[\min \sum w_i\|r_i\|^2.\]

理想情况下,只要令外点对应的权值为 0 即可从优化中忽略掉它们。

然而真实情况中,我们可能不知道外点都有哪些,因此我们需要一个可以自动降低外点影响的加权方案。”M-估计”便是一种常用的方法,它在优化中引入了“鲁棒函数”来根据误差的大小自动调整权重:

\[\min \sum \rho_i(\|r_i\|).\]

这里的 $\rho_i$ 便是鲁棒函数,如果 $\rho_i(x) = x^2$ ,我们就得到了一般的最小二乘问题。为了减小外点(带有巨大误差的点)的影响,一般要取 $\lim_{x\to\infty}\frac{\rho(x)}{x^2} \to 0$ 的函数。也就是当误差越大,它对整体优化的影响就越弱于普通的平方项。但同时,$\rho(x)$ 还应该满足从 0 向两个方向增长,否则优化可能会收敛到错误的结果(比如直接拒绝掉全部的数据点从而得到更小的能量)。

常见的鲁棒函数选择很多,比如 Huber 函数、Cauchy 函数等等。网上资料很丰富,这里就不列举了。

加入了鲁棒函数后,我们的问题就不是最小二乘了。但前面非线性最小二乘中,我们通过线性化得到了线性最小二乘子问题。这里我们也尝试一下转化它成为我们熟悉的形式。

我们比较带鲁棒函数的能量和一般的最小二乘能量的导数,根据极值的一阶条件我们知道:

\[\begin{aligned} \frac{\partial \|r\|^2}{\partial x} &= 2r^T\frac{\partial r}{\partial x} = 0 \\ \frac{\partial \rho(\|r\|)}{\partial x} &= \frac{d\rho}{dx} \frac{1}{\|r\|}r^T\frac{\partial r}{\partial x} = 0. \end{aligned}\]

如果我们记 $w(x) = \frac{d\rho}{dx}\frac{1}{x}$,进一步地令常系数 $w’= w(\|r\|)$ ,然后构造加权的最小二乘优化 $\min \frac{w’}{2}\|r\|^2$,它的 Jacobian 是:

\[\frac12\frac{\partial w’\|r\|^2}{\partial r} = w’r^T\nabla r.\]

可以发现,它与前面带有鲁棒函数的最小二乘优化形式是相同的。

也就是说,我们可以在每一轮迭代时,采用当前的 $\|r\|$ 计算加权值。固定加权值后,将当前迭代中的优化问题看做一个加权最小二乘问题进行优化。这时,带鲁棒函数的 M-估计与带更新的加权最小二乘是等价的。

Cholesky 分解(番外二)

为了模拟 Cholesky 分解的 fill-in 现象,写了一个矩阵消元小游戏。

它是 Web 界面,地址是 http://epipolar.tumblr.com/matrix-elimination 。这也是我第一个用 javascript 和 processing (p5.js) 编写的程序。

欢迎尝试~

最小二乘问题(五)

上回末有人搞事儿,希望加入非高斯的噪声……但我们偏不!我们要先看一些更重要的扩展。

我们记 $r(x) = Ax-b$ ,这是我们的(线性)残差函数。那么标准的最小二乘问题就是 $\min \|r(x)\|^2$ 。残差函数关于变量 $x$ 的导数(Jacobian)记作 $J_r$ ,我们知道 $J_r = A$,于是最小二乘的标准方程还可以写成 $J_r^TJ_rx = J_r^Tb$ 。

如果 $r(x)$ 并非关于 $x$ 的线性函数呢?此时我们得到了年轻人的第一个非线性最小二乘问题!

它与线性最小二乘最大的区别便在于:非线性最小二乘通常不是凸优化问题。这是什么意思?恩我们不在这里解释了。那这个区别有什么用?恩它没什么用,它是来给我们制造新的困难的,这个困难就是:对于线性最小二乘问题,你可以求解标准方程得到全局最优解;然而对于非线性最小二乘问题……

  • 它有标准方程么?是什么形式?
  • 它的解如何寻找?可以求解标准方程么?
  • 它有全局最优解么?

恩……对于一般的非线性最小二乘问题,通常我们只能找到一个局部最优解,至于局部最优是不是全局最优,这就要靠手段加运气了。它的解要如何寻找呢?答案是从一个初始的猜测出发,一步步走到更优的解,直到收敛——这也叫做迭代优化方法。至于标准方程,恩我们下面来看看非线性最小二乘的标准方程。

如果残差函数有具体的表达式,那当然是可以更好地分析,不过这里我们从它的一般形式入手。一个问题……它不是线性的,怎么办?线性化咯!

于是我们将残差函数 $r(x)$ 用它的一阶近似表达,也就是

\[r(x_0+\Delta x) = r(x_0)+J_r\Delta x.\]

这样,当我们从一个初始值 $x_0$ 出发时,我们可以寻找一个让一阶近似更优的变化量 $\Delta x$,这样我们便“大约”将原始问题优化了一些。

这个思路写成最优化形式,就是

\[\min_{\Delta x} \|J_r\Delta x + r(x_0)\|^2.\]

它是个关于 $\Delta x$ 的线性最小二乘问题,套用我们前面的知识,我们知道它的标准方程是

\[J_r^TJ_r\Delta x = -J_r^Tr(x_0).\]

前面我们已经提到了,整个问题需要一步步不停地迭代,所以说这里我们只是完成了一个小目标,接下来我们要更新变量的值并进一步优化。

所以完整的非线性最小二乘优化算法是这样的:

  1. $x \gets x_0;$
  2. $\Delta x \gets \arg\min_{\Delta x} |J_r\Delta x+r(x)|^2;$
    • 求解标准方程 $J_r^TJ_r\Delta x = -J_r^Tr(x);$
  3. $x \gets x+\gamma \Delta x;$
  4. 回到 2 ,除非 $|\Delta x|<\epsilon$ 。

这个算法叫做 Gauss-Newton 算法。而其中关于 $\Delta x$ 的标准方程也叫作 Gauss-Newton 标准方程。$\gamma > 0 $ 是一个控制步长的系数。

那么如果有多项呢?如果残差遵循不同的高斯分布呢?线性化,然后回忆前面的故事吧!