Epipolar a personal journal

最小二乘问题(四)

前面我们发现各种类型的线性最小二乘最后都可以变成一个最朴素的标准型。那么接下来是不是该讲怎么解标准型最小二乘问题了呢?很遗憾我们不会介绍怎么求解最小二乘问题。

为什么?前面我们提到了,使最小二乘问题最优的解,就是它的标准方程的解。也就是说,我们应该转而研究如何求解线性方程。这就在本系列话题之外了,会另外开辟内容介绍。

今天我们介绍最小二乘问题的概率背景,即最大似然估计(MLE)。

首先解释什么叫最大似然估计。白话说就是:这么多事儿里头就这个最容易发生,那就是这个了。概率上讲就是求使得概率 $P(x)$ 最大的 $x$ 。还有一个比较突出重点的解释,就是最大似然估计是求解随机分布中的众数(mode)。

设 $r\equiv Ax-b$ ,一般我们称 $r$ 为残差。这个名字是有来由的,通常我们认为它是个误差项,并且满足均值为零的高斯分布。在很多问题中,我们认为 $r$ 的每一维都是独立同分布的,在这样的条件下,我们知道 $P\left(r\right) \propto \exp-\frac{r^Tr}{2\sigma^2}$ 。

那么最大化上述概率就是:

由此可见两者是等价的。

接下来,我们来看一个重要的变形:如果 $r$ 的每一维不是独立同分布的,但都是高斯的。此时 $r$ 符合一个一般的多变量高斯分布。如果设 $r$ 的协方差矩阵为 $\Sigma$ ,那么此时 $r$ 的最大似然估计则是:

如果我们引入 Mahalanobis 距离 $\|x\|_P = \sqrt{x^TPx}$ 。令矩阵 $H = \Sigma^{-1}$ ,那么上面的问题可以简写成 $\min \|r\|_H^2$ 。

引入这个记法除了简化记号外有啥别的帮助么?范数都变了呀,是不是现在我们终于得到了一个不能写成线性最小二乘标准型的问题呢?

别着急,如果我们设矩阵 $U$ 满足 $U^TU = H$ ,那么我们便知道

怎么样,又回到最小二乘了吧?前面我们介绍过了 Cholesky 分解,因此我们知道当 $H$ 是对称正定矩阵的时候,对应的 $U$ 肯定能找得到。

