正交化与最小二乘法

隐藏目录

Householder反射

Jacobi-Givens旋转

$QR$ 分解

Householder $QR$ 分解
Givens $QR$ 方法
Gram-Schmidt方法

满秩的最小二乘问题

法方程组方法
用 $QR$ 分解求解LS问题

病态矩阵的正交分解

选主列的 $QR$ 分解
完全正交分解

双对角化

秩亏损的LS问题

完全正交分解,SVD与 $x_{\text{LS}}$
选主列的 $QR$ 分解与基本解

欠定方程组的LS解

广义逆
$LU$ 分解
正交分解法
消去法
消去与投影法

秩亏损情形

广义最小二乘问题

本章主要讨论超定方程组与不定方程组的最小二乘解,即 $\|Ax-b\|_2$ 的最小化,其中 $A\in \mathbb{R}^{n\times n}$ .解此问题的可靠方法是通过正交变换将 $A$ 约化为各种典型形式.我们在以下的讨论中主要考虑实矩阵.但容易看出,复矩阵的 $QR$ 分解完全可以仿照实矩阵进行.

Householder反射

通过正交变换将矩阵约化为典型形式 (如Hessenberg形式)的典型方法是通过一系列正交变换,逐个地将矩阵列向量的某些指定分量化为 $0$ .下面我们就先讨论这一过程的两种实现:Householder反射与Jacobi-Givens旋转.

Householder反射的几何意义是取向量 $x$ 关于指定超平面 $S\equiv\mathrm{span}\{v\},\quad v\neq 0$ 的反射. $v^Tx/||v||$ 为向量 $x$ 关于平面 $S$ 的法线的投影,从而 $x-2\frac{v^Tx}{||v||^2}v=x-2\frac{vv^T}{v^Tv}x$ 为 $x$ 关于该平面的投影.因此,设 $v\in \mathbb{R}^n$ 是非零向量,称形如 $P=I-\frac{2}{v^Tv}vv^T$ 的 $n\times n$ 矩阵 $P$ 为Householder反射.如果用 $P$ 去乘向量 $x$ ,则得到 $Px$ 是 $x$ 关于超平面 $\mathrm{span}\{v\}^{\bot}$ 的反射.

利用Householder反射可将一个向量的某些选定分量变为 $0$ .设给定 $0\neq x\in \mathbb{R}^n$ ,欲使 $Px$ 沿 $e_1$ 的方向.由几何直观容易知道 $\begin{equation*} v=x\pm \|x\|_2e_1 \Longrightarrow Px= \left(I-\frac{2vv^T}{v^Tv}\right)x=\mp \|x\|_2e_1. \end{equation*}$ 有一些细节需要注意:

由于 $v_1=x-\|x\|_2e_1$ 使得 $Px$ 是 $e_1$ 的正倍数,故常常采用.但当 $x$ 接近 $e_1$ 的一个正倍数时,则会出现严重的相消.由 $v_1=x_1-\|x\|_2=\frac{x_1^2-\|x\|_2^2}{x_1+\|x\|_2}=\frac{- (x_2^2+\cdots +x_n^2)}{x_1+\|x\|_2}$ 在 $x_1>0$ 的情况下则可避免相消.
实际中常常将Householder向量规范化使得 $v(1)=1$ ,这就允许将 $v(2:n)$ 存储到 $x$ 已化零的位置,即 $x(2:n)$ ,而无需增加额外的存储开销.在这种情况下,我们将 $v(2:n)$ 称为Householder向量的基本部分.
在进行Householder矩阵与向量的乘法运算时,我们并不直接采用矩阵-向量乘法,而是采用 $Px=x-2v(v^Tx)/(v^Tv)$ 的运算顺序,可以显著减少运算量.
在基于Householder反射的矩阵分解算法中,常用到形如 $Q=Q_1Q_2\cdots Q_r$ , $Q_j=I-\beta_jv^{(j)}v^{(j)T}$ 的若干个Householder矩阵的乘积,其中 $v^{(j)}=\left(\underbrace{0,\cdots,0 }_{j-1\text{个}},1,v_{j+1}^{(j)},\cdots,v_n^{(j)}\right)^T$ .通常我们并不显式地求出 $Q$ ,而是将其保持因子形式,逐个作用到给定矩阵上去,在不需显式求出变换矩阵时可显著减少运算量.
若要将 $Q$ 显式求出(或部分求出),则采用向后累积的形式,即先将 $Q$ 赋值为单位矩阵,依次计算 $Q\leftarrow Q_jQ, j=r,r-1,\cdots,1$ .这样可以减少运算所需的flop数.

Jacobi-Givens旋转

Jacobi-Givens旋转是指在特定二维坐标坐标平面内进行的旋转,其一般形式为 $\begin{equation*} G (i,k,\theta)= \begin{pmatrix} 1& \\ & \ddots& \\ & & c& \cdots& s& \\ & & \vdots& \ddots&\vdots& \\ & & -s&\cdots&c&\\ & & & & &\ddots&\\ & & & & & &1 \end{pmatrix}. \end{equation*}$ 其中 $c=\cos\theta$ , $s=\sin\theta$ , $c$ 和 $s$ 等分别位于第 $i$ 行和第 $k$ 行.

用 $G (i,j,\theta)^T$ 左乘即产生在 $(i,k)$ 坐标平面的 $\theta$ 角的旋转.运用Givens旋转可有选择地消去一些矩阵元素.

