数值计算之 最小二乘法（3）最小二乘的矩阵解法

前言

之前将最小二乘法与线性方程组求解关联，得到了线性最小二乘的矩阵形式，以及线性最小二乘的几何意义。

本篇将介绍线性最小二乘的矩阵解法。

回顾最小二乘的线性解

对于线性超定方程组 A x = b , A ∈ R m × n , m > n Ax=b,A\in R^{m\times n},m>n Ax=b,A∈Rm×n,m>n，其最小二乘解可表示为： arg min ⁡ x ∣ ∣ A x − b ∣ ∣ 2 2 = arg min ⁡ x ( A x − b ) T ( A x − b ) → A T A x = A T b i f r a n k ( A ) = n , x = ( A T A ) − 1 A T b \argmin_x ||Ax-b||^2_2=\argmin_x (Ax-b)^T(Ax-b) \\ \to A^TAx=A^Tb \\ \quad \\ if \quad rank(A)=n, x=(A^TA)^{-1}A^Tb xargmin​∣∣Ax−b∣∣22​=xargmin​(Ax−b)T(Ax−b)→ATAx=ATbifrank(A)=n,x=(ATA)−1ATb

以上解法具有两个问题：① A T A A^TA ATA求逆运算的效率；② r a n k ( A ) < n rank(A)n,rank(A)=n Ax=b,A∈Rm×n,m>n,rank(A)=n，有 A T A A^TA ATA对称且可逆。回顾矩阵分解的知识，可知 A T A A^TA ATA矩阵能够进行Cholesky分解和QR分解。

Cholesky分解求线性最小二乘解

通过将 A T A A^TA ATA分解为下三角矩阵 L L L的乘积 L L T LL^T LLT： A T A = L L T A^TA=LL^T ATA=LLT 则可以将复杂的线性方程组求解： A T A x = A T b A^TAx=A^Tb ATAx=ATb 转化为两个简单的三角线性方程组求解： L L T x = A T b → L T x = y , L y = A T b LL^Tx=A^Tb \to L^Tx=y, Ly=A^Tb LLTx=ATb→LTx=y,Ly=ATb

Cholesky分解速度很快，但精度一般，稳定性差。适合在限定时间内的大规模超定线性方程计算求解。

QR分解求线性最小二乘解

对于列满秩矩阵 A A A而言，可以唯一分解为正交矩阵与对角元为正的上三角矩阵的乘积： A = Q R = Q [ R 0 ] A=QR=Q\begin{bmatrix} R \\ 0 \end{bmatrix} A=QR=Q[R0​] 现在考虑最小二乘法的QR分解法： arg min ⁡ x ∣ ∣ A x − b ∣ ∣ 2 2 = arg min ⁡ x ∣ ∣ Q [ R 0 ] x − b ∣ ∣ = arg min ⁡ x ∣ ∣ [ R 0 ] x − Q T b ∣ ∣ 2 2 = arg min ⁡ x ∣ ∣ [ R x 0 ] − [ β 1 β 2 ] ∣ ∣ 2 2 = arg min ⁡ x ∣ ∣ R x − β 1 ∣ ∣ 2 2 + ∣ ∣ β 2 ∣ ∣ 2 2    ⟺    arg min ⁡ x ∣ ∣ R x − β 1 ∣ ∣ 2 2 → R x = β 1 , x = R − 1 β 1 [ β 1 β 2 ] = Q T b \argmin_x ||Ax-b||^2_2=\argmin_x ||Q\begin{bmatrix} R \\ 0 \end{bmatrix}x-b|| \\ = \argmin_x ||\begin{bmatrix} R \\ 0 \end{bmatrix}x-Q^Tb||^2_2 \\ = \argmin_x ||\begin{bmatrix} Rx \\ 0 \end{bmatrix} - \begin{bmatrix} \beta_1 \\ \beta_2 \end{bmatrix}||^2_2 \\ = \argmin_x ||Rx-\beta_1||^2_2 + ||\beta_2||^2_2 \\ \iff \argmin_x ||Rx-\beta_1||^2_2 \\ \to Rx=\beta_1,x=R^{-1}\beta_1 \\ \quad \\ \begin{bmatrix} \beta_1 \\ \beta_2 \end{bmatrix} = Q^Tb xargmin​∣∣Ax−b∣∣22​=xargmin​∣∣Q[R0​]x−b∣∣=xargmin​∣∣[R0​]x−QTb∣∣22​=xargmin​∣∣[Rx0​]−[β1​β2​​]∣∣22​=xargmin​∣∣Rx−β1​∣∣22​+∣∣β2​∣∣22​⟺xargmin​∣∣Rx−β1​∣∣22​→Rx=β1​,x=R−1β1​[β1​β2​​]=QTb

