5.11 加权Gram-Schmidt 分解

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

假设矩阵 A=QRA=QRA=QR ，其中 QmnQ_{mn}Qmn 为列满秩矩阵，RnnR_{nn}Rnn 为上三角方阵，带入上式得
x^=((QR)TW(QR))−1(QR)TWb=(RTQTWQR)−1RTQTWb=R−1(QTWQ)−1R−TRTQTWb=R−1(QTWQ)−1QTWb\mathbf{\hat{x}} = ((QR)^TW(QR))^{-1}(QR)^TW\mathbf{b}\\ = (R^TQ^TWQR)^{-1}R^TQ^TW\mathbf{b} \\ = R^{-1}(Q^TWQ)^{-1}R^{-T}R^TQ^TW\mathbf{b} \\ = R^{-1}(Q^TWQ)^{-1}Q^TW\mathbf{b} x^=((QR)TW(QR))−1(QR)TWb=(RTQTWQR)−1RTQTWb=R−1(QTWQ)−1R−TRTQTWb=R−1(QTWQ)−1QTWb

其中矩阵 QTWQQ^TWQQTWQ 求逆，为了简化求逆，我们希望其为对角阵，又当权重矩阵为单位阵时，QTEQ＝QTQQ^TEQ＝Q^TQQTEQ＝QTQ 必须为单位阵，以便和一般 Gram-Schmidt 分解保持一致。在这个约束下，各矩阵尺寸为 Qmn,WmmQ_{mn},W_{mm}Qmn,Wmm ，(QTWQ)nn(Q^TWQ)_{nn}(QTWQ)nn ，所以可以使 QTWQ=WnQ^TWQ=W_nQTWQ=Wn ，对角阵 WnW_nWn 为对角阵 WWW 前 nnn 个对角元素构成的对角阵，即 Wn=diag(w1,w2,⋯,wn)W_n = diag(w_1,w_2,\cdots,w_n)Wn=diag(w1,w2,⋯,wn) 。注意此时矩阵 QQQ 一般不再是正交阵，即 QTQ=EQ^TQ=EQTQ=E 不再成立。

我们先研究矩阵乘法 QTWQQ^TWQQTWQ ，首先 WQ=(Wq1,⋯,Wqn)WQ=(W\mathbf{q}_1,\cdots,W\mathbf{q}_n)WQ=(Wq1,⋯,Wqn) ，QTWQ=([q1T⋮,qnT])(Wq1,⋯,Wqn)=[q1TWq1,⋯,q1TWqn⋮,qnTWq1,⋯,qnTWqn]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)=⎣⎢⎡q1TWq1,⋯,q1TWqn⋮,qnTWq1,⋯,qnTWqn⎦⎥⎤ ，所以 QTWQQ^TWQQTWQ 第 iii 行第 jjj 列的元素为 qiTWqj\mathbf{q}^T_iW\mathbf{q}_jqiTWqj 。根据 QTWQ=WnQ^TWQ=W_nQTWQ=Wn 等式，则 qiTWqj=wifori=j,else0i≠j\mathbf{q}^T_iW\mathbf{q}_j = w_i \quad for \quad i=j, else \quad 0 \quad i \ne jqiTWqj=wifori=j,else0i=j ，从这个意义上说向量组 qi\mathbf{q}_iqi 关于矩阵 WWW 正交，因为不同向量的广义内积 qiTWqj\mathbf{q}^T_iW\mathbf{q}_jqiTWqj 为 000 ，相同向量广义内积为 wiw_iwi 。

根据 A=QRA=QRA=QR 和 QTWQ=WnQ^TWQ=W_nQTWQ=Wn 两个等式，可以推导出加权Gram-Schmidt 分解公式。

