5.11 加权Gram-Schmidt 分解_施普林格改进的schmidt分割算法

5.11 加权Gram-Schmidt 分解

前面介绍了加权最小二乘法的解为 $\hat{\mathbf{x}}=(A^{T}DA{)}^{-1}{A}^{T}D\mathbf{b}$ ，其中 $D$ 为权重矩阵，最常用的是对角阵，每个对角元素表示对应测量点的权重，均为正值，元素值越大表示该测量点越重要。本节仅研究权重矩阵为对角阵的情况，权重矩阵记为 $W=diag(w_1,w_2,\cdots ,w_m)$ ，$w_i$ 为测量点 $i$ 的权重。

假设矩阵 $A=QR$ ，其中 $Q_{mn}$ 为列满秩矩阵，$R_{nn}$ 为上三角方阵，带入上式得
$$\begin{gathered} \hat{\mathbf{x}}=((QR{)}^{T}W(QR){)}^{-1}(QR{)}^{T}W\mathbf{b} \\ =(R^T{Q}^{T}WQR{)}^{-1}{R}^T{Q}^{T}W\mathbf{b} \\ =R^{-1}(Q^{T}WQ{)}^{-1}{R}^{-T}{R}^T{Q}^{T}W\mathbf{b} \\ =R^{-1}(Q^{T}WQ{)}^{-1}{Q}^{T}W\mathbf{b} \end{gathered}$$

其中矩阵 $Q^{T}WQ$ 求逆，为了简化求逆，我们希望其为对角阵，又当权重矩阵为单位阵时，$Q^{T}EQ＝Q^{T}Q$ 必须为单位阵，以便和一般 Gram-Schmidt 分解保持一致。在这个约束下，各矩阵尺寸为 $Q_{mn},W_{mm}$ ，$(Q^{T}WQ{)}_{nn}$ ，所以可以使 $Q^{T}WQ=W_n$ ，对角阵 $W_n$ 为对角阵 $W$ 前 $n$ 个对角元素构成的对角阵，即 $W_n=diag(w_1,w_2,\cdots ,w_n)$ 。注意此时矩阵 $Q$ 一般不再是正交阵，即 $Q^{T}Q=E$ 不再成立。

我们先研究矩阵乘法 $Q^{T}WQ$ ，首先 $WQ=(W{\mathbf{q}}_1,\cdots ,W{\mathbf{q}}_n)$ ， Q T W Q = ( [ q 1 T ⋮ , q n T ] ) ( W q 1 , ⋯   , W q n ) = [ q 1 T W q 1 , ⋯   , q 1 T W q n ⋮ , q n T W q 1 , ⋯   , q n T W q n ] Q^TWQ=(\left[

$$\begin{matrix} \mathbf{q}^T_1 \\ \vdots,\\ \mathbf{q}^T_n \end{matrix}$$ \right])(W\mathbf{q}_1,\cdots,W\mathbf{q}_n)=\left[ $$\begin{matrix} \mathbf{q}^T_1W\mathbf{q}_1,\cdots, \mathbf{q}^T_1W\mathbf{q}_n\\ \vdots,\\ \mathbf{q}^T_nW\mathbf{q}_1,\cdots, \mathbf{q}^T_nW\mathbf{q}_n \end{matrix}$$ \right] QTWQ=(⎣⎢⎡​q1T​⋮,qnT​​⎦⎥⎤​)(Wq1​,⋯,Wqn​)=⎣⎢⎡​q1T​Wq1​,⋯,q1T​Wqn​⋮,qnT​Wq1​,⋯,qnT​Wqn​​⎦⎥⎤​ ，所以 $Q^{T}WQ$ 第 $i$ 行第 $j$ 列的元素为 ${\mathbf{q}}_i^{T}W{\mathbf{q}}_j$ 。根据 $Q^{T}WQ=W_n$ 等式，则 ${\mathbf{q}}_i^{T}W{\mathbf{q}}_j=w_{i}fori=j,else0i\ne j$ ，从这个意义上说向量组 ${\mathbf{q}}_i$ 关于矩阵 $W$ 正交，因为不同向量的广义内积 ${\mathbf{q}}_i^{T}W{\mathbf{q}}_j$ 为 $0$ ，相同向量广义内积为 $w_i$ 。

