线性最小二乘问题的数值解法

线性最小二乘问题的数值解法

这里我们讨论的问题主要是Ax=b，$b\notin R(A)$时的处理方法

最小二乘问题的描述

方程组的分类

对于方程组
$$Ax=b,A\in C^{m\times n},b\in C^m\text{\hspace{1em}(1)}$$

• 超定方程组：m>n
• 欠定方程组：m<n

线性最小二乘问题的基本概念

最小二乘问题的特征及一般表示

• 最小二乘问题

  设$A\in C^{m\times n},b\in C^m$，确定$x\in C^n$使得
  $$\Vert Ax-b{\Vert }_2=\underset{x\in C^n}{\min }\Vert Ax-b{\Vert }_2\text{\hspace{1em}(2)}$$

  问题(2)称为线性最小二乘问题（LS问题），而x称为最小二乘解。称$r(x)=b-Ax$为残差向量

• 所有最小二乘解的集合记为$S_{LS}$
  S L S = { x ∈ C n : x 满 足 ( 1 ) } (3) S_{LS}=\{x\in C^n:x满足(1)\}\tag{3} SLS​={x∈Cn:x满足(1)}(3)

  $S_{LS}$中2-范数最小者称为极小范数解，记为$x_{LS}$
  ∥ x L S ∥ 2 = min ⁡ { ∥ x ∥ 2 : x ∈ S L S } (4) \lVert x_{LS}\rVert_2=\min\{\lVert x\rVert_2:x\in S_{LS}\}\tag{4} ∥xLS​∥2​=min{∥x∥2​:x∈SLS​}(4)

• 极小范数解：方程组(1)有解的情况
  $$\Vert x_0{\Vert }_2=\underset{Ax=b}{\min }\Vert x{\Vert }_2\text{\hspace{1em}(5)}$$

• 极小范数最小二乘解：方程组无解的情况
  $$\Vert x_0{\Vert }_2=\underset{min\Vert Ax-b{\Vert }_2}{\min }\Vert x{\Vert }_2\text{\hspace{1em}(6)}$$

• 线性方程组(1)有解的充要条件是
  $$A{A}^{†}b=b\text{\hspace{1em}(7)}$$
  并且在有解时，其通解为
  $$\begin{gathered} x=A^{†}b+(I-A^{†}A)z \\ 其中z\in C^n任意\text{\hspace{1em}(8)} \end{gathered}$$
  证明：

  “必要性”：

  若$A{A}^{†}b=b$，则$x=A^{†}b$就是其解

  “充分性”

  若$Ax=b$有解，则 N ( A ) = { ( I n − A † A ) z : z ∈ C n } N(A)=\{(I_n-A^{\dagger}A)z:z\in C^n\} N(A)={(In​−A†A)z:z∈Cn}

• 如果方程组(1)有解，则它的极小范数解$x_0$唯一，并且$x_0=A^{†}b$

线性最小二乘问题的等价性问题

法方程

• x是最小二乘问题(2)的极小解的，即$x\in S_{LS}$的充分必要条件是x为方程

  $$A^{H}Ax=A^{H}b\text{\hspace{1em}(9)}$$

​ 的解，(9)式称为最小二乘问题的法方程。

• 证明：

  ${\min }_x\Vert Ax-b{\Vert }_2$等价于${\min }_x\Vert Ax-b{\Vert }_2^2={\min }_x(Ax-b{)}^T(Ax-b)$

  化简得， min ⁡ x { x T A T A x − x T A T b − b T A x + b T b } \min_x\{x^TA^TAx-x^TA^Tb-b^TAx+b^Tb\} minx​{xTATAx−xTATb−bTAx+bTb}。设$y=x^T{A}^{T}Ax-x^T{A}^{T}b-b^{T}Ax+b^{T}b$，
  令$\frac{dy}{dx}=0$得，
  $$\begin{gathered} 2A^{T}Ax-2A^{T}b=0 \\ {A}^{T}Ax=A^{T}b \end{gathered}$$
  所以，(2)式的极小解为$A^{H}Ax=A^{H}b$的解

KKT方程

• 设$A\in C^{m\times n}$，$b\in R^n$，则$x$和$r=b-Ax$分别是最小二乘问题(2)的极小解和残量的充分必要条件是x和r为鞍点系统
  $$\left[\begin{array}{cc}I& A\\ A^H& O\end{array}\right]\left[\begin{array}{c}r\\ x\end{array}\right]=\left[\begin{array}{c}b\\ 0\end{array}\right]\text{\hspace{1em}(10)}$$
  的解。上述线性系统称为最小二乘问题的KKT方程

