介绍一下如何求解最小二乘法问题。

线性最小二乘

介绍一下如何求解最小二乘法问题。于是查看了一下Numpy中是如何求解最小二乘法问题的，也即numpy.polyfit函数的代码，一路追下来之后，发现最终使用了LAPACK的ZGELSD方法。它们之间的调用顺序如下：

1
numpy.polyfit -> numpy.lstsq -> numpy.linalg.lapack_lite.zgelsd

而实际上，在LAPACK中，存在四种求取实数最小二乘法的方式：

  • 使用 QR分解 或 LR分解 ： ZGELS
  • 使用 完全正交分解 : ZGELSY
  • 使用 SVD : ZGELSS
  • 使用 分治SVD : ZGELSD

问题描述

首先让我们来简单描述一下问题。我们需要解的是这样一个方程：

$$
Ax = b
$$

在这个式子中，我们可以认为$A$是一个$m\times n$的矩阵，$x$是一个$n\times 1$的矩阵，那么$b$就是一个 $m\times 1$的矩阵了。

我们期望根据给出的一组A和b，求出相应的x。

  • 对于$m = n$的情况。方程有唯一的解，也就是 $x = inv(A) * b$。（要求A非奇异）
  • 对于$m > n$的情况。称为超定问题。此时不存在解。
  • 对于$m < n$的情况，称为负定问题

我们现在所要面临的就是一个超定问题，由于原方程可以理解为用A的列向量来线性地表示b，而对于超定问题，b根本不在A的列向量张成的线性空间中，则这个方程根本没有办法获得精确解，在这时我们就需要求得最小二乘解，也即求使得误差$ J = \sum_{i=1}^N{(Ax_i - b)^2} $ 最小的$x$。因为2-范数$ |x|2 = \sqrt{\sum{i=1}^N{x_i^2} }$，所以我们的求解目标常常写作：

$$
min_x|Ax-b|_2^2
$$

其实也就是解上面这个式子，得到的结果称为最小二乘解。

最主要的问题就在于解决不同维度矩阵相乘除的问题。

基本解法

最初的想法肯定是左右同时乘除一个矩阵使得某一边成为方阵可以进行求逆处理：

$$
A^HAx = A^Hb \
x = inv(A^HA)*A^Hb
$$

我们可以照着这个用Python写一下代码：

1
2
3
4
5
6
7
8
import numpy as np

A = np.mat("[1,1; 2,0; 0,2]")
x = np.mat("[2;3]")
b = A*x

xp = (A.transpose() * A).I * A.transpose() * b
print(xp - x)
1
2
3
4
[[5 1]
    [1 5]]
[[4.4408921e-16]
    [0.0000000e+00]]

上面代码的结果应该是两个小于1e-15的数字，这样说明在这个程序中我们的结果还算比较好。

但是这种直接通过$A^HA$的方式求解问题的解法有很大的不稳定性。（我们可以随便换换A和x，就能看到结果的波动了）

ZGELS

这个方法是指使用QR分解求最小二乘解。

对上面问题中的A进行QR分解：

$$
A = QR
$$

其中Q是一个半正定矩阵，满足 $ Q^HQ = I $，其大小应为 $m\times n$。R是一个 $n \times n$ 的上三角矩阵所以我们可以把原本的式子写成这样:

$$
Q^H A x = R x = Q^Hb
$$

然后：

$$
x = inv(R) * Q^H *b
$$

简单的Python代码：

1
2
q, r = np.linalg.qr(A)
print(r.I * q.T *b - x) 
1
2
[[-4.4408921e-16]
 [-4.4408921e-16]]

ZGELSS

所谓SVD指下面这种分解，将矩阵分成：

$$
A = U \Sigma V^H
$$

其中$ U, V $分别是$m \times m,n \times n$ 的酉矩阵， 而$ \Sigma $是$m\times n$的除了对角线元素外全为0的矩阵。

为了简便起见我们把矩阵表示为：

$$
\Sigma = \left[\begin{matrix} S \ 0 \end{matrix} \right] \
U = \left[ \begin{matrix} U_1, U_2 \end{matrix} \right]
$$

其中S为 $n \times n$的对角矩阵, $U_1$为 $m\times n$ 的矩阵，$U_2$ 为 $ m \times (m-n)$的矩阵。

如此我们可以列出这样的式子：

$$
\begin{aligned}
J &= |Ax - b|_2^2 = \left|U \left[\begin{matrix} S \ 0 \end{matrix} \right] V^Hx - b \right|_2^2 \
&= \left|\left[\begin{matrix} S \ 0 \end{matrix} \right] V^Hx - U^H b \right|_2^2 \
&= \left|\left[\begin{matrix} SV^Hx \ 0 \end{matrix} \right] - \left[\begin{matrix} U_1^Hb \ U_2^Hb \end{matrix} \right] \right |_2^2 \
&= \left|\left[\begin{matrix} SV^Hx - U_1^Hb \ -U_2^Hb \end{matrix} \right] \right |_2^2 \
&= \left|SV^Hx - U_1^Hb\right|_2^2 + \left| -U_2^Hb\right|_2^2 \ge \left| U_2^Hb\right|_2^2
\end{aligned}
$$

因为2-范数大于等于0，由上式可以得到误差$J$最小只能为$ \left| U_2^Hb\right|_2^2 $，且仅在：$ \left|SV^Hx - U_1^Hb\right|_2^2 = 0$，也即$ SV^Hx = U_1^Hb$ 时成立。

可以根据上面求出$x$为：

$$
x = \left(SV^H\right)^{-1}U_1^Hb
$$

这篇博客的MathJax花了我比较多时间来Debug…最后发现少写了一个下划线