【代数之美】线性方程组Ax=0的求解方法

在3D视觉中，我们常常会遇到这样一个问题：求解线性方程组 A x = 0 Ax=0 Ax=0，从矩阵映射的角度来说，所有解组成了矩阵 A A A的零空间。一个典型的场景比如用八点法求解本质矩阵 E E E，参见我前面的博文：立体视觉入门指南（2）：关键矩阵（本质矩阵，基础矩阵，单应矩阵）。这是一个基础且常见的线性代数问题，本篇我们来讨论下此类问题的解法，也算是一个入门课程。

对于 A x = 0 Ax=0 Ax=0，很明显， x = 0 x=0 x=0必然是其中一个解，但是对我们实际应用来说，得到一个零解往往没有什么意义，所以不做讨论。

如前所述，我们实际要得到的，是非零解，它们显然比零解要隐藏的深。

在讨论非零解之前，我们必须介绍一下特解的概念。一个显而易见的点是， A x = 0 Ax=0 Ax=0是尺度不变的，即若 x s x_s xs​是其中一个解，则 k s x s k_sx_s ks​xs​必然也是其解（ k s k_s ks​是实数），而 k s x s k_sx_s ks​xs​和 x s x_s xs​是线性相关的；再进一步，如果 x t x_t xt​是另一个与 x s x_s xs​线性无关的解，即 A x s = 0 , A x t = 0 , x s 和 x t 线 性 无 关 Ax_s=0,Ax_t=0,x_s和x_t线性无关 Axs​=0,Axt​=0,xs​和xt​线性无关

则显然 A ( k s x s + k t x t ) = 0 A(k_sx_s+k_tx_t)=0 A(ks​xs​+kt​xt​)=0也是成立的（ k s , k t k_s,k_t ks​,kt​是实数），即 k s x s + k t x t k_sx_s+k_tx_t ks​xs​+kt​xt​也是方程的解。这揭示一个现象：如果方程组 A x = 0 Ax=0 Ax=0有若干个线性无关解 x 1 , x 2 , . . . , x n x_1,x_2,...,x_n x1​,x2​,...,xn​，则 x 1 , x 2 , . . . , x n x_1,x_2,...,x_n x1​,x2​,...,xn​的任意线性组合也是其解，换句话说，所有的解都可以通过 x 1 , x 2 , . . . , x n x_1,x_2,...,x_n x1​,x2​,...,xn​来线性组合。 我们便把这些线性无关解 x 1 , x 2 , . . . , x n x_1,x_2,...,x_n x1​,x2​,...,xn​称为方程组 A x = 0 Ax=0 Ax=0的特解。

所以，我们的目的，实际上是要计算所有的非零特解。

那么首先我们想搞清楚，到底有多少非零特解呢？

我们从消元回代解法来分析，假设一个矩阵 A ∈ R m × n ( m = 3 , n = 4 ) A\in R^{m\times n}(m=3,n=4) A∈Rm×n(m=3,n=4)： A = [ 1 3 6 6 2 2 4 8 4 4 8 16 ] A=\left[\begin{matrix}1&3&6&6\\2&2&4&8\\4&4&8&16\end{matrix}\right] A=⎣⎡​124​324​648​6816​⎦⎤​

首先对 A A A进行消元， [ 1 3 6 6 2 2 4 8 4 4 8 16 ] → [ 1 3 6 6 2 2 4 8 0 0 0 0 ] → [ 1 3 6 6 0 2 4 2 0 0 0 0 ] \left[\begin{matrix}1&3&6&6\\2&2&4&8\\4&4&8&16\end{matrix}\right]\rightarrow\left[\begin{matrix}1&3&6&6\\2&2&4&8\\0&0&0&0\end{matrix}\right]\rightarrow\left[\begin{matrix}\boxed1&3&6&6\\0&\boxed2&4&2\\0&0&0&0\end{matrix}\right] ⎣⎡​124​324​648​6816​⎦⎤​→⎣⎡​120​320​640​680​⎦⎤​→⎣⎡​1​00​32​0​640​620​⎦⎤​

在消元后的矩阵中，方框内的1和2称为主变量（pivot variable），也称为主元。主元的个数称为矩阵的秩，这里主元的个数为2，那么矩阵A的秩（rank）也为2，即 r a n k ( A ) = 2 rank(A)=2 rank(A)=2。

主元所在的列称为主列（pivot column），这里的主列分别为第1列和第2列，其余列第3列和第4列为自由列（free column）。自由列中的变量为自由变量（free variable），自由变量的个数为 n − r a n k ( A ) = 4 − 2 = 2 n-rank(A)=4-2=2 n−rank(A)=4−2=2。

按照消元解法，在消元后我们给自由变量赋值，回代去求主列变量的值。分别给自由变量 [ x 3 x 4 ] \left[\begin{matrix}x_3\\x_4\end{matrix}\right] [x3​x4​​]赋值为 [ 1 0 ] \left[\begin{matrix}1\\0\end{matrix}\right] [10​]和 [ 0 1 ] \left[\begin{matrix}0\\1\end{matrix}\right] [01​]，回代方程得到以下两个解向量： [ 0 − 2 1 0 ] , [ − 3 − 1 0 1 ] \left[\begin{matrix}0\\-2\\1\\0\end{matrix}\right],\left[\begin{matrix}-3\\-1\\0\\1\end{matrix}\right] ⎣⎢⎢⎡​0−210​⎦⎥⎥⎤​,⎣⎢⎢⎡​−3−101​⎦⎥⎥⎤​