……那如果不是高斯分布呢?(你闭嘴!

Cholesky 分解(番外)

前面我们介绍了 Schur 补,它可以看做是块状高斯消元的产物。Cholesky 分解同样可以看做一种特殊的高斯消元,因此它们之间存在着一些联系。

我们将对称正定矩阵 $A$ 按满足 Schur 补的条件分块。得到

与此同时,我们按块进行 Cholesky 分解,将 $A$ 分解为 $LL^T$ 。参考我们系列文章的第一篇可以知道 $L$ 的对应块满足:

如果将 $L_{11}$ 、$L_{21}$ 代入到 $L_{22}$ 中,整理可得:

也就是说在 Cholesky 分解中,每一步消去后得到的剩余方阵恰好是原始矩阵关于消去行列的 Schur 补。因此我们可以直接利用已知的知识来分析和处理计算 Schur 补时产生的 fill-in 情况。

最小二乘问题(三)

上回介绍了三种(线性)最小二乘问题的变形,我们采用了相同的手段,就是构造出新的系数矩阵,然后将原始的问题变成新的系数矩阵下的单项的最小二乘问题。

从上文的三个问题稍加总结,便可以看出原始问题中的系数矩阵在新的系数矩阵里有这样的关系:

  • 每一项的系数矩阵单独占据一行
  • 每个变量的系数矩阵单独占据一列

同时,原始问题中的常数向量直接纵向连接便得到等价问题中的常数向量。

根据这个规则,下面我们介绍一般的多个二次项,多变量的最小二乘问题,以及如何将它化为单个二次项,单变量的最小二乘标准形式。

设一系列残差函数 $r_i(x_1, x_2, \dots, x_n) = \sum_j^n A_{ij} x_j - b_i$ ,$i = 1\dots m$ 。一个一般的线性最小二乘问题具有如下的形式:

定义系数矩阵 $A$ 、常数向量 $b$ 以及向量 $x$,分别为:

那么上面的最小二乘一般形式可以转化为如下的标准形式:

于是我们兜了个大圈子又回来了……

Pointer to Implementation

在 C++ 开发中有一个重要的技巧叫做 Pointer to Implementation (pimpl) 。它被用来隔离接口类和实现类。通常实现类中会包含一些对第三方二进制代码的依赖。而接口类的头文件中不会出现任何与第三方二进制有关的成员。这样一来便使得包含了头文件的用户代码永远不需要知道第三方二进制代码的存在。当这些第三方二进制发生改动时,只要接口类本身的接口(以及头文件)没有任何变动,用户代码侧就不需要重编译。这样就一定程度实现了 ABI 层的兼容。

举个例子,我们可能会有下面的代码:

struct Interface {
    void foo();
    OtherClass m_other;
};

这个接口类在析构时需要考虑到 OtherClass 的内存布局,那么当 OtherClass 发生了二进制级别的变化时,Interface 就不再兼容这些变化了。

我们可以采用下面的方法隔离开它:

class OtherClass;
struct Interface {
    void foo();
    OtherClass *m_pother;
};

由于 m_pother 是一个指针,它的大小是固定的,内存布局也是确定的,因此就无需关心 OtherClass 的具体实现。这样的指针就叫做 pimpl 。

这个实现类需要进行分配和释放,所以问题就多了一些,我们需要在接口类中添加构造和析构函数,并且它们要定义在接口类的实现(.cpp)中。这是因为如果在头文件里直接定义,就需要添加到第三方代码的头文件,于是用户代码里就会出现依赖,这是我们不希望看到的。

随着 C++ 的发展,智能指针也开始普及。利用智能指针进行资源管理有着简单安全的优点,所以一定会有人尝试用智能指针来管理这个 pimpl ,就像下面这样:

class OtherClass;
struct Interface {
    void foo();
    std::unique_ptr<OtherClass> m_pother;
};

但是这么做存在一个问题:std::unique_ptr 需要知道 OtherClass 的析构函数。上面这个接口类会交由编译器自动生成析构函数,在用户的代码里,我们没有 OtherClass 的完整定义,因此 std::unique_ptr 的析构无法生成,也就导致 Interface 不能自动生成析构函数。因此我们不能让编译器自动生成析构函数。

如果要使用智能指针完成 pimpl ,需要像下面这样:

class OtherClass;
struct Interface {
    ~Interface();
    void foo();
    std::unique_ptr<OtherClass> m_pother;
};

在类的接口定义中显式声明析构函数,然后在它的实现代码里:

Interface::~Interface() = default; // 或者加入你需要的析构过程

后记(2017-11-12):当然构造函数也要类似处理,否则 C++ 找不到实现类的定义,会产生编译错误。当然,这比运行时出现隐蔽的问题要容易发现。

最小二乘问题(二)

上一回我们介绍了最朴素的最小二乘问题。现在我们对它进行一些变形,来得到“新”的最小二乘问题。

前面的问题中,只有一项,如果我们关于相同的变量添加一个新项,得到:

我们引入个新的矩阵 $M = \begin{pmatrix}A\\C\end{pmatrix}$,把 $b$ 和 $d$ 连接成 $w=\begin{pmatrix}b\\d\end{pmatrix}$ 。那么上面的问题等价于

如果我们增加了含有新的变量的项呢?例如

我们同样构造一个新的分块矩阵 $M=\begin{pmatrix}A & 0 \\ 0 & C\end{pmatrix}$ ,将变量 $x$ 和 $y$ 也连接成一个向量 $v=\begin{pmatrix}x \\ y\end{pmatrix}$ ,类似地把 $b$ 和 $d$ 也连接成 $w=\begin{pmatrix}b\\d\end{pmatrix}$ 。那么这个问题等价于

除了上述的情形,我们还可能面对同时含有两个变量(或者多个)的项,比如

同样连接 $x$ $y$ 为 $v$,记矩阵 $M = \begin{pmatrix}A & B\end{pmatrix}$ ,于是它等价于 $\min_v \|Mv-c\|^2$ 。