根据 $A=QR$ 和 $Q^{T}WQ=W_n$ 两个等式，可以推导出加权Gram-Schmidt 分解公式。

首先，${\mathbf{a}}_1={\mathbf{q}}_1{r}_{11}$ 两边左乘 ${\mathbf{q}}_1^{T}W$ 得，${\mathbf{q}}_1^{T}W{\mathbf{a}}_1={\mathbf{q}}_1^{T}W{\mathbf{q}}_1{r}_{11}=w_1{r}_{11}$ ，又 ${\mathbf{q}}_1^{T}W{\mathbf{a}}_1={\mathbf{a}}_1^{T}W{\mathbf{a}}_1/r_{11}$ ，所以 ${\mathbf{a}}_1^{T}W{\mathbf{a}}_1/r_{11}=w_1{r}_{11}$ 得 $r_{11}=\sqrt{{\mathbf{a}}_1^{T}W{\mathbf{a}}_1}/\sqrt{w_1}$ ，令 $a_1^2={\mathbf{a}}_1^{T}W{\mathbf{a}}_1$ 为广义内积，则 $r_{11}=a_1/\sqrt{w_1}$ ，${\mathbf{q}}_1={\mathbf{a}}_1/r_{11}=\sqrt{w_1}{\mathbf{a}}_1/a_1$ 。

其次，${\mathbf{a}}_2={\mathbf{q}}_1{r}_{12}+{\mathbf{q}}_2{r}_{22}$ 两边左乘 ${\mathbf{q}}_1^{T}W$ 得，${\mathbf{q}}_1^{T}W{\mathbf{a}}_2={\mathbf{q}}_1^{T}W{\mathbf{q}}_1{r}_{12}+{\mathbf{q}}_1^{T}W{\mathbf{q}}_2{r}_{22}$ ，根据 ${\mathbf{q}}_i^{T}W{\mathbf{q}}_j=w_{i}fori=j,else0i\ne j$ ，得 ${\mathbf{q}}_1^{T}W{\mathbf{a}}_2={\mathbf{q}}_1^{T}W{\mathbf{q}}_1{r}_{12}+0=w_1{r}_{12}$ 得 $r_{12}={\mathbf{q}}_1^{T}W{\mathbf{a}}_2/w_1$ 为广义投影坐标值。

${\mathbf{a}}_2={\mathbf{q}}_1{r}_{12}+{\mathbf{q}}_2{r}_{22}$ 两边左乘 ${\mathbf{q}}_2^{T}W$ 得，${\mathbf{q}}_2^{T}W{\mathbf{a}}_2={\mathbf{q}}_2^{T}W{\mathbf{q}}_1{r}_{12}+{\mathbf{q}}_2^{T}W{\mathbf{q}}_2{r}_{22}$ ，得 ${\mathbf{q}}_2^{T}W{\mathbf{a}}_2=0+w_2{r}_{22}$ ，得 $w_2{r}_{22}={\mathbf{q}}_2^{T}W({\mathbf{a}}_2-{\mathbf{q}}_1{r}_{12})$ ，又 ${\mathbf{q}}_2=({\mathbf{a}}_2-{\mathbf{q}}_1{r}_{12})/r_{22}$ 带入得 $w_2{r}_{22}^2=({\mathbf{a}}_2-{\mathbf{q}}_1{r}_{12})W({\mathbf{a}}_2-{\mathbf{q}}_1{r}_{12})$ ，令 $a_2^2=({\mathbf{a}}_2-{\mathbf{q}}_1{r}_{12})W({\mathbf{a}}_2-{\mathbf{q}}_1{r}_{12})$ 为广义内积，所以 $r_{22}=a_2/\sqrt{w_2}$ 。所以${\mathbf{q}}_2=\sqrt{w_2}({\mathbf{a}}_2-{\mathbf{q}}_1{r}_{12})/a_2$ 。