• 证明：

  若$x$为最小二乘问题(2)的极小解，而$r=b-Ax$为其残量，则
  $$\begin{gathered} x=A^{†}b+(I-A^{†}A)z \\ r=b-Ax=b-A{A}^{†}b=(I-A{A}^{†})b \end{gathered}$$
  由广义逆$A^{†}$的性质$A{A}^{†}=(A{A}^{†}{)}^T,A{A}^{†}A=A$及等式$r+Ax=b$，得
  $$A^{H}r=A^H(I-A{A}^{†})b=[(I-A{A}^{†})A{]}^{H}b=[A-A{A}^{†}A{]}^{H}b=0$$
  所以(10)式是相容的线性系统，且x和r满足式(10).

  反之，设
  $$B=\left[\begin{array}{cc}I& A\\ A^H& O\end{array}\right]$$
  可以验证
  $$B^{†}=\left[\begin{array}{cc}I-A{A}^{†}& (A^{†}{)}^H\\ A^H& -A^{†}(A^{†}{)}^H\end{array}\right]$$
  所以
  $$\left[\begin{array}{c}r\\ x\end{array}\right]=B^{†}\left[\begin{array}{c}b\\ 0\end{array}\right]+(I-B^{†}B)\left[\begin{array}{c}y\\ z\end{array}\right]=\left[\begin{array}{c}(I-A{A}^{†})b\\ A^{†}b-(I-A^{†}A)z\end{array}\right]$$
  所以r,x分别是最小二乘问题的残量和极小解

线性最小二乘问题的正则化

设矩阵$A\in C_r^{m\times n}$的奇异值分解为
$$A=U\Sigma V^H$$
其中$U=[u_1,u_2,\cdots ,u_m]$为m阶酉矩阵；$V=[v_1,v_2,\cdots ,v_n]$为n阶酉矩阵；$\Sigma =\left[\begin{array}{cc}{\Sigma }_r& 0\\ 0& 0\end{array}\right]$

则最小二乘问题(2)的极小范数最小二乘解$x_{LS}$可表示为
$$x_{LS}=A^{†}b=\sum _{i=1}^r\frac{u_i^{H}b}{{\sigma }_i}{v}_i$$

最小二乘问题的求解

求解列满秩最小二乘问题的数值方法

法方程方法

分析

将最小二乘问题化为上文所述的法方程
$$A^{H}Ax=A^{H}b$$
由于$A\in C^{m\times n}$是列满秩矩阵，所以$rank(A^{H}A)=n$，且$A^{H}A$正定，所以可以用Cholesky分解求解方程

算法：

• 对n阶Hermite正定矩阵$B=A^{H}A$作Cholesky分解$B=L{L}^H$，其中L为下三角矩阵
• 依次解$Ly=A^{H}b$，$L^{H}x=y$得到最小二乘问题的解$x_{LS}$

QR分解法

分析

由于正交矩阵保持向量二范数不变，所以可以对最小二乘问题做变换
$$\begin{gathered} \underset{x}{\min }\Vert Ax-b\Vert \Rightarrow \underset{x}{\min }\Vert Q^T(Ax-b)\Vert \\ 其中Q为正交矩阵 \end{gathered}$$
所以可以找一个合适的Q使原来的最小二乘问题便于求解。

设A有QR分解：
$$\begin{gathered} A=Q\left[\begin{array}{c}R\\ 0\end{array}\right]=Q_{1}R \\ 其中Q\in R^{m\times n}为正交矩阵 \\ {Q}_1为Q的前n列组成的矩阵 \\ R为对角元均为正的上三角矩阵 \end{gathered}$$
取上式中的Q对b做变换得
$$f=Q^{T}b=\left[\begin{array}{c}{Q}_1^T\\ Q_2^T\end{array}\right]b=\left[\begin{array}{c}{f}_1\\ f_2\end{array}\right]$$
则
$$\begin{gathered} \Vert Q^T(Ax-b){\Vert }_2^2=\Vert \left[\begin{array}{c}R\\ 0\end{array}\right]x-Q^{T}b{\Vert }_2^2 \\ =\Vert \left[\begin{array}{c}R\\ 0\end{array}\right]x-\left[\begin{array}{c}{f}_1\\ f_2\end{array}\right]{\Vert }_2^2=\Vert Rx-f_1{\Vert }_2^2+\Vert f_2{\Vert }_2^2 \end{gathered}$$
由此可知，x是最小二乘问题的解，当且仅当x是上三角矩阵方程组$Rx=f_1$的解

在实施过程中，采取对$\hat{A}=[A,b]$进行QR分解
$$H_n\dots H_1[A,b]=[R,\hat{b}]\Rightarrow Q^T[A,b]=[R,\hat{b}]$$
得
$$A=QR,\hat{b}=Q^{T}b$$

算法

给定$A\in R^{m\times n}$为列满秩矩阵，$b\in R^m$

• 计算增广矩阵$\hat{A}=[A,b]$的QR分解
• 用回代法求解$Rx=f_1$