首先，a1=q1r11\mathbf{a}_1=\mathbf{q}_1 r_{11}a1=q1r11 两边左乘 q1TW\mathbf{q}^T_1Wq1TW 得，q1TWa1=q1TWq1r11=w1r11\mathbf{q}^T_1W\mathbf{a}_1=\mathbf{q}^T_1W\mathbf{q}_1 r_{11}=w_1r_{11}q1TWa1=q1TWq1r11=w1r11 ，又 q1TWa1=a1TWa1/r11\mathbf{q}^T_1W\mathbf{a}_1=\mathbf{a}^T_1W\mathbf{a}_1/r_{11}q1TWa1=a1TWa1/r11 ，所以 a1TWa1/r11=w1r11\mathbf{a}^T_1W\mathbf{a}_1/r_{11}=w_1r_{11}a1TWa1/r11=w1r11 得 r11=a1TWa1/w1r_{11}=\sqrt{\mathbf{a}^T_1W\mathbf{a}_1}/\sqrt{w_1}r11=a1TWa1/w1 ，令 a12=a1TWa1a^2_1=\mathbf{a}^T_1W\mathbf{a}_1a12=a1TWa1 为广义内积，则 r11=a1/w1r_{11}=a_1/\sqrt{w_1}r11=a1/w1 ，q1=a1/r11=w1a1/a1\mathbf{q}_1=\mathbf{a}_1/r_{11}=\sqrt{w_1}\mathbf{a}_1/a_1q1=a1/r11=w1a1/a1 。

其次，a2=q1r12+q2r22\mathbf{a}_2 = \mathbf{q}_1 r_{12} + \mathbf{q}_2 r_{22}a2=q1r12+q2r22 两边左乘 q1TW\mathbf{q}^T_1Wq1TW 得，q1TWa2=q1TWq1r12+q1TWq2r22\mathbf{q}^T_1W\mathbf{a}_2 = \mathbf{q}^T_1W\mathbf{q}_1 r_{12} + \mathbf{q}^T_1W\mathbf{q}_2 r_{22}q1TWa2=q1TWq1r12+q1TWq2r22 ，根据 qiTWqj=wifori=j,else0i≠j\mathbf{q}^T_iW\mathbf{q}_j = w_i \quad for \quad i=j, else \quad 0 \quad i \ne jqiTWqj=wifori=j,else0i=j ，得 q1TWa2=q1TWq1r12+0=w1r12\mathbf{q}^T_1W\mathbf{a}_2 = \mathbf{q}^T_1W\mathbf{q}_1 r_{12} + 0 = w_1 r_{12}q1TWa2=q1TWq1r12+0=w1r12 得 r12=q1TWa2/w1r_{12} = \mathbf{q}^T_1W\mathbf{a}_2/w_1r12=q1TWa2/w1 为广义投影坐标值。

a2=q1r12+q2r22\mathbf{a}_2 = \mathbf{q}_1 r_{12} + \mathbf{q}_2 r_{22}a2=q1r12+q2r22 两边左乘 q2TW\mathbf{q}^T_2Wq2TW 得，q2TWa2=q2TWq1r12+q2TWq2r22\mathbf{q}^T_2W\mathbf{a}_2 = \mathbf{q}^T_2W\mathbf{q}_1 r_{12} + \mathbf{q}^T_2W\mathbf{q}_2 r_{22}q2TWa2=q2TWq1r12+q2TWq2r22 ，得 q2TWa2=0+w2r22\mathbf{q}^T_2W\mathbf{a}_2 = 0 + w_2 r_{22}q2TWa2=0+w2r22 ，得 w2r22=q2TW(a2−q1r12)w_2 r_{22} = \mathbf{q}^T_2W(\mathbf{a}_2- \mathbf{q}_1 r_{12})w2r22=q2TW(a2−q1r12) ，又 q2=(a2−q1r12)/r22\mathbf{q}_2 = (\mathbf{a}_2 - \mathbf{q}_1 r_{12})/r_{22}q2=(a2−q1r12)/r22 带入得 w2r222=(a2−q1r12)W(a2−q1r12)w_2 r^2_{22} = (\mathbf{a}_2 - \mathbf{q}_1 r_{12})W(\mathbf{a}_2 - \mathbf{q}_1 r_{12})w2r222=(a2−q1r12)W(a2−q1r12) ，令 a22=(a2−q1r12)W(a2−q1r12)a^2_2=(\mathbf{a}_2 - \mathbf{q}_1 r_{12})W(\mathbf{a}_2 - \mathbf{q}_1 r_{12})a22=(a2−q1r12)W(a2−q1r12) 为广义内积，所以 r22=a2/w2r_{22} = a_2/\sqrt{w_2}r22=a2/w2 。所以q2=w2(a2−q1r12)/a2\mathbf{q}_2 = \sqrt{w_2}(\mathbf{a}_2 - \mathbf{q}_1 r_{12})/a_2q2=w2(a2−q1r12)/a2 。