最后同理对任意列向量 ${\mathbf{a}}_i={\mathbf{q}}_1{r}_{1i}+{\mathbf{q}}_2{r}_{2i}+\cdots +{\mathbf{q}}_i{r}_{ii}$ ，两边左乘 ${\mathbf{q}}_j^{T}W,j<i$ 得
${\mathbf{q}}_j^{T}W{\mathbf{a}}_i={\mathbf{q}}_j^{T}W{\mathbf{q}}_1{r}_{1i}+{\mathbf{q}}_j^{T}W{\mathbf{q}}_2{r}_{2i}+\cdots +{\mathbf{q}}_j^{T}W{\mathbf{q}}_i{r}_{ii}={\mathbf{q}}_j^{T}W{\mathbf{q}}_j{r}_{ji}=w_j{r}_{ji}$ 所以 $r_{ji}={\mathbf{q}}_j^{T}W{\mathbf{a}}_i/w_j$ 。两边左乘 ${\mathbf{q}}_i^{T}W$ 得
${\mathbf{q}}_i^{T}W{\mathbf{a}}_i={\mathbf{q}}_i^{T}W{\mathbf{q}}_1{r}_{1i}+{\mathbf{q}}_i^{T}W{\mathbf{q}}_2{r}_{2i}+\cdots +{\mathbf{q}}_i^{T}W{\mathbf{q}}_i{r}_{ii}={\mathbf{q}}_i^{T}W{\mathbf{q}}_i{r}_{ii}=w_i{r}_{ii}$ 所以 $r_{ii}={\mathbf{q}}_i^{T}W{\mathbf{a}}_i/w_i$ ，根据 ${\mathbf{a}}_i={\mathbf{q}}_1{r}_{1i}+{\mathbf{q}}_2{r}_{2i}+\cdots +{\mathbf{q}}_i{r}_{ii}$ 最后化简得，令 $a_i^2=({\mathbf{a}}_i-({\mathbf{q}}_1{r}_{1i}+{\mathbf{q}}_2{r}_{2i}+\cdots ))W({\mathbf{a}}_i-({\mathbf{q}}_1{r}_{1i}+{\mathbf{q}}_2{r}_{2i}+\cdots ))$ 为广义内积，所以 $r_{ii}=a_i/\sqrt{w_i}$ 。所以 ${\mathbf{q}}_i=\sqrt{w_i}({\mathbf{a}}_i-({\mathbf{q}}_1{r}_{1i}+{\mathbf{q}}_2{r}_{2i}+\cdots ))/a_i$ 。

上面介绍的方法就是经典的加权Gram-Schmidt分解，过程和经典的Gram-Schmidt分解一样，差别在于上三角矩阵 $R$ 的元素 $r_{ji}={\mathbf{q}}_j^{T}W{\mathbf{a}}_i/w_j$ 为广义投影 ${\mathbf{q}}_j^{T}W{\mathbf{a}}_i$ 除以权重 $w_j$ ，$r_{ii}=a_i/\sqrt{w_i}$ 为广义内积 $a_i$ 除以权重 $\sqrt{w_i}$ ，当广义内积趋近 $0$ 时，则方程是病态。${\mathbf{q}}_i=\sqrt{w_i}({\mathbf{a}}_i-({\mathbf{q}}_1{r}_{1i}+{\mathbf{q}}_2{r}_{2i}+\cdots ))/a_i$ 是垂直分量的标准化 $({\mathbf{a}}_i-({\mathbf{q}}_1{r}_{1i}+{\mathbf{q}}_2{r}_{2i}+\cdots ))/a_i$ 乘以权重 $\sqrt{w_i}$ 。

改进加权 Gram-Schmidt 分解