若 $x\in\mathbb{R}^n$ , $y=G (i,k,\theta)^T$ ,则 $\begin{equation*} y_j= \begin{cases} cx_i-sx_k& j=i,\\ sx_i+cx_k& j=k,\\ x_j& j\neq i,k. \end{cases} \end{equation*}$ 若令 $\begin{align*} c=\frac{x_i}{\sqrt{x_i^2+x_k^2}},\\ s=-\frac{x_k}{\sqrt{x_i^2+x_k^2}}, \end{align*}$ 则可将 $y_k$ 化为 $0$ .通过合理运算顺序安排,可以避免计算中溢出发生.本算法共需 $5$ 个flop和 $1$ 次平方根运算.

在应用Givens旋转时,同样有一些细节问题需要特别注意:

当用Givens矩阵进行左乘或右乘修正时,应注意到它只影响矩阵的两行或两列.
同Householder变换类似,多个Givens矩阵的乘积通常并不显式地出现.
可将每一个旋转变换对应于一个浮点数 $\rho$ ,从而将其存储在向量已消元的位置,通常在正弦较小时存储 $\mathrm{sgn}(c) s/2$ ,余弦较小时存储 $\mathrm{sign} (s) 2/c$ .可以由 $\rho$ 重新构造Givens矩阵.

除标准Givens变换外,还有所谓快速Givens变换,又称"免平方根 (square root-free) Givens变换".将在应用时介绍.

$QR$ 分解

一个 $m\times n$ 矩阵 $A$ 的 $QR$ 分解为 $A=QR$ ,其中 $Q\in \mathbb{R}^{m\times m}$ 是正交矩阵, $R\in \mathbb{R}^{n\times n}$ 是上三角矩阵.本节将讨论 $m\geqslant n$ 且 $A$ 为列满秩情形下的 $QR$ 分解,主要工具为Householder变换,Givens变换以及Gram-Schmidt正交化方法.

Householder $QR$ 分解

利用Householder变换进行 $QR$ 分解是直截了当的,可给出算法描述如下:

算法1(Householder QR分解) 给定 $A\in \mathbb{R}^{m\times m}$ , $m\geqslant n$ .在算法的第 $j$ 步,我们计算Householder矩阵 $H_j$ 满足:左乘 $A\leftarrow H_j^TA$ 使 $A(j+1:m,j)$ 部分化为 $0$ .如果 $Q=H_1\cdots H_n$ ,则 $Q^TA=R$ 是上三角阵.特别地, $A$ 的上三角部分被 $R$ 的上三角部分覆盖,第 $j$ 个Householder向量的基本部分存储于 $A(j+1:m,j)$ , $\forall j<m$ .

本算法共需要 $2n^2 (m-n/3)$ 个flop.

Givens $QR$ 方法

选择一个合适的消去的顺序,Givens $QR$ 分解也是简单有效的.如通过自左而右逐列考察,每一列中通过自下而上地进行Givens变换对 $A(j+1:m,j)$ 引入零元,共需 $3n^2 (m-n/3)$ 个flop即可将 $A$ 化为上三角形式.Givens变换矩阵累积起来即得到正交变换矩阵 $Q$ .我们也可以将 $(c,s)$ 用单个浮点数 $\rho$ 表示,从而将其存储在已化零的 $A(i,j)$ 中.

下面讨论快速Givens $QR$ 方法.

快速Givens方法思想就是当 $Q$ 是一系列Givens旋转的乘积时巧妙地表达它.具体地说, $Q$ 用一个矩阵对 $(M,D)$ 来表示,其中 $M^TM=D=\mathrm{diag} (d_i)$ 且 $d_i>0,\forall{i}$ .矩阵 $Q$ , $M$ 和 $D$ 通过公式 $Q=MD^{-1/2}$ 联系起来.易知 $Q$ 是正交的.若 $F$ 使 $F^TDF\equiv D_{\text{new}}$ 是对角阵,则取 $M_{\text{new}}=MF$ 满足 $M_{\text{new}}^TM_{\text{new}}=D_{\text{new}}$ .从而由快速Givens的表示形式 $(M,D)$ 做变换即可得到 $(M_{\text{new}},D_{\text{new}})$ .

以 $2$ 维的情形来说明以上思想.令 $x= (x_1,x_2)^T$ ,给定 $D=\mathrm{diag} (d_1,d_2),(d_1,d_2>0)$ .定义 $F_1=\left[\begin{smallmatrix}\beta_1& 1\\1& \alpha_1\end{smallmatrix}\right]$ ,则 $\begin{align*} F_1^Tx&=\begin{bmatrix} \beta_1x_1+x_2\\ x_1+\alpha_1x_2 \end{bmatrix},\\ F_1^TDF_1&= \begin{bmatrix} d_2+\beta_1^2d_1& d_1\beta_1+d_2\alpha_1\\ d_1\beta_1+d_2\alpha_1& d_1+\alpha_1^2d_2 \end{bmatrix} \equiv D_1. \end{align*}$ 若 $x_2\neq 0$ ,取 $\alpha_1=-x_1/x_2$ , $\beta_1=-\alpha_1\alpha_2/d_1$ ,则 $\begin{align*} F_1^Tx&= \begin{bmatrix} x_2 (1+\gamma_1)\\ 0 \end{bmatrix},\\ F_1^TDF_1&= \begin{bmatrix} d_2 (1+\gamma_1)\\ & d_1 (1+\gamma_1) \end{bmatrix}. \end{align*}$ 其中 $\gamma_1\equiv -\alpha_1\beta_1= (d_2/d_1) (x_1/x_2)^2$ .