这两个解向量就是线性方程组 A x = 0 Ax=0 Ax=0的两个非零特解，其他解都可以通过这两个解来线性表达： x = a [ 0 − 2 1 0 ] + b [ − 3 − 1 0 1 ] x=a\left[\begin{matrix}0\\-2\\1\\0\end{matrix}\right]+b\left[\begin{matrix}-3\\-1\\0\\1\end{matrix}\right] x=a⎣⎢⎢⎡​0−210​⎦⎥⎥⎤​+b⎣⎢⎢⎡​−3−101​⎦⎥⎥⎤​。也可以说特解的任意线性组合构造出矩阵 A A A的整个零空间。

可知，线性方程组 A x = 0 Ax=0 Ax=0的非零特解个数为 n − r a n k ( A ) n-rank(A) n−rank(A)。

以上只是从纯数学角度来分析，但在计算机科学应用中，我们最常遇到的情况是，通过大量具有一定噪声的观测值，构造方程数远大于未知量的超定线性方程组 A x = 0 , A ∈ R m × n , m ≫ n Ax=0,A\in R^{m\times n},m\gg n Ax=0,A∈Rm×n,m≫n，这时候往往没有严格意义上的解。

比如在求解两张图像的本质矩阵 E E E时，我们会匹配大量的特征点对，然后基于对极约束构造线性方程组 A e = 0 Ae=0 Ae=0求解本质矩阵的元素向量 e e e，特征点对的个数往往远超过本质矩阵的元素数（即 m ≫ n m\gg n m≫n），且因为点位噪声的存在，对极约束不是严格成立的（即 A x ≈ 0 Ax\approx 0 Ax≈0）。面对这种情况，我们通常的做法是求最小二乘解。即： x ^ = arg ⁡ min ⁡ x ∣ ∣ A x ∣ ∣ 2 2 \hat {x}=\arg\min_x||Ax||_2^2 x^=argxmin​∣∣Ax∣∣22​

但我们立马发现一个问题，如前所述的尺度不变性，如果 x s x_s xs​是一个解，则 k s x s k_sx_s ks​xs​也是解，则当我们求出一个解后，无限减小解向量的模长， ∣ ∣ A x ∣ ∣ 2 2 ||Ax||_2^2 ∣∣Ax∣∣22​不就会无限小吗，到底该取哪个呢？这显然这陷入一个僵局。

所以，我们要给 x x x一个约束，即限定它的模长为一个常量，比如最常用的单位模长1，即 x x x必须满足， ∣ ∣ x ∣ ∣ 2 2 = 1 ||x||_2^2=1 ∣∣x∣∣22​=1。这就重新构建了一个带约束的最小二乘问题： x ^ = arg ⁡ min ⁡ x ∣ ∣ A x ∣ ∣ 2 2 , subject to ∣ ∣ x ∣ ∣ 2 2 = 1 \hat {x}=\arg\min_x||Ax||_2^2,\text{subject to} ||x||_2^2=1 x^=argxmin​∣∣Ax∣∣22​,subject to∣∣x∣∣22​=1

根据拉格朗日乘数法，另 L ( x , λ ) = ∣ ∣ A x ∣ ∣ 2 2 + λ ( 1 − ∣ ∣ x ∣ ∣ 2 2 ) L(x,\lambda)=||Ax||_2^2+\lambda(1-||x||_2^2) L(x,λ)=∣∣Ax∣∣22​+λ(1−∣∣x∣∣22​)

要求 L ( x , λ ) L(x,\lambda) L(x,λ)最小值，则首先求极值，即一阶偏导为0的点，我们分别对 x , λ x,\lambda x,λ求偏导并另其等于0： L ′ ( x ) = 2 A T A x − 2 λ x = 0 L ′ ( λ ) = 1 − x T x = 0 \begin{aligned} L^{'}(x)&=2A^TAx-2\lambda x=0\\ L^{'}(\lambda)&=1-x^Tx=0 \end{aligned} L′(x)L′(λ)​=2ATAx−2λx=0=1−xTx=0​

则有 A T A x = λ x A^TAx=\lambda x ATAx=λx，这熟悉的公式告诉我们， x x x是矩阵 A T A A^TA ATA特征值为 λ \lambda λ的特征向量。

但是矩阵 A T A A^TA ATA最多有 n n n个特征值，对应着 n n n个特征向量，哪一个是最小值点呢？

我们再作一下推导：当 x x x是矩阵 A T A A^TA ATA特征值为 λ \lambda λ的特征向量时，有 ∣ ∣ A x ∣ ∣ 2 2 = x T A T A x = x T λ x = λ x T x = λ ||Ax||_2^2=x^TA^TAx=x^T\lambda x=\lambda x^Tx=\lambda ∣∣Ax∣∣22​=xTATAx=xTλx=λxTx=λ

因此，当 λ \lambda λ取最小时， ∣ ∣ A x ∣ ∣ 2 2 ||Ax||_2^2 ∣∣Ax∣∣22​取最小值，即最小二乘解 x ^ \hat{x} x^是矩阵 A T A A^TA ATA最小特征值对应的特征向量。这就是超定线性方程组 A x = 0 Ax=0 Ax=0的最小二乘解法。