如同改进Gram-Schmidt分解，可以按行计算矩阵 $R$ 和先减去投影分量，得到改进加权 Gram-Schmidt 分解：令 $a_i^2={\mathbf{a}}_i^{T}W{\mathbf{a}}_i$ 为广义内积。

1、$r_{11}=a_1/\sqrt{w_1}$ ，${\mathbf{q}}_1=\sqrt{w_1}{\mathbf{a}}_1/a_1$ ，$r_{1i}={\mathbf{q}}_1^{T}W{\mathbf{a}}_i/w_1$ ，${\mathbf{a}}_i={\mathbf{a}}_i-{\mathbf{q}}_1{r}_{1i},i>1$ 。令 ${\mathbf{b}}_0=\mathbf{b}$ ，注意对向量 $\mathbf{b}$ 也要同样处理，即令 ${\delta }_1={\mathbf{q}}_1^{T}W\mathbf{b}/w_1$ ，$\mathbf{b}=\mathbf{b}-{\mathbf{q}}_1{\delta }_1$ 。

2、$r_{22}=a_2/\sqrt{w_2}$ ，${\mathbf{q}}_2=\sqrt{w_2}{\mathbf{a}}_2/a_2$ ，$r_{2i}={\mathbf{q}}_2^{T}W{\mathbf{a}}_i/w_2$ ，${\mathbf{a}}_i={\mathbf{a}}_i-{\mathbf{q}}_2{r}_{2i},i>2$ 。注意对向量 $\mathbf{b}$ 也要同样处理，即令 ${\delta }_2={\mathbf{q}}_2^{T}W\mathbf{b}/w_2$ ，$\mathbf{b}=\mathbf{b}-{\mathbf{q}}_2{\delta }_2$ 。

3、$r_{ii}=a_i/\sqrt{w_i}$ ，${\mathbf{q}}_i=\sqrt{w_i}{\mathbf{a}}_i/a_i$ ，$r_{ij}={\mathbf{q}}_i^{T}W{\mathbf{a}}_j/w_i$ ，${\mathbf{a}}_j={\mathbf{a}}_j-{\mathbf{q}}_i{r}_{ij},j>i$ 。注意对向量 $\mathbf{b}$ 也要同样处理，即令 ${\delta }_i={\mathbf{q}}_i^{T}W\mathbf{b}/w_i$ ，$\mathbf{b}=\mathbf{b}-{\mathbf{q}}_i{\delta }_i$ 。

一直计算，直到 $i=n$ 结束。

令向量 $\mathbf{d}=({\delta }_1,\cdots ,{\delta }_m)$ ，因为 $\mathbf{d}=(Q^{T}WQ{)}^{-1}{Q}^{T}W{\mathbf{b}}_0$ ，则最优近似解为 $\hat{\mathbf{x}}=R^{-1}\mathbf{d}$ ，即 ${\hat{x}}_i=({\delta }_i-{\sum }_{j=i+1}^n({\delta }_j{\hat{x}}_j))/r_{ii}$。

阻尼倒数法

计算 ${\mathbf{q}}_i=\sqrt{w_i}{\mathbf{a}}_i/a_i$ 涉及到除以广义内积 $a_i$ ，为了增加数值稳定性，可以采用阻尼倒数法，具体方法和QR分解的阻尼倒数法一样，这里省略。

权重排序

计算上三角阵 $R$ 元素 $r_{ij}={\mathbf{q}}_i^{T}W{\mathbf{a}}_j/w_i$ 涉及除以权重 $w_i$ ，如果其为或者趋近 $0$ ，则也会导致数值不稳定。为了解决这个问题，可采用权重排序计算。即按照权重从大到小排序 $w_i>w_{i+1}$ ，根据权重排序结果调整方程 $A\mathbf{x}=\mathbf{b}$ 中子方程的顺序，使之与权重顺序对应。这样只要前 $n$ 个权重不趋近 $0$ 即可保证数值稳定。如果所有权重大小差不多，则不必排序，可直接计算。注意此时不能采用阻尼倒数法。