最后同理对任意列向量 ai=q1r1i+q2r2i+⋯+qirii\mathbf{a}_i = \mathbf{q}_1 r_{1i} + \mathbf{q}_2 r_{2i} + \cdots + \mathbf{q}_i r_{ii}ai=q1r1i+q2r2i+⋯+qirii ，两边左乘 qjTW,j<i\mathbf{q}^T_jW,j < iqjTW,j<i 得
qjTWai=qjTWq1r1i+qjTWq2r2i+⋯+qjTWqirii=qjTWqjrji=wjrji\mathbf{q}^T_jW\mathbf{a}_i = \mathbf{q}^T_jW\mathbf{q}_1 r_{1i} + \mathbf{q}^T_jW\mathbf{q}_2 r_{2i} + \cdots + \mathbf{q}^T_jW\mathbf{q}_i r_{ii} = \mathbf{q}^T_jW\mathbf{q}_j r_{ji} = w_j r_{ji}qjTWai=qjTWq1r1i+qjTWq2r2i+⋯+qjTWqirii=qjTWqjrji=wjrji 所以 rji=qjTWai/wjr_{ji} = \mathbf{q}^T_jW\mathbf{a}_i/w_jrji=qjTWai/wj 。两边左乘 qiTW\mathbf{q}^T_iWqiTW 得
qiTWai=qiTWq1r1i+qiTWq2r2i+⋯+qiTWqirii=qiTWqirii=wirii\mathbf{q}^T_iW\mathbf{a}_i = \mathbf{q}^T_iW\mathbf{q}_1 r_{1i} + \mathbf{q}^T_iW\mathbf{q}_2 r_{2i} + \cdots + \mathbf{q}^T_iW\mathbf{q}_i r_{ii} = \mathbf{q}^T_iW\mathbf{q}_i r_{ii} = w_i r_{ii}qiTWai=qiTWq1r1i+qiTWq2r2i+⋯+qiTWqirii=qiTWqirii=wirii 所以 rii=qiTWai/wir_{ii} = \mathbf{q}^T_iW\mathbf{a}_i/w_irii=qiTWai/wi ，根据 ai=q1r1i+q2r2i+⋯+qirii\mathbf{a}_i = \mathbf{q}_1 r_{1i} + \mathbf{q}_2 r_{2i} + \cdots + \mathbf{q}_i r_{ii}ai=q1r1i+q2r2i+⋯+qirii 最后化简得，令 ai2=(ai−(q1r1i+q2r2i+⋯))W(ai−(q1r1i+q2r2i+⋯))a^2_i=(\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))ai2=(ai−(q1r1i+q2r2i+⋯))W(ai−(q1r1i+q2r2i+⋯)) 为广义内积，所以 rii=ai/wir_{ii} = a_i/\sqrt{w_i}rii=ai/wi 。所以 qi=wi(ai−(q1r1i+q2r2i+⋯))/ai\mathbf{q}_i = \sqrt{w_i}(\mathbf{a}_i - (\mathbf{q}_1 r_{1i} + \mathbf{q}_2 r_{2i} + \cdots))/a_iqi=wi(ai−(q1r1i+q2r2i+⋯))/ai 。

上面介绍的方法就是经典的加权Gram-Schmidt分解，过程和经典的Gram-Schmidt分解一样，差别在于上三角矩阵 RRR 的元素 rji=qjTWai/wjr_{ji} = \mathbf{q}^T_jW\mathbf{a}_i/w_jrji=qjTWai/wj 为广义投影 qjTWai\mathbf{q}^T_jW\mathbf{a}_iqjTWai 除以权重 wjw_jwj ，rii=ai/wir_{ii} = a_i/\sqrt{w_i}rii=ai/wi 为广义内积 aia_iai 除以权重 wi\sqrt{w_i}wi ，当广义内积趋近 000 时，则方程是病态。qi=wi(ai−(q1r1i+q2r2i+⋯))/ai\mathbf{q}_i = \sqrt{w_i}(\mathbf{a}_i - (\mathbf{q}_1 r_{1i} + \mathbf{q}_2 r_{2i} + \cdots))/a_iqi=wi(ai−(q1r1i+q2r2i+⋯))/ai 是垂直分量的标准化 (ai−(q1r1i+q2r2i+⋯))/ai(\mathbf{a}_i - (\mathbf{q}_1 r_{1i} + \mathbf{q}_2 r_{2i} + \cdots))/a_i(ai−(q1r1i+q2r2i+⋯))/ai 乘以权重 wi\sqrt{w_i}wi 。