类似地,如果 $x_1\neq 0$ ,则可定义 $F_2=\left[\begin{smallmatrix}1& \alpha_2\\ \beta_2& 1\end{smallmatrix}\right]$ ,其中 $\alpha_2\equiv -x_2/x_1,\beta_2 \equiv - (d_1/d_2)\alpha_2$ .则 $\begin{align*} F_2^Tx&= \begin{bmatrix} x_1 (1+\gamma_2)\\ 0 \end{bmatrix},\\ F_2^TDF_2&= \begin{bmatrix} d_1(1+\gamma_2)\\ & d_2 (1+\gamma_2) \end{bmatrix}. \end{align*}$ 其中 $\gamma_2=-\alpha_2\beta_2= (d_1/d_2) (x_2/x_1)^2$ .

这样,我们就完成了对该算法思想的描述.实际中选择变化类型使矩阵元素的"增长因子" $(1+\gamma_i)$ 以 $2$ 为界.

借助快速Givens变换,我们可以计算非奇异的 $M\in\mathbb{R}^{m\times m}$ 和正的 $d(1:m)$ 使得 $M^TA=T$ 为上三角阵,且 $M^TM=\mathrm{diag} (d_1,\cdots,d_m)$ .则 $A= (MD^{-1/2}) (D^{1/2}T)$ 是 $A$ 的 $QR$ 分解.这种算法需要 $2n^2 (m-n/3)$ 个flop,且不涉及平方根运算.

Givens方法对于稀疏矩阵问题,尤其是带状矩阵问题有特别的优势,例如对Hessenberg矩阵进行 $QR$ 分解只需要 $3n^2$ 个flop.

Gram-Schmidt方法

经典Gram-Schmidt方法(CGS)是我们熟知的,通过对 $A$ 的各列施行正交化手续,每步求出 $Q$ 和 $R$ 的一列.但CGS的数值特性非常差, $Q$ 的正交性经常会严重损失.通过改变计算次序,得到修正Gram-Schmidt方法(MGS),则可靠得多.

在MGS的第 $k$ 步,求出 $Q$ 的第 $k$ 列 $q_k$ 和 $R$ 的第 $k$ 行 $r_k^T$ .定义矩阵 $A^{(k)}\in\mathbb{R}^{m\times(n-k+1)}$ 为 $\begin{equation*} A-\sum\limits_{i=1}^{k-1}q_ir_i^T=\sum\limits_{i=k}^nq_ir_i^T\equiv (0, A^{(k)}). \end{equation*}$ 因此,如果 $A^{(k)}= (z,B)$ ,则有 $\begin{align*} r_{kk}&=\|z\|_2,\\ q_k&=z/r_{kk},\\ (r_{k,k+1},\cdots,r_{kn})&=q_k^TB. \end{align*}$ 随后计算外积 $\begin{equation*} A^{(k+1)}=B-q_k (r_{k,k+1},\cdots,r_{kn}) \end{equation*}$ 并进行下一步.这样我们就完成了对MGS算法第 $k$ 步的描述.这一算法共需 $2mn^2$ 个flop.但MGS的缺点在于 $A$ 病态时,不能保证良好的正交性.故仅当在被正交化的向量独立性很强时,才可用MGS求正交基.

满秩的最小二乘问题

要求找到向量 $x\in \mathbb{R}^n$ 使 $Ax=b$ ,其中 $A\in\mathbb{R}^{m\times n}$ 和 $b\in\mathbb{R}^m$ .若 $m\geqslant n$ ,即方程个数多于未知数个数,则称方程组 $Ax=b$ 是超定 (overdeterminant) 的,通常这样的方程组没有严格解.在实际中,我们考虑对适当选取的 $p$ ,极小化 $\|Ax-b\|_p$ .不同的范数给出不同的最优解,其中应用较多的是 $p=1,2,\infty$ 的情形.其中,2-范数情形便于解析求解.此时,该问题称为最小二乘(Least Square,LS)问题.

容易证明,若 $A$ 是列满秩的,则存在唯一的LS解 $x_{\text{LS}}$ ,它是对称正定线性方程组 $A^TAx_{\text{LS}}=A^Tb$ 的解,这一方程组称为法方程组.称 $r_{\text{LS}}=b-Ax_{\text{LS}}$ 为最小剩余,用 $L_{\text{LS}}\equiv \|b-Ax_{\text{LS}}\|_2$ 表示其大小.

当评估一个LS解 $\hat{x}_{\text{LS}}$ 的质量时,常考虑两方面的问题:

$\hat{x}_{\text{LS}}$ 有多靠近 $x_{\text{LS}}$ ?
与 $L_{\text{LS}}$ 相比, $\hat{L}_{\text{LS}}\equiv \|\hat{r}_{\text{LS}}\|_2=\|b-A\hat{x}_{\text{LS}}\|_2$ 有多小?

这两条标准在实际应用中常要考虑.

法方程组方法

