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

前言

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

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

回顾最小二乘的线性解

对于线性超定方程组$Ax=b,A\in R^{m\times n},m>n$，其最小二乘解可表示为：
$$\begin{gathered} \underset{x}{\mathrm{argmin}}||Ax-b|{|}_2^2=\underset{x}{\mathrm{argmin}}(Ax-b{)}^T(Ax-b) \\ \to A^{T}Ax=A^{T}b \\ ifrank(A)=n,x=(A^{T}A{)}^{-1}{A}^{T}b \end{gathered}$$

以上解法具有两个问题：①$A^{T}A$求逆运算的效率；②$rank(A)<n$时的最小二乘解。

首先讨论问题①。

列满秩矩阵的最小二乘解法

对于线性方程组$Ax=b,A\in R^{m\times n},m>n,rank(A)=n$，有$A^{T}A$对称且可逆。回顾矩阵分解的知识，可知$A^{T}A$矩阵能够进行Cholesky分解和QR分解。

Cholesky分解求线性最小二乘解

通过将$A^{T}A$分解为下三角矩阵$L$的乘积$L{L}^T$：
$$A^{T}A=L{L}^T$$
则可以将复杂的线性方程组求解：
$$A^{T}Ax=A^{T}b$$
转化为两个简单的三角线性方程组求解：
$$L{L}^{T}x=A^{T}b\to L^{T}x=y,Ly=A^{T}b$$

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

QR分解求线性最小二乘解

对于列满秩矩阵$A$而言，可以唯一分解为正交矩阵与对角元为正的上三角矩阵的乘积：
A = Q R = Q [ R 0 ] A=QR=Q

$$\begin{bmatrix} R \\ 0 \end{bmatrix}$$ A=QR=Q[R0​]
现在考虑最小二乘法的QR分解法：
$$\begin{gathered} \underset{x}{\mathrm{argmin}}||Ax-b|{|}_2^2=\underset{x}{\mathrm{argmin}}||Q\left[\begin{array}{c}R\\ 0\end{array}\right]x-b|| \\ =\underset{x}{\mathrm{argmin}}||\left[\begin{array}{c}R\\ 0\end{array}\right]x-Q^{T}b|{|}_2^2 \\ =\underset{x}{\mathrm{argmin}}||\left[\begin{array}{c}Rx\\ 0\end{array}\right]-\left[\begin{array}{c}{\beta }_1\\ {\beta }_2\end{array}\right]|{|}_2^2 \\ =\underset{x}{\mathrm{argmin}}||Rx-{\beta }_1|{|}_2^2+||{\beta }_2|{|}_2^2 \\ ⟺\underset{x}{\mathrm{argmin}}||Rx-{\beta }_1|{|}_2^2 \\ \to Rx={\beta }_1,x=R^{-1}{\beta }_1 \\ \left[\begin{array}{c}{\beta }_1\\ {\beta }_2\end{array}\right]=Q^{T}b \end{gathered}$$

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

亏秩矩阵的最小二乘解法

对于线性方程组$Ax=b,A\in R^{m\times n},m>n,rank(A)<n$，称为亏秩超定方程组。此时$A^{T}A$不是正定矩阵，QR分解也不是唯一的。通常使用SVD来解决亏秩问题。

SVD分解求亏秩最小二乘解

对于矩阵$A_{m\times n},m>n,rank(A)=r<n$，可分解为正交矩阵与奇异值矩阵的乘积：
A = U m × m [ Σ r × r 0 r × ( n − r ) 0 ( m − r ) × r 0 ( m − r ) × ( n − r ) ] V n × n T U T U = E , V T V = E A=U_{m\times m}

$$\begin{bmatrix}\Sigma_{r\times r} & 0_{r\times (n-r)} \\ 0_{(m-r)\times r} & 0_{(m-r)\times (n-r)}\end{bmatrix}$$ V_{n\times n}^T \\ \quad \\ U^TU=E,V^TV=E A=Um×m​[Σr×r​0(m−r)×r​​0r×(n−r)​0(m−r)×(n−r)​​]Vn×nT​UTU=E,VTV=E
现在考虑最小二乘法的SVD解法：
$$\begin{gathered} \underset{x}{\mathrm{argmin}}||Ax-b|{|}_2^2=\underset{x}{\mathrm{argmin}}||U^{T}Ax-U^{T}b|{|}_2^2 \\ =\underset{x}{\mathrm{argmin}}||U^{T}U\left[\begin{array}{cc}\Sigma & 0\\ 0& 0\end{array}\right]V^{T}x-U^{T}b|{|}_2^2 \\ =\underset{x}{\mathrm{argmin}}||\left[\begin{array}{cc}\Sigma & 0\\ 0& 0\end{array}\right]V^{T}x-U^{T}b|{|}_2^2 \\ =\underset{x}{\mathrm{argmin}}||\left[\begin{array}{cc}\Sigma & 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{\alpha }_1\\ {\alpha }_2\end{array}\right]-\left[\begin{array}{c}{\beta }_1\\ {\beta }_2\end{array}\right]|{|}_2^2 \\ =\underset{x}{\mathrm{argmin}}||\left[\begin{array}{c}\Sigma {\alpha }_1\\ 0\end{array}\right]-\left[\begin{array}{c}{\beta }_1\\ {\beta }_2\end{array}\right]|{|}_2^2 \\ =\underset{x}{\mathrm{argmin}}||\Sigma {\alpha }_1-{\beta }_1|{|}_2^2+||{\beta }_2|{|}_2^2 \\ \to \Sigma {\alpha }_1={\beta }_1 \\ \left[\begin{array}{c}{\alpha }_1\\ {\alpha }_2\end{array}\right]=V^{T}x,\left[\begin{array}{c}{\beta }_1\\ {\beta }_2\end{array}\right]=U^{T}b \\ {\alpha }_1={\Sigma }^{-1}{\beta }_1,x=V\left[\begin{array}{c}{\alpha }_1\\ {\alpha }_2\end{array}\right] \end{gathered}$$
在上面的解法中，由于${\alpha }_2$是一个任意向量，因此求出的$x$不是唯一的。即亏秩方程的最小二乘解不唯一。可以选择范数最小的最小二乘解作为亏秩方程的解，即${\alpha }_2=\vec{0}$。