QR分解的速度较快，精度一般。不过由于存在高精度的QR分解方式，因此适合需要高精度解的小规模超定方程组计算。

亏秩矩阵的最小二乘解法

对于线性方程组 A x = b , A ∈ R m × n , m > n , r a n k ( A ) < n Ax=b,A\in R^{m\times n},m>n,rank(A)n,rank(A) n , r a n k ( A ) = r < n A_{m\times n},m>n,rank(A)=rn,rank(A)=r n A_{m\times n}x=0,m>n Am×n​x=0,m>n而言，其最小二乘解就是 A A A的SVD分解后的 V V V的最后一个列向量。

证明： A = U [ Σ 0 ] V T A x = U [ Σ 0 ] V T x = 0 y = V T x , Σ y = 0 也 就 是 说 ， 求 A x = 0 ， 转 换 为 求 Σ y = 0 的 最 小 二 乘 解 arg min ⁡ y ∣ ∣ Σ y ∣ ∣ 2 2 = arg min ⁡ y y T Σ T Σ y = arg min ⁡ y ∑ i = 1 n σ i 2 y i 2 s . t . ∣ ∣ y ∣ ∣ > 0 Σ 对 角 线 元 素 由 大 到 小 排 列 ， 则 满 足 y = [ 0 , 0 , … , 0 , y n ] T , y n ≠ 0 时 ， 获 得 Σ y 的 最 小 二 乘 解 V T x = y , x = V y = y n v n i f ∣ ∣ y ∣ ∣ 2 = 1 t h e n x = v n A=U\begin{bmatrix} \Sigma \\ 0 \end{bmatrix} V^T \\ \quad \\ Ax= U\begin{bmatrix} \Sigma \\ 0 \end{bmatrix} V^Tx=0 \\ \quad \\ y=V^Tx,\Sigma y=0 \\ \quad \\ 也就是说，求Ax=0，转换为求\Sigma y=0的最小二乘解 \\ \quad \\ \argmin_y ||\Sigma y||^2_2 \\ =\argmin_y y^T\Sigma^T\Sigma y \\ =\argmin_y \sum_{i=1}^n \sigma_i^2 y_i^2 \\ s.t. \quad ||y||>0 \\ \quad \\ \Sigma 对角线元素由大到小排列，则满足 \\ \quad \\ y=[0,0,\dots,0,y_n]^T,y_n\ne 0 \\ \quad \\ 时，获得\Sigma y的最小二乘解 \\ \quad \\ V^Tx=y,x=Vy=y_nv_n \\ \quad \\ if \quad ||y||_2=1 \\ then \quad x=v_n \\ A=U[Σ0​]VTAx=U[Σ0​]VTx=0y=VTx,Σy=0也就是说，求Ax=0，转换为求Σy=0的最小二乘解yargmin​∣∣Σy∣∣22​=yargmin​yTΣTΣy=yargmin​i=1∑n​σi2​yi2​s.t.∣∣y∣∣>0Σ对角线元素由大到小排列，则满足y=[0,0,…,0,yn​]T,yn​​=0时，获得Σy的最小二乘解VTx=y,x=Vy=yn​vn​if∣∣y∣∣2​=1thenx=vn​ 进一步，超定齐次线性方程组 A x = 0 Ax=0 Ax=0的最小二乘解还等于 A T A A^TA ATA的最小特征值对应的特征向量：

证明： A = U [ Σ 0 ] V T A T A = V [ Σ 0 ] [ Σ 0 ] V T = V Σ 2 V T A T A x = V Σ 2 V T x = 0 y = V T x , Λ = Σ 2 , Λ y = 0 A=U\begin{bmatrix} \Sigma \\ 0 \end{bmatrix} V^T \\ \quad \\ A^TA = V \begin{bmatrix} \Sigma & 0 \end{bmatrix}\begin{bmatrix} \Sigma \\ 0 \end{bmatrix}V^T=V\Sigma^2V^T \\ \quad \\ A^TAx = V\Sigma^2V^Tx=0 \\ y=V^Tx,\Lambda=\Sigma^2,\Lambda y=0 \\ \quad \\ A=U[Σ0​]VTATA=V[Σ​0​][Σ0​]VT=VΣ2VTATAx=VΣ2VTx=0y=VTx,Λ=Σ2,Λy=0 后面的证明就与上面一模一样了。

补充2：欠定行满秩方程组的线性最小二乘解法

还需要补充的是欠定方程组的最小二乘解法。虽然用的少，但却实实在在的碰到了这个问题：对极几何中本质矩阵的求解。

首先需要确定，线性齐次欠定方程组必然存在无穷组解析解。

对于非齐次欠定方程组 A m × n x = b , m < n , r a n k ( A ) = r A_{m\times n}x=b,m