算法2(法方程组方法) 给定 $A\in \mathbb{R}^{m\times n}$ 具有性质 $\mathrm{rank} (A)=n$ 和 $b\in \mathbb{R}^m$ .本算法计算LS问题 $\min\|Ax-b\|_2$ 的解 $x_{\text{LS}}$ .

计算 $C=A^TA$ 的下三角部分.
计算 $d=A^Tb$ .
计算Cholesky分解 $C=GG^T$ .
解 $Gy=d$ 和 $G^Tx_{\text{LS}}=y$ .

这一算法基于许多标准算法,共需要 $(m+n/3)n^2$ 个flop.经过分析,该算法计算解的精度依赖于 $A$ 的条件数的平方: $\begin{equation*} \frac{\|\hat{x}_{\text{LS}}-x_{\text{LS}}\|_2}{\|x_{\text{LS}}\|_2}\simeq \mu\kappa_2 (A)^2, \end{equation*}$ 其中 $\kappa_2 (A)$ 为 $A$ 的条件数.

用 $QR$ 分解求解LS问题

设 $A\in \mathbb{R}^{m\times n}$ , $m\geqslant n$ 且 $b\in \mathbb{R}^m$ 给定.若以计算得到 $A$ 的 $QR$ 分解 $\begin{equation*} Q^TA=R= \begin{bmatrix} R_1\\ O \end{bmatrix} \begin{array}[c]{l} \}n\\ \}m-n \end{array}, \end{equation*}$ 其中 $R_1$ 是 $n\times n$ 上三角阵.记 $\begin{equation*} \begin{bmatrix} c\\ d \end{bmatrix} \begin{array}[c]{l} \}n\\ \}m-n \end{array} \end{equation*}$ 则 $\begin{equation*} \|Ax-b\|_2=\|Q^TAx-Q^Tb\|_2=\|R_1x-c\|_2^2+\|d\|_2^2. \end{equation*}$ 于是, $\hat{x}_{\text{LS}}$ 由三角方程组 $R_1x_{\text{LS}}=c$ 定义,且 $L_{\text{LS}}=\|d\|_2$ .于是,只要求得 $A$ 的 $QR$ 分解,LS问题就很容易了.如采用Householder算法求得 $A$ 的 $QR$ 分解,则求解满秩LS问题,共需 $2n^2 (m-n/3)$ 个flop.

经分析,这种方法在 $\kappa_2 (A)\simeq 1/u$ 时失败.事实上,其解的灵敏性大致与 $\kappa_2 (A)+\rho_{\text{LS}}\kappa_2 (A)^2$ 成正比.因此, $QR$ 分解的方法适用范围比法方程组方法大一些.

基于MGS以及快速Givens方法的 $QR$ 分解也可用于求解LS问题,且解的稳定性相当.

病态矩阵的正交分解

如果 $A$ 是秩亏损的,则 $QR$ 分解不一定能给出 $\mathrm{span} (A)$ 的一组正交基.这问题可通过计算经过置换后的 $A$ 的 $QR$ 分解来解决,即 $A\Pi=QR$ ,其中 $\Pi$ 是置换阵.

如果右乘一个一般正交阵 $Z$ , $A$ 中的数据可被进一步压缩 $Q^TAZ=T$ . $Q$ 和 $Z$ 的选取以及相应的列主元的 $QR$ 分解正是本节要讨论的.

选主列的 $QR$ 分解

设 $A\in\mathbb{R}^{m\times n},(m\geqslant n)$ 且 $r\equiv\mathrm{rank} (A)<n$ ,则 $QR$ 分解不一定产生 $\mathrm{span}(A)$ 的一组正交基.但通过简单的修改,即计算分解 $\begin{equation*} Q^TA\Pi= \begin{bmatrix} R_{11}& R_{12}\\ O& O \end{bmatrix}. \end{equation*}$ 其中 $Q$ 为正交矩阵, $R_{11}$ 为 $r\times r$ 非奇异的上三角矩阵.于是 $Q$ 的前 $r$ 列给出 $\mathrm{span} (A)$ 的正交基.

$Q$ 和 $\Pi$ 分别是Householder矩阵和初等置换阵之积.假定对某个 $k$ ,我们已计算了Householder矩阵 $H_1,\cdots,H_{k-1}$ 和置换阵 $\Pi_1,\cdots,\Pi_{k-1}$ 使得 $\begin{equation*} (H_{k-1}\cdots H_1)A (\Pi_1\cdots\Pi_{k-1})=R^{(k-1)}= \begin{bmatrix} R_{11}^{(k-1)}& R_{12}^{(k-1)}\\ O& R_{22}^{(k-1)} \end{bmatrix}_{m\times n}, \end{equation*}$ 其中 $R_{11}^{(k-1)}$ 为 $(k-1)\times (k-1)$ 阶非奇异上三角阵.假定 $\begin{equation*} R_{22}^{(k-1)}= \begin{bmatrix} z_k^{(k-1)},\cdots,z_n^{(k-1)} \end{bmatrix}, \end{equation*}$ 且令 $z_p$ 为其中2-范数最大者.注意若 $k-1=\mathrm{rank} (A)$ ,则该最大模是 $0$ ,计算至此结束.否则,令 $\Pi_k$ 是互换第 $p$ 列和第 $k$ 列的置换阵,并确定Householder矩阵 $H_k$ 使 $R^{(k)}=H_kR^{(k-1)}\Pi_k$ 有 $R^{(k)} (k+1:m,k)=0$ .换句话说, $\Pi_k$ 把最大列移到前面, $H_k$ 把它的非对角元化为 $0$ .