可以通过正交矩阵分块，将SVD最小二乘解进一步化简：
U = [ U 1 , U 2 ] , V = [ V 1 , V 2 ] U 1 ∈ R m × r , V 1 ∈ R n × r α 1 = Σ − 1 β 1 [ α 1 α 2 ] = [ V 1 T V 2 T ] x [ β 1 β 2 ] = [ U 1 T U 2 T ] b V 1 T x = Σ − 1 U 1 T b x = V 1 Σ − 1 U 1 T b U=[U_1,U_2],V=[V_1,V_2] \\ U_1\in R^{m\times r},V_1\in R^{n\times r} \\ \quad \\ \alpha_1=\Sigma^{-1}\beta_1 \\ \quad \\

$$\begin{bmatrix} \alpha_1 \\ \alpha_2 \end{bmatrix}$$ = $$\begin{bmatrix} V_1^T \\ V_2^T \end{bmatrix}$$ x \\ \quad \\ $$\begin{bmatrix} \beta_1 \\ \beta_2 \end{bmatrix}$$ = $$\begin{bmatrix} U_1^T \\ U_2^T \end{bmatrix}$$ b \\ \quad \\ V_1^Tx=\Sigma^{-1} U_1^Tb \\ \quad \\ x=V_1\Sigma^{-1} U_1^Tb \\ U=[U1​,U2​],V=[V1​,V2​]U1​∈Rm×r,V1​∈Rn×rα1​=Σ−1β1​[α1​α2​​]=[V1T​V2T​​]x[β1​β2​​]=[U1T​U2T​​]bV1T​x=Σ−1U1T​bx=V1​Σ−1U1T​b

需要说明的是，SVD也适用于列满秩矩阵的最小二乘求解。实际上SVD分解是最常用的线性最小二乘解法之一。

补充1：超定齐次方程组的线性最小二乘解法

之前的最小二乘法都在考虑$Ax=b$的最小非齐次方程组问题。现在补充齐次方程组的最小二乘解法。

对于齐次方程组$A_{m\times n}x=0,m>n$而言，其最小二乘解就是$A$的SVD分解后的$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​
进一步，超定齐次线性方程组$Ax=0$的最小二乘解还等于$A^{T}A$的最小特征值对应的特征向量：

证明：
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\times n}x=b,m<n,rank(A)=r$，其最小二乘也是通过SVD求解：
A = [ U 1 m × r , U 2 m × ( m − r ) ] [ Σ r × r 0 0 0 ( m − r ) × ( n − r ) ] [ V 1 n × r , V 2 n × ( n − r ) ] T A x = U 1 Σ V 1 T x = b y = V 1 T x , U 1 Σ y = b 特 别 的 ， 如 果 r a n k ( A ) = m ， 有 U Σ y = b , y = Σ − 1 U T b , x = V Σ − 1 U T b A=[U_1^{m\times r},U_2^{m\times (m-r)}]

$$\begin{bmatrix} \Sigma^{r\times r} & 0 \\ 0 & 0^{(m-r)\times (n-r)} \end{bmatrix}$$ [V_1^{n\times r},V_2^{n\times (n-r)}]^T \\ \quad \\ Ax = U_1\Sigma V_1^Tx=b \\ y=V_1^Tx,U_1\Sigma y=b \\ \quad \\ 特别的，如果rank(A)=m，有 \\ \quad \\ U\Sigma y=b,y=\Sigma^{-1}U^Tb,x=V\Sigma^{-1}U^Tb A=[U1m×r​,U2m×(m−r)​][Σr×r0​00(m−r)×(n−r)​][V1n×r​,V2n×(n−r)​]TAx=U1​ΣV1T​x=by=V1T​x,U1​Σy=b特别的，如果rank(A)=m，有UΣy=b,y=Σ−1UTb,x=VΣ−1UTb

后记

本篇介绍了最小二乘的矩阵解法，包括列满秩超定方程组的Cholesky分解法和QR分解法，以及亏秩超定方程组的SVD解法。

后续将继续学习非线性方程组的最小二乘问题。