改进加权 Gram-Schmidt 分解

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

1、r11=a1/w1r_{11}=a_1/\sqrt{w_1}r11=a1/w1 ，q1=w1a1/a1\mathbf{q}_1=\sqrt{w_1}\mathbf{a}_1/a_1q1=w1a1/a1 ，r1i=q1TWai/w1r_{1i} = \mathbf{q}^T_1W\mathbf{a}_i/w_1r1i=q1TWai/w1 ，ai=ai−q1r1i,i>1\mathbf{a}_i = \mathbf{a}_i - \mathbf{q}_1 r_{1i}, i>1ai=ai−q1r1i,i>1 。令 b0=b\mathbf{b}_0=\mathbf{b}b0=b ，注意对向量 b\mathbf{b}b 也要同样处理，即令 δ1=q1TWb/w1\delta_1 = \mathbf{q}^T_1W\mathbf{b}/w_1δ1=q1TWb/w1 ，b=b−q1δ1\mathbf{b} = \mathbf{b} - \mathbf{q}_1\delta_1b=b−q1δ1 。

2、r22=a2/w2r_{22}=a_2/\sqrt{w_2}r22=a2/w2 ，q2=w2a2/a2\mathbf{q}_2=\sqrt{w_2}\mathbf{a}_2/a_2q2=w2a2/a2 ，r2i=q2TWai/w2r_{2i} = \mathbf{q}^T_2W\mathbf{a}_i/w_2r2i=q2TWai/w2 ，ai=ai−q2r2i,i>2\mathbf{a}_i = \mathbf{a}_i - \mathbf{q}_2 r_{2i}, i>2ai=ai−q2r2i,i>2 。注意对向量 b\mathbf{b}b 也要同样处理，即令 δ2=q2TWb/w2\delta_2 = \mathbf{q}^T_2W\mathbf{b}/w_2δ2=q2TWb/w2 ，b=b−q2δ2\mathbf{b} = \mathbf{b} - \mathbf{q}_2\delta_2b=b−q2δ2 。

3、rii=ai/wir_{ii}=a_i/\sqrt{w_i}rii=ai/wi ，qi=wiai/ai\mathbf{q}_i=\sqrt{w_i}\mathbf{a}_i/a_iqi=wiai/ai ，rij=qiTWaj/wir_{ij} = \mathbf{q}^T_iW\mathbf{a}_j/w_irij=qiTWaj/wi ，aj=aj−qirij,j>i\mathbf{a}_j = \mathbf{a}_j - \mathbf{q}_i r_{ij}, j>iaj=aj−qirij,j>i 。注意对向量 b\mathbf{b}b 也要同样处理，即令 δi=qiTWb/wi\delta_i = \mathbf{q}^T_iW\mathbf{b}/w_iδi=qiTWb/wi ，b=b−qiδi\mathbf{b} = \mathbf{b} - \mathbf{q}_i\delta_ib=b−qiδi 。

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

令向量 d=(δ1,⋯,δm)\mathbf{d} = (\delta_1,\cdots,\delta_m)d=(δ1,⋯,δm) ，因为 d=(QTWQ)−1QTWb0\mathbf{d} = (Q^TWQ)^{-1}Q^TW\mathbf{b}_0d=(QTWQ)−1QTWb0 ，则最优近似解为 x^=R−1d\mathbf{\hat{x}} = R^{-1}\mathbf{d}x^=R−1d ，即 x^i=(δi−∑j=i+1n(δjx^j))/rii\hat{x}_i = (\delta_i - \sum^{n}_{j=i+1} (\delta_j\hat{x}_j))/r_{ii}x^i=(δi−∑j=i+1n(δjx^j))/rii。

阻尼倒数法

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

权重排序

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