利用对正交阵 $Q\in \mathbb{R}^{s\times s}$ $\begin{equation*} Q^Tz= \begin{bmatrix} \alpha\\ w \end{bmatrix} \Longrightarrow \|w\|_2^2=\|z\|_2^2-\alpha^2 \end{equation*}$ 的性质,我们只需修正旧的列范数得到新的列范数,即 $\begin{equation*} \|z^{(j)}\|_2^2=\|z^{(j-1)}\|_2^2-r_{kj}^2. \end{equation*}$ 从而使选主列的工作量从 $O(mn^2)$ 减少到 $O(mn)$ .这样,我们就完成了对该算法的描述.

本算法共需 $4mnr-2r^2 (m+n)+4r^3/3$ 个flop,其中 $r=\mathrm{rank} (A)$ .正交矩阵 $Q$ 可以以因子形式存储与 $A$ 的对角线之下.

完全正交分解

如果对上述算法产生的 $R$ ,从右边用一组适当的Householder矩阵相乘,则它可以进一步缩小.具体地说,我们可通过列满秩的 $QR$ 分解来计算 $\begin{equation*} Z_r\cdots Z_1 \begin{bmatrix} R_{11}^T\\ R_{12}^T \end{bmatrix} = \begin{bmatrix} T_{11}^T\\ O \end{bmatrix} \begin{array}[c]{l} \}r\\ \}n-r \end{array}. \end{equation*}$ 其中, $Z_i$ 为Householder变换, $T_{11}^T$ 是上三角矩阵,于是有 $\begin{equation*} Q^TAZ=T= \begin{bmatrix} T_{11}& O\\ O& O \end{bmatrix}, \end{equation*}$ 其中, $Z=\Pi Z_1\cdots Z_r$ .我们称其为完全正交分解.在这种情况下, $\mathrm{null} (A)=\mathrm{span} (Z (1:n,r+1:n))$ .

双对角化

矩阵的双对角化即通过正交变换将矩阵化为上双对角形 (上带宽为一的带状阵). $A$ 的双对角化与 $A^TA$ 的三对角化密切相关.通过对双对角阵 $B$ 的迭代运算,可以得到 $A$ 的奇异值分解(SVD).我们将其放在"对称特征值问题"一章中讨论.

假定 $A\in \mathbb{R}^{m\times n}$ 且 $m\geqslant n$ .我们要计算正交阵 $U_B\in \mathbb{R}^{m\times m}$ 和 $V_B\in\mathbb{R}^{n\times n}$ 使 $\begin{equation*} U_B^TAV_B= \begin{bmatrix} F\\ O \end{bmatrix}, \end{equation*}$ 其中, $F\in\mathbb{R}^{n\times n}$ 为上双对角阵,即上带宽为 $1$ .很容易通过直接的Householder变换,使 $\begin{align*} U_B&=U_1\cdots U_n,\\ V_B&=V1\cdots V_{n-2}. \end{align*}$ 其中,算法进行的第 $k$ 步, $U_k$ 将 $A (k+1:m,k)$ 化为 $0$ ,而 $V_k$ 将 $A (k,k+2:n)$ 化为 $0$ .本算法共需 $4mn^2-4n^3/3$ 个flop.

当 $m\gg n$ 时,若在进行双对角化之前先进行上三对角化,则可得到一个更为快速的双对角化算法,称为R双对角化.具体地说,假定我们计算一个正交阵 $Q\in \mathbb{R}^{m\times m}$ 使得 $\begin{equation*} Q^TA= \begin{bmatrix} R_1\\ O \end{bmatrix} \end{equation*}$ 是上三角阵.然后对 $R_1$ 进行双对角化 $\begin{equation*} U_R^TR_1V_B=B_1, \end{equation*}$ 则 $\begin{equation*} U=Q \begin{bmatrix} U_R\\ & I_{m-n} \end{bmatrix}, \end{equation*}$ 而 $\begin{equation*} U^TAV= \begin{bmatrix} B_1\\ O \end{bmatrix} \equiv B \end{equation*}$ 是 $A$ 的双对角化.它需要 $2mn^2+2n^3$ 个flop.可以看出, $m\geqslant \frac{5}{3}n$ 时,它的计算量要少一些.

秩亏损的LS问题

当 $A$ 是秩亏损的,则LS问题就有无穷多个解.可以证明,所有极小解的集合: $\chi \equiv\{x\in \mathbb{R}^n:\|Ax-b\|_2\text{极小}\}$ 是凸集,从而存在唯一元素具有极小2-范数,记为 $x_{\text{LS}}$ .

完全正交分解,SVD与 $x_{\text{LS}}$