QR分解可以用

• Householder方法
• Givens正交化方法
• 修正的Gram-Schmit方法

求解秩亏最小二乘问题的数值解法

考虑 $\min \Vert Ax-b{\Vert }_2$，其中$A\in R_r^{m\ast n}(m\ge n)$,$rank(A)=r<n$。此时最小二乘问题有无穷个解。

列主元QR分解法

分析

• A是秩亏损时，QR分解不一定能给出列空间R(A)的一组标准正交基，但是可以通过列置换之后的矩阵$\hat{A}=AP$的QR分解来求解，即AP=QR，其中P是置换矩阵。

• 设$b=b_1+b_2,b_1\in R(A),b_2\in R^{\perp }(A)$
  可以证明：$Ax=b_1$
  假定$R(A)=R(Q_1)$，其中$Q_1\in R^{m\ast r}$且$Q_1^T{Q}_1=I_r$，即$Q_1$构成R(A)的一组标准正交基。

• 则存在矩阵$S\in R^{r\ast n}$（A在$Q_1$基下的坐标）和向量$h\in R^r$（$b_1$在$Q_1$基下的坐标），使得

  $$A=Q_{1}S,b_1=Q_{1}h$$
  $又\because Ax=b_1\therefore Q_{1}Sx=Q_{1}h\therefore Sx=h$

  所以问题转换为求解$Sx=h$，而方程组$Sx=h$总是有解的

• $$\begin{gathered} S=Q_1^{T}A, \\ h=Q_1^T{b}_1=Q_1^T(b-b_2)=Q_1^{T}b \\ {b}_2\in R^{\perp }(Q_1) \end{gathered}$$

• 现在问题变成求解R(A)的一组标准正交基$Q$

• 由于rank(A)=r < n,一般不能直接产生R(A)的一组标准正交基，但是如果对A进行适当的排列使其前r列线性无关，然后再进行QR分解，则仍然可产生R(A)的一组标准正交基

• 设有分解式

  $$AP=Q\left[\begin{array}{cc}{R}_{11}& R_{12}\\ 0& 0\end{array}\right]$$
  其中P为列置换矩阵，Q为正交矩阵，$R_{11}$为非奇异上三角阵，则Q的前r列就是R(A)的一组标准正交基

• 一旦求出分解式，就有
  $$(Q^{T}AP)(P^{T}x)=Q^T{b}_1$$
  令$P^{T}x=\left[\begin{array}{c}w\\ z\end{array}\right]$，得
  $$\begin{gathered} \left[\begin{array}{cc}{R}_{11}& R_{12}\\ 0& 0\end{array}\right]\left[\begin{array}{c}w\\ z\end{array}\right]=Q^T{b}_1=\left[\begin{array}{c}{Q}_1^T\\ Q_2^T\end{array}\right]b_1=\left[\begin{array}{c}h\\ g\end{array}\right] \\ \Rightarrow \left[\begin{array}{c}{R}_{11}w+R_{12}z\\ 0\end{array}\right]=\left[\begin{array}{c}h\\ g\end{array}\right] \\ \therefore w=R_{11}^{-1}(h-R_{12}z),g=0 \end{gathered}$$
  所以
  $$x=P\left[\begin{array}{c}w\\ z\end{array}\right]=P\left[\begin{array}{c}{R}_{11}^{-1}(h-R_{12}z)\\ z\end{array}\right]=P\left[\begin{array}{c}{R}_{11}^{-1}h\\ 0\end{array}\right]+\left[\begin{array}{c}-R_{11}^{-1}{R}_{12}\\ I_{n-r}\end{array}\right]z$$
  $x$就是最小二乘问题的通解
  $$令x_b=P\left[\begin{array}{c}{R}_{11}^{-1}h\\ 0\end{array}\right]$$
  $x_b$为最小二乘问题的基本解
  $$A{x}_b=\left[\begin{array}{c}{b}_1\\ 0\end{array}\right]$$

实际求解方法

• 综述：一边用Householder变换，一边做列变换将A中列向量二范数非零的向量提到下一个Householder变换列，直到A中没有二范数非零向量