任何完全正交分解都可用来计算 $x_{\text{LS}}$ .设 $\begin{equation*} Q^TAZ= \begin{bmatrix} T_{11}& O\\ O& O \end{bmatrix}\in \mathbb{R}^{n\times n} \end{equation*}$ 为 $A$ 的完全正交分解,其中 $T_{11}$ 为 $r\times r$ 上三角矩阵, $r$ 为 $A$ 的秩.则 $\begin{equation*} \|Ax-b\|_2^2=\| (Q^TAZ)Z^Tx-Q^Tb\|_2^2=\|T_{11}w-c\|_2^2+\|d\|_2^2, \end{equation*}$ 其中 $\begin{align*} Z^Tb&= \begin{bmatrix} w\\ y \end{bmatrix} \begin{array}[c]{l} \}r\\ \}n-r \end{array},\\ Q^Tb&= \begin{bmatrix} c\\ d \end{bmatrix} \begin{array}[c]{l} \}r\\ \}m-r \end{array}. \end{align*}$ 于是 $\begin{equation*} x_{\text{LS}}=Z \begin{bmatrix} T_{11}^{-1}c\\ O \end{bmatrix}. \end{equation*}$

SVD是非常"显露"的完全正交分解,当然可以用它来计算 $x_{\text{LS}}$ ,并有以下定理,它提供了 $x_{\text{LS}}$ 的一个简洁表达式和最小剩余量的范数 $L_{\text{LS}}$ .

定理1 假定 $U^TAV=\Sigma$ 是 $A\in R^{m\times n}$ 的SVD,且 $r=\mathrm{rank} (A)$ .如果 $U=[u_1,\cdots,u_m]$ , $V=[v_1,\cdots,v_n]$ 是按列划分的, $b\in \mathbb{R}^m$ .则 $\begin{equation*} x_{\text{LS}}=\sum\limits_{i=1}^r\frac{u_i^Tb}{\sigma_i}v_i \end{equation*}$ 使 $\|Ax-b\|_2$ 极小化,且是所有极小值点中2-范数最小的.而且 $\begin{equation*} L_{\text{LS}}^2\equiv \|Ax_{\text{LS}}-b\|_2^2=\sum\limits_{i=r+1}^m (u_i^Tb)^2. \end{equation*}$

LS问题的极小2-范数解 $x_{\text{LS}}$ 与 $A$ 的广义逆矩阵有如下关系:

设 $\begin{equation*} U^TAV=\Sigma=\mathrm{diag} (\sigma_1,\cdots,\sigma_r,0,\cdots,0) \end{equation*}$ 是 $A\in \mathbb{R}^{m\times n}$ 的SVD,则 $A$ 的广义逆(又称Moore-Penrose广义逆或加号广义逆) $A^+\in \mathbb{R}^{n\times m}$ 为 $\begin{equation*} A^+=V\Sigma^+U^T, \end{equation*}$ 其中 $\begin{equation*} \Sigma^+\equiv \mathrm{diag} (1/\sigma_1,\cdots,1/\sigma_r,0,\cdots,0)\in \mathbb{R}^{n\times m}. \end{equation*}$ 则 $\begin{equation*} x_{\text{LS}}=A^+b, \end{equation*}$ 且 $\begin{equation*} \rho_{\text{LS}}=\|(I-AA^+)b\|_2. \end{equation*}$

$A^+$ 是 $\min\limits_{X\in \mathbb{R}^{n\times m}}\|AX-I_m\|$ 的唯一的最小Frobenius范数解.

LS问题在 $A$ 秩亏损时变得相当困难, $x_{\text{LS}}$ 甚至不是 $A$ 中元素的的连续函数. $A$ 和 $b$ 的微小变化可能引起 $x_{\text{LS}}=A^+b$ 的任意大的变化.

选主列的 $QR$ 分解与基本解

假定 $A\in \mathbb{R}^{m\times n}$ 秩为 $r$ ,选主列的 $QR$ 分解给出 $\begin{equation*} A\Pi=QR, \end{equation*}$ 其中 $\begin{equation*} R= \begin{bmatrix} R_{11}& R_{12}\\ O& O \end{bmatrix}\in \mathbb{R}^{m\times n}, \end{equation*}$ $R_{11}$ 为 $r\times r$ 满秩上三角阵.则对任意 $x\in \mathbb{R}^n$ 有 $\begin{equation*} \|Ax-b\|_2^2=\| (Q^TA\Pi) (\Pi^Tx)- (Q^Tb)\|_2^2=\|R_{11}y- (c-R_{12}z)\|_2^2+\|d\|_2^2, \end{equation*}$ 其中 $\begin{align*} \Pi^Tx= \begin{bmatrix} y\\ z \end{bmatrix} \begin{array}[c]{l} \}r\\ \}n-r \end{array},\\ Q^Tb= \begin{bmatrix} c\\ d \end{bmatrix} \begin{array}[c]{l} \}r\\ \}m-r \end{array}. \end{align*}$ 于是,若 $x$ 是一个极小解,则有 $\begin{equation*} x=\Pi \begin{bmatrix} R_{11}^{-1} (c-R_{12}z)\\ z \end{bmatrix}. \end{equation*}$ 若取 $z=0$ ,则得到所谓基本解 $\begin{equation*} x_B=\Pi \begin{bmatrix} R_{11}^{-1}c\\ O \end{bmatrix}. \end{equation*}$ 可以看到, $x$ 至多只有 $r$ 个非零分量.但除非 $R_{12}=0$ ,基本解一般不是最小2-范数解.

欠定方程组的LS解

当 $m<n$ 时,我们称线性方程组 $Ax=b,A\in\mathbb{R}^{m\times n},b\in \mathbb{R}^m$ 是欠定的(underdetermined),它要么无解,要么有无穷多解.我们希望求出它的最小2-范数解.以下我们先讨论系数矩阵 $A$ 行满秩情形的算法,随后对秩亏损(但非矛盾)的情形加以考察.本节的讨论主要来自[1].