• 下面描述求解的递推过程：

  考虑分解式(8)的具体实施过程。假设对对某一正整数k，已求k-1个Householder变换$H_1,H_2,\cdots ,H_{k-1}$和k-1个初等变换，$P_1,P_2,\cdots ,P_{k-1}$使得

  $R_{k-1}=H_{k-1}\cdots H_2{H}_{1}A{P}_1{P}_2\cdots P_{k-1}=\left[\begin{array}{cc}{R}_{11}^{(k-1)}& R_{12}^{(k-1)}\\ 0& R_{22}^{(k-1)}\end{array}\right]$，其中$R_{11}^{(k-1)} $为非奇异上三角矩阵，记$R_{22}^{(k-1)}=[V_k{(k-1)},V_{k+1}^{(k-1)}，\cdots ,V_n^{(k-1)}] $，其中$V_i^{{(k-1)}$表示$R_{22}}{(k-1)}$的第（i-k+1)列，下一步是确定指标P($k\le P \le n$)满足

  ∥ V p ( k − 1 ) ∥ 2 = max ⁡ { ∥ V k ( k − 1 ) ∥ 2 , ∥ V k + 1 ( k − 1 ) ∥ 2 , ⋯   , ∥ V n ( k − 1 ) ∥ 2 } \lVert V_p^{(k-1)}\rVert_2=\max \{\lVert V_k^{(k-1)}\rVert_2,\lVert V_{k+1}^{(k-1)}\rVert_2,\cdots,\lVert V_ n^{(k-1)}\rVert_2\} ∥Vp(k−1)​∥2​=max{∥Vk(k−1)​∥2​,∥Vk+1(k−1)​∥2​,⋯,∥Vn(k−1)​∥2​}

  若$\Vert V_p^{(k-1)}{\Vert }_2=0$，停止计算。否则取$P_k$为第k列与第p列交换的初等列变换矩阵，并确定一个Householder变换，${\hat{H}}_k\in R^{(m-k+1)\times (m-k+1)}$使得${\hat{H}}_k{V}_p^{(k-1)}={\gamma }_{kk}{e}_1$

  令$H_k=diag(I_{k-1},{\hat{H}}_k)$，则有$R_k=\left[\begin{array}{cc}{R}_{11}^{(k)}& R_{12}^{(k)}\\ 0& R_{22}^{(k)}\end{array}\right]$，其中$R_{11}^{(k-1)} $为非奇异上三角矩阵$R_k=H_k\begin{bmatrix}R_{11}^{(k-1)} & R_{12}^{(k-1)} \0 & R_{22}^{(k-1)}\end{bmatrix}P_k$

  $R_{12}$约化为0矩阵，$z_1,z_2,\cdots ,z_r$和置换，$\hat{P}$，使得$[R_{11},\cdots ,R_{1r}]z_1\cdots z_r\hat{P}=[T,0]$，T为上三角矩阵，令

  $Z=P{z}_1\cdots z_r$，则有$Q^{T}AZ=\left[\begin{array}{cc}T& 0\\ 0& 0\end{array}\right]$，由此可得$x_{LS}=z\left[\begin{array}{c}{T}^{-1}C\\ 0\end{array}\right]$，C为由$Q^{T}b$前r个分量构成的r维向量

算法

秩亏最小二乘问题的列主元QR分解法

• Step1：输入矩阵A，计算$P(j):=j，{\sigma }_j=A(:j{)}^{T}A(:j)$，置k=1
• Step2： 确定r使满足${\sigma }_r:{\max }_{k\le j\le n}{\sigma }_j$，若${\sigma }_r=0$停算
• Step3：交换$P(r)$与$P(k)$，${\sigma }_r$与${\sigma }_k$以及A(:,r)与A(:,k)
• Step4：确定一个Householder矩阵${\hat{H}}_k$使得${\hat{H}}_{r}A(k:m,k)={\gamma }_{kk}{e}_1$并计算$[A,b]:=diag(I_{k-1},{\hat{H}}_{k-1})[A,b],{\sigma }_j:={\sigma }_j-a_{kj}^2$
• Step5：置k:=k+1，转步骤2

变换技巧
$$\begin{gathered} \Vert Q^{T}x{\Vert }_2^2=\Vert [\alpha ,z{]}^T{\Vert }_2^2={\alpha }^2+\Vert z{\Vert }_2^2 \\ \Vert z{\Vert }_2^2=\Vert x{\Vert }_2^2-{\alpha }^2 \\ \therefore \Vert x^{(j)}{\Vert }_2^2=\Vert x^{(j-1)}{\Vert }_2^2-{\alpha }_{kj}^2 \end{gathered}$$

二、奇异值分解法

设$A\in R^{m\times n}(m>n)$的奇异值分解为$A=U\Sigma V^T$,$\sigma =\left[\begin{array}{cc}{\Sigma }_r& 0\\ 0& 0\end{array}\right]$，其中$U=[u_1,\cdots ,u_m],V=[v_1,\cdots ,v_n],{\Sigma }_r=diag({\sigma }_1,\cdots ,{\sigma }_r)$，则有
$$\begin{gathered} x_{LS}=A^{+}b=V\left[\begin{array}{cc}{\Sigma }_r^{-1}& 0\\ 0& 0\end{array}\right]U^T \\ =\sum _{i=1}^r\frac{u_i^{T}b}{{\sigma }_i}{v}_i \end{gathered}$$

参考文献

马昌凤 柯艺芬.数值线性代数与算法[M].国防工业出版社:北京,2017:247-279.