广义逆

设欠定方程组 $Ax=y$ ,其中 $A\in\mathbb{R}^{m\times n},(m<n)$ 为行满秩的系数矩阵.于是 $m\times m$ 矩阵 $AA^T$ 非奇异,且 $x=A^T (AA^T)^{-1}y$ 正是方程组的最小2-范数解,其中 $A^T (AA^T)^{-1}$ 是 $A$ 的Moore-Penrose广义逆 $A^+$ .

一个直接的算法是显式形成 $AA^T$ ,用Cholesky分解或 $QR$ 分解求解得到 $(AA^T)^{-1}y$ .但误差理论指出,若显式计算 $AA^T$ ,因条件数增大 ( $\kappa (AA^T)=\kappa (A)^2$ ) 可能会造成精度的严重丢失.因而这种方法并不可取.

$LU$ 分解

算法3 设通过全选主元的 $LU$ 分解得到 $P_1AP_2=LDU$ ,其中 $P_1$ , $P_2$ 为置换矩阵, $L$ 为 $m\times m$ 的单位下三角阵, $D$ 为 $m\times m$ 的的非奇异对角阵, $U$ 为 $m\times n$ 的单位上梯形阵.则 $Ax=y$ 的最小2-范数解可由以下步骤求出:

计算 $LU$ 分解 $P_1AP_2=LDU$ .
由向前消去法求下三角方程 $Lw=P_1y$ 的解.
计算 $DUz=w$ 的最小2-范数解,其详细做法将随后叙述.
取 $\tilde{x}=P_2\tilde{z}$ ,即要求的最小2-范数解.

在计算 $LU$ 分解时,只借助行选主元也能保证计算的稳定性,同时能够大大减少比较的次数.

$\tilde{z}$ 的计算可以通过显式形成 $UU^T$ 来进行,因 $U$ 是单位上梯形阵,其条件数会小于 $\kappa (A)$ .以下还会介绍其他的方法.

正交分解法

通过运用Householder变换或MGS对系数矩阵 $A$ 进行选主元的 $LQ$ 分解 (即 $QR$ 分解的列形式,其中 $Q$ 是正交列),得到 $\begin{equation*} PA=LQ, \end{equation*}$ 其中 $L$ 为 $m\times m$ 单位下三角阵, $Q$ 为 $m\times n$ 的矩阵满足 $QQ^T=D$ 为非奇异的对角阵.则最小2-范数解 $\begin{equation*} \tilde{z}=Q^TD^{-1}w, \end{equation*}$ 其中 $w$ 满足 $Lw=P_1y$ .

消去法

最小二乘解 $\tilde{z}$ 也可通过如下的正交分解过程得到.

算法4 借助右乘Householder变换或行MGS计算行正交的矩阵 $Q$ ,满足 $UQ=R$ 为 $m\times m$ 的单位下三角阵,且 $QQ^T=D_1$ 为 $m\times m$ 的非奇异对角阵.随后按如下步骤计算 $\tilde{z}$ :

用向后消去法求解 $Rz=D_1^{-1}w$ .
令 $\tilde{z}=Q^TD_1^{-1}z$ .

消去与投影法

我们考虑一个直接构造欠定方程组的最小2-范数解的方法.注意到方程组的所有解的 $w$ 都可表为 $\begin{equation*} w=u+v, \end{equation*}$ 其中 $u\in \mathrm{span} (A^T)$ , $v\in \mathrm{null} (A)$ 且 $\mathrm{span} (A^T)$ 与 $\mathrm{null} (A)$ 互为正交补.故 $\begin{equation*} \|w\|_2^2=\|u\|_2^2+\|v\|_2^2\geqslant \|u\|_2^2. \end{equation*}$ 从而 $\begin{equation*} \tilde{x}\equiv u \end{equation*}$ 是 $w$ 到 $\mathrm{span} (A^T)$ 的正交投影.于是我们采用下述步骤就可得到 $\tilde{x}$ :

求出方程组的任意非零解 $w$ .
求得 $\mathrm{null} (A)$ 的一组正交基 $\{s_1,\cdots,s_{n-m}\}$ .记 $S\equiv (s_1,\cdots,s_{n-m})$ .
做正交投影 $\begin{equation*} \tilde{x}=\left[I-S (S^T S)^{-1}S^T\right]w. \end{equation*}$

下面具体说明:

做 $LU$ 分解 $\begin{equation*} \hat{A}\equiv P_1AP_2=LDU, \end{equation*}$ 记 $U= [U_1,U_2]$ ,其中 $U_1$ 为 $m\times m$ 上三角阵.令 $\begin{equation*} z= \begin{bmatrix} U_1^{-1}D^{-1}L^{-1}P_1y\\ O \end{bmatrix}, \end{equation*}$ 则 $w\equiv P_2z$ 为方程组的一个解, $z$ 可通过两次消去法得到.
由于 $\mathrm{null} (A)=\mathrm{null} (UP_1^T)$ ,故向量 $s\in \mathrm{null} (A)$ 当且仅当 $P_2s\in \mathrm{null} (U)$ .于是我们只要计算 $\mathrm{null} (U')$ 的一组正交基.首先注意到通过对 $\begin{equation*} U_1s_1^{(1)}=-U_2s_1^{(2)} \end{equation*}$ 实施向后消去法,即可得到 $\begin{equation*} s_1= \begin{bmatrix} s_1^{(1)}\\ s_1^{(2)} \end{bmatrix} \in \mathrm{null} (U), \end{equation*}$ 其中 $s_1^{(2)}$ 为任意取定的 $n-m$ 维非零向量.得到 $s_1$ 后,通过(非选主元的)Gauss消去法将 $\begin{equation*} \begin{bmatrix} U\\ s_1^T \end{bmatrix} \end{equation*}$ 化为上梯形阵(实际只要对最后一行实施初等变换) $\begin{equation*} \begin{bmatrix} U\\ s_1'^T \end{bmatrix}. \end{equation*}$ 而后求得 $\begin{equation*} \begin{bmatrix} U\\ s_1'^T \end{bmatrix}s_2=0 \end{equation*}$ 的解 $s_2$ .类似可依次由 $s_1',\cdots,s_k'$ 得到 $s_{k+1}'$ .在此过程中一般地要进行后 $n-m$ 列的交换 $Q$ .记 $\begin{equation*} Q_1= \begin{bmatrix} I\\ & Q^T \end{bmatrix}, \end{equation*}$ 则 $\begin{gather*} \mathrm{null} (U')=\mathrm{span} (Q_1S),\\ \tilde{x}=P_2Q_1\left[I-S (S^T S)^{-1}S^T\right]z. \end{gather*}$
计算 $\begin{equation*} w=z-\sum\limits_{i=1}^{n-m}\frac{s_i^T z}{\|s_i\|^2}s_i. \end{equation*}$ 令 $\begin{equation*} \tilde{z}=P_2Q_1w. \end{equation*}$

秩亏损情形

在 $A$ 为行不满秩且方程组不矛盾时,我们称之为秩亏损情形.设 $A$ 为 $r$ 秩的 $m\times n$ 矩阵,且有分解 $\begin{equation*} \hat{A}\equiv P_1AP_2 =LDU. \end{equation*}$ 其中 $L$ 为 $m\times r$ 单位下梯形阵, $D$ 为 $r\times r$ 对角阵, $U$ 为 $r\times n$ 单位上梯形阵.记 $\begin{equation*} L= \begin{bmatrix} L_1\\ L_2 \end{bmatrix}, \end{equation*}$ 其中 $L_1$ 为 $r\times r$ 下三角阵,则 $Ax=y$ 可化为等价的 $\begin{equation*} UP_2^Tx=D^{-1} (L_1^{-1},O)P_1y, \end{equation*}$ 其中 $y ' = D^{-1} (L_1^{-1},O)P_1y$ 可通过向后消去法求解 $\begin{equation*} L_1 w= \hat{y} \end{equation*}$ 得到 $\begin{equation*} y'=D^{-1}w. \end{equation*}$ 其中 $\hat{y}$ 表示 $y$ 的前 $r$ 维分量.这样原方程组就化为 $\begin{equation*} Ux'=y', \end{equation*}$ 可通过之前的方法求解.

综合考虑算法的稳定性与复杂度,一般认为基于Householder变换的正交分解法较好,而投影方法在 $m$ 略小于 $n$ 时较为快速.

广义最小二乘问题

若已知 $A\in \mathbb{R}^{m\times n}$ , $B\in \mathbb{R}^{m\times m}$ , $b\in\mathbb{R}^m$ ,要求 $v_{\text{LS}}\in \mathbb{R}^m$ 使 $\begin{equation*} \|v_{\text{LS}}\|_2^2=\min\limits_{b=Ax+Bv}v^Tv. \end{equation*}$ 其中 $x\in \mathbb{R}^n$ ,这一类问题通常称为广义最小二乘问题.Paige (1979)[2]在假定 $A$ , $B$ 为满秩的情况下,提出了如下算法:

算法5 首先计算 $A$ 的 $QR$ 分解 $\begin{equation*} Q^TA= \begin{bmatrix} R_1\\ O \end{bmatrix}, \end{equation*}$ 其中, $Q=[Q_1,Q_2]$ .

然后确定正交阵 $Z\in \mathbb{R}^{m\times m}$ 使得 $\begin{equation*} Q_2^TBZ= \begin{bmatrix} O& S \end{bmatrix}, \end{equation*}$ 其中 $S\in \mathbb{R}^{n\times (m-n)}$ 是上三角阵.记 $Z=[Z_1\quad Z_2]$ .于是约束条件化为 $\begin{equation*} \begin{bmatrix} Q_1^T& b\\ Q_2^T& b \end{bmatrix} = \begin{bmatrix} R_1\\ O \end{bmatrix} + \begin{bmatrix} Q_1^TBZ_1& Q_1^TBZ_2\\ O& S \end{bmatrix} \begin{bmatrix} Z_1^T& O\\ Z_2^T& O \end{bmatrix}. \end{equation*}$ 方程下半部分可解出 $v$ : $\begin{gather*} Su=Q_2^T,\\ v=Z_2u. \end{gather*}$ 而其上半部确定了 $x$ : $\begin{equation*} R_1x=Q_1^Tb- (Q_1^TBZ_1Z_1^T+Q_1BZ_2Z_2^T)v=Q_1^Tb-Q_1^TBZ_2u. \end{equation*}$