30 正交三角分解
考虑如下线性模型的估计问题: \begin{align} \boldsymbol y = X \boldsymbol\beta + \boldsymbol\varepsilon, \tag{30.1} \end{align} 其中X是n \times p已知矩阵(n>p), \boldsymbol y是n维因变量观测值向量, \boldsymbol\beta是p维未知模型系数向量, \boldsymbol\varepsilon是n维未知的模型随机误差向量。 用最小二乘法估计\boldsymbol\beta, 即求\boldsymbol\beta使得 \begin{align} g(\boldsymbol\beta) = \| \boldsymbol y - X \boldsymbol\beta \|^2 \tag{30.2} \end{align} 最小(本节中的向量范数\| \boldsymbol x \|都表示欧式空间中的长度\sqrt{\boldsymbol x^T \boldsymbol x}), 归结为求解如下正规方程 \begin{align} X^T X \boldsymbol\beta = X^T \boldsymbol y. \tag{30.3} \end{align} 此方程一定有解, 当X列满秩时有唯一解\hat{\boldsymbol\beta} = (X^T X)^{-1} X^T \boldsymbol y。
考虑正规方程(30.3)的数值计算问题。 我们需要解出\boldsymbol\beta(记作\hat{\boldsymbol\beta}), 并计算残差平方和 \begin{align} \text{SSE} =& \| \boldsymbol y - X \hat{\boldsymbol\beta} \|^2. \end{align} 在X列满秩时, 正规方程(30.3)的系数矩阵X^T X是正定阵, 可以用Cholesky分解方法高效地计算出\hat{\boldsymbol\beta}和SSE, 并且所需的存储空间也只有O(p^2)阶。 但是,X^T X相当于做了平方, 当X^T X条件数\kappa_2(X^T X)很大时, 直接求解正规方程(30.3)会引入很大误差。 有结果表明(参见 Stewart (1973) 定理5.2.4), \hat{\boldsymbol\beta}的相对误差有如下控制式: \begin{align} \frac{ \| \Delta \hat{\boldsymbol\beta} \|}{ \| \hat{\boldsymbol\beta} \|} \leq& \sqrt{\kappa_2(X^T X)} \frac{ \| \Delta \hat{\boldsymbol y} \|}{\| \hat{\boldsymbol y}\| }, \end{align} 其中\hat{\boldsymbol y} = X \hat{\boldsymbol\beta}, \Delta \hat{\boldsymbol y}是观测值\boldsymbol y中的误差引起的\hat{\boldsymbol y}的误差, \Delta \hat{\boldsymbol\beta}是观测值\boldsymbol y中的误差引起的\hat{\boldsymbol\beta}的误差。 所以, 应该寻找计算更稳定的算法来求解最小二乘问题(30.2)。 另外,在样本量n很大时, 计算X^T X时的累积误差可能导致其实际不正定, 使得Cholesky分解算法失败。 新的算法最好能够在X非列满秩时也可以求出适当的最小二乘解。 直接利用X计算而不是对X^T X计算可以提高稳定性。
注意,在统计数据分析中, 数据中的随机误差和舍入误差往往超过计算中的误差, 但不稳定的算法会放大这些数据中的误差。 一般只要控制计算中造成的误差比随机误差小一个数量级就可以满足要求。
30.1 Gram-Schmidt正交化方法
定义30.1 (QR分解) 对n \times p矩阵X, 若有n \times p矩阵Q满足Q^T Q = I_p, 以及p阶上三角矩阵R使得 \begin{align} X = QR, \tag{30.4} \end{align} 则称(30.4)为矩阵X的QR分解, 或正交—三角分解。
假设最小二乘问题(30.2)中的X有QR分解X = QR且R满秩。 设Q^* = (Q\,|\,Q_2)是n阶正交阵,则 \begin{align} \| \boldsymbol y - X \boldsymbol\beta \|^2 =& \| (Q^*)^T (\boldsymbol y - Q R \boldsymbol\beta) \|^2 = \| Q^T \boldsymbol y - R \boldsymbol\beta \|^2 + \| Q_2^T \boldsymbol y \|^2, \tag{30.5} \end{align} \| \boldsymbol y - X \boldsymbol\beta \|^2与\| Q^T \boldsymbol y - R \boldsymbol\beta \|^2有相同的最小值点, 用回代法求解 \begin{align} R \boldsymbol\beta = Q^T \boldsymbol y \tag{30.6} \end{align} 就可以得到最小二乘估计\hat{\boldsymbol\beta}, 而(30.5)式右边的\| Q_2^T \boldsymbol y \|^2就是残差平方和。 由于X^T X = R^T R, R^T是下三角形矩阵, 所以从QR分解也可以得到Cholesky分解(可能符号有差别)。
(30.4)实际是把X的各列进行了正交化。 把X的第一列标准化为Q的第一列, 把X的第一、二列正交化和标准化得到Q的第二列, 以此类推,这正是线性代数中的Gram-Schmidt正交化过程, 算法用伪代码表示如下:
QR 分解的Gram-Schmidt算法(RGS) |
---|
令r_{11} \leftarrow \| X_{\cdot 1} \|, Q_{\cdot 1} \leftarrow r_{11}^{-1} X_{\cdot 1}. |
for(j in 2:p) { # 求解Q_{\cdot j}和R的第j列 |
\qquad Q_{\cdot j} \leftarrow X_{\cdot j} |
\qquad for(k in 1:(j-1)){ |
\qquad\qquad r_{kj} \leftarrow X_{\cdot j}^T Q_{\cdot k} |
\qquad\qquad Q_{\cdot j} \leftarrow Q_{\cdot j} - r_{kj} Q_{\cdot k} |
\qquad } |
\qquad r_{jj} \leftarrow \| \boldsymbol Q_{\cdot j} \|, Q_{\cdot j} \leftarrow r_{jj}^{-1} Q_{\cdot j} |
} # end for j |
输出Q和R |
以上的Gram-Schmidt算法(记作RGS)虽然得到了X^T X的QR分解, 但是该算法得到的Q矩阵由于计算误差影响可能会正交性不够好, 使得用这样的到Q矩阵以及(30.6)求解\hat{\boldsymbol\beta}的条件数与直接求解正规方程(30.3)的条件数相同。 对Gram-Schmidt算法略作修改就可以得到更高精度的分解结果。
修正的Gram-Schmidt算法(记作MGS)仅仅调整了计算的次序, 并把Q的分解结果保存在了X的存储空间中。 注意 r_{ij} = Q_{\cdot i}^T X_{\cdot j} = Q_{\cdot i}^T \left( X_{\cdot j} - \sum_{k=1}^{i-1} b_{ik} Q_{\cdot k}\right), \ j=i+1, \dots, p 其中b_{ik}为任意实数。 在MGS的第i步, Q的前i-1列已经求得并保存在X的前i-1列中, 并且R的前i-1行作为Q的前i-1列与X的所有列的内积已经求得, 另外, X的第i到第p列中已经扣除了Q的前i-1列的投影。 所以,第i步只需要求当前保存在X的第i列中的向量长度作为r_{ii}, 然后将当前保存在X第i列中的向量标准化作为Q的第i列, 并计算Q_{\cdot i}(保存在X_{\cdot i}中)与当前X的第i+1, \dots, p列的内积, 作为r_{i, i+1}, r_{i, i+2}, \dots, r_{i,p}的值, 然后从X的第i+1, \dots, p列中扣除Q_{\cdot i}(保存在X_{\cdot i}中)的投影。 所以,算法的第i步将求得Q的第i列(保存在X_{\cdot i}中)和R的第i行, 并把X的第i+1, \dots, p列减去其在Q_{\cdot i}(保存在X_{\cdot i}中)上的投影。 这样,循环到某后续列j'>j时, X的第j'列中该减去的投影就都已经减去了, 只需要标准化第j'列。 算法用伪代码表示如下:
QR分解的修正GS算法(MGS) |
---|
for(i in 1:p) { # 求解Q_{\cdot i}和R的第i行 |
\qquad r_{ii} \leftarrow \| X_{\cdot i} \|, X_{\cdot i} \leftarrow r_{ii}^{-1} X_{\cdot i} |
\qquad for(j in (i+1):p){ |
\qquad\qquad r_{ij} \leftarrow X_{\cdot i}^T X_{\cdot j} |
\qquad\qquad X_{\cdot j} \leftarrow X_{\cdot j} - r_{ij} X_{\cdot i} |
\qquad } |
} # end for i |
输出Q=X和R |
上面的伪代码用了R语言的语法,
但是没有特别针对R语言特点进行定制,
所以伪代码不能作为真正的R语言程序。
在R语言中,for j in (i+1):p
在运行到i=p
时会出错,
因为这时相当于for j in (p+1):p
,
在其他语言中这会执行0次循环,
但R语言会执行j=p+1
和j=p
两步。
在以上的RGS和MGS算法过程中, r_{jj}是X的第j列对X的第1,2,\dots,j-1列作线性回归(无截距项)得到的残差平方和的平方根, 所以某一步中如果r_{jj}=0说明X的第j列与前面的列共线。 在MGS算法过程中, 第j步以后, X的第j+1, j+2, \dots, p列中保存的是原始的X的那些列关于 Q的前j列回归(无截距项)后的残差, 这时类似于解线性方程组时的主元方法, 可以优先选后续列中向量范数最大的列对应的自变量进入模型。
在X列满秩时, RGS和MGS算法得到的上三角阵R的主对角线元素都为正值, 所以X^T X = R^T R是X^T X的Cholesky分解。
如果X不满秩, \text{rank}(X)=r < \min(n,m), 则存在V_{n \times r}, U_{m\times r}, V^T V=I_r, U^T U=I_r, 使得 X = V D U^T, 其中D是n\times m对角阵, 前r个主对角线元素为正,其余元素均为零, 这称为矩阵的奇异值分解, 见特征值和奇异值。
30.2 Householder变换(*
)
Householder变换方法逐次对X进行正交变换把X变成上三角形, 这样就得到了X的QR分解。
设\boldsymbol x = (x_1, x_2, \dots, x_m)^T是一个m维向量, 找一个正交变换把\boldsymbol x最后m-1个元素都变成零。 取U = I_m - d \boldsymbol u \boldsymbol u^T, 其中d = 2 / \boldsymbol u^T \boldsymbol u, 易见对任意的m维非零向量\boldsymbol u,U都是对称的正交方阵(U^T U = U U^T = I_m), 且对任意非零常数b,b \boldsymbol u和\boldsymbol u对应相同的U矩阵。 给定向量\boldsymbol x后,来求\boldsymbol u使得U \boldsymbol x仅有第一个元素非零, 即U \boldsymbol x = s \boldsymbol e_1(\boldsymbol e_1是仅有第一个元素为1,其它元素都等于0的m维向量)。
由于正交变换是使得向量长度不变的, 有\| \boldsymbol x \|^2 = \| U \boldsymbol x \|^2 = s^2 \| \boldsymbol e_1 \|^2 = s^2, 即s^2 = \| \boldsymbol x \|^2 = \boldsymbol x^T \boldsymbol x。 由U \boldsymbol x = s \boldsymbol e_1得 (d \boldsymbol u^T \boldsymbol x) \boldsymbol u = \boldsymbol x - s \boldsymbol e_1, 即\boldsymbol u是\boldsymbol x - s \boldsymbol e_1的数乘结果, 因为\boldsymbol u和\boldsymbol u的非零数乘结果对应相同的矩阵U, 所以不妨取\boldsymbol u = \boldsymbol x - s \boldsymbol e_1, 其中s = \| \boldsymbol x \|。 容易验证, 这样选取的u使得U \boldsymbol x = \| \boldsymbol x \| \boldsymbol e_1。 对向量\boldsymbol x的这种正交变换叫做Householder变换。 因为矩阵\boldsymbol u \boldsymbol u^T秩为1, 这种变换也称为秩1的更新。
可以看出,对任何的m维向量\boldsymbol x都可以用Householder变换把最后m-1个元素都变成零。 对一个n \times p矩阵X(n > p), 可以左乘X_{\cdot 1}的Householder变换阵, 把X_{\cdot 1}的最后n-1个元素变成零; 然后对变换后的X, 左乘一个变换矩阵, 使得X的第1行不变, 而使得X_{\cdot 2}的最后n-2个元素变成零, 即仅对X的第2到n行作n-1维的Householder变换, 而且第一列中的最后n-1个零不变。 以此类推, 在第j次变换时仅对X的第j到n行作n-j+1维的Householder变换。 于是,变换p次后, X变成了一个上三角矩阵R (n \times p的上三角矩阵是指当i>j时总有r_{ij}=0)。
令X^{(0)} = X, \begin{align*} U_j =& \left(\begin{array}{cc} I_{j-1} & \boldsymbol 0 \\ \boldsymbol 0 & U^{(j)} \end{array}\right), \quad X^{(j)} = U_j X^{(j-1)}, \end{align*} 其中U_j是X^{(j-1)}_{\cdot j}的最后n-j+1个元素的Householder变换矩阵。 用\boldsymbol x^{(j)}表示X^{(j-1)}_{\cdot j}的最后n-j+1个元素, 令s_j = \| \boldsymbol x^{(j)} \|, \boldsymbol u^{(j)} = \boldsymbol x^{(j)} - s_j \boldsymbol e_1^{(n-j+1)}, U^{(j)} = I_{n-j+1} - \frac{2}{\| \boldsymbol u^{(j)} \|^2} \boldsymbol u^{(j)} (\boldsymbol u^{(j)})^T, 其中\boldsymbol e_1^{(n-j+1)}表示仅有第一个元素等于1,所有其它元素等于0的n+1-j维向量。 在p步变换后得到 \begin{align*} R^* = U_p U_{p-1} \dots U_1 X \end{align*} 是一个n \times p的上三角阵, 设R^*的前p行组成的矩阵为R, 则R是p阶上三角方阵, 而R^*的后n-p行均为零。 令U^* = U_p U_{p-1} \cdot U_1, 则U^*是一个n阶对称正交阵, 设U^*的前p列组成的n \times p矩阵为U, 则 \begin{align*} X =& U^* R^* = (U \ \ *) \left(\begin{array}{c} R \\ \boldsymbol 0 \end{array}\right) = U R, \end{align*} 于是用p次Householder变换得到了X的QR分解。 用Householder变换作QR分解来求解最小二乘问题是计算精度较高的方法。
编程实现时, 在X的存储空间中保存每次得到的X^{(j)}, 因为第j列的最后n-j个元素为零, 第j个元素为s_j, 所以可以单独保存s_1, s_2, \dots, s_p, 并把\boldsymbol u^{(j)}保存在X的第j列的最后n-j+1个元素的存储空间中。
当X列满秩时,R的主对角线元素都为正值, X^T X = R^T R是X^T X的Cholesky分解。
如果X存在共线问题,比如, X_{\cdot j}可以用X的前j-1列线性表示, 则Householder变换过程进行了j-1次以后, 因为X的前j-1列的最后n-j+1行的元素都变成了零, 所以X的第j列的最后n-j+1个元素也都变成了零, 这时s_j = \| \boldsymbol x^{(j)} \| = 0。 这样在回归计算中可以判断共线是否发生, 也可以用类似选主元的方法, 在第j步看后面的n-j+1列的后n-j+1行组成的子矩阵中哪一列的向量范数最大, 就把那一列对应的自变量与第j个自变量交换位置。
30.3 Givens变换(*
)
Givens变换也是一种正交变换, 是二维向量的旋转变换的推广, 可以把向量的指定元素变成零。
对二维非零向量\boldsymbol x = (x_1, x_2)^T, 把\boldsymbol x看成直角坐标平面上的向量, 作旋转变换把新的x轴旋转到\boldsymbol x的方向上, 设旋转角度为\theta, 则 \begin{aligned} \cos\theta =& \frac{x_1}{\sqrt{x_1^2 + x_2^2}}, \quad \sin\theta = \frac{x_2}{\sqrt{x_1^2 + x_2^2}}, \end{aligned} 旋转变换矩阵和变换结果为 \begin{aligned} R =& \left(\begin{array}{cc} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{array}\right) = (x_1^2 + x_2^2)^{-\frac12} \left(\begin{array}{cc} x_1 & x_2 \\ -x_2 & x_1 \end{array}\right), \\ R \boldsymbol x =& \left(\begin{array}{c} \sqrt{x_1^2 + x_2^2} \\ 0 \end{array}\right). \end{aligned}
对n维向量\boldsymbol x, 为了把x_k变成零, 可以对(x_i, x_k)(i<k)作旋转使得x_k变成零, 保持其它元素不变。 这相当于左乘了一个正交变换矩阵U_{ik}, U_{ik}与单位阵I_n仅在4个元素上有差别: 其(i,i), (i,k), (k,i), (k,k)元素组成的子矩阵恰好为 (x_i, x_k)对应的旋转阵 \begin{align} R_{ik} =& (x_i^2 + x_k^2)^{-\frac12} \left(\begin{array}{cc} x_i & x_k \\ -x_k & x_i \end{array}\right), \end{align} 则U_{ik} \boldsymbol x除去第i, k号元素以外不变, 第i号元素变成\sqrt{x_i^2 + x_k^2}, 第k号元素变成0。
仿照用Householder变换进行QR分解的方法, 利用Givens变换也可以对n\times p矩阵X(n>p)作QR分解。 首先,用n-1次Givens变换可以把X的第1列的最后n-1个元素变成零, 每次用该列主对角线元素与其下面的元素构成Givens变换。 然后,用n-2次Givens变换可以把X的第2列的最后n-2个元素变成零, 每次用该列主对角线元素与其下面的一个元素构成Givens变换。 设已经把X的前j-1列变成了上三角形(前j-1列主对角线下方的元素都变成了零), 在第j步, 用n-j次Givens变换把主对角线下方的元素都变成零, 每次用该列主对角先元素与其下面的一个元素构成Givens变换, 注意此变换仅影响到第j行和第k行(k>j), 所以不会影响到所有列的前j-1行元素, 也就不会影响到前j-1列的上三角部分, 前j-1列已经变成零的那些元素也保持为零, 第j列中已经变成零的元素也保持为零。 注意,第j步关于第j和第k元素作Givens变换时, 仅需对当前X的第j, k两行和第j, j+1, \dots, p列组成的两行的子矩阵作二维的Givens变换即可, 并且对第j列的两个元素的变换结果分别为\sqrt{x_j^2 + x_k^2}和0, 不用重复计算。 如此进行p步后就把X变成了上三角形, 共需要进行 (n-1) + (n-2) + \dots + (n-p) = np - \frac12 p(p+1)次Givens变换。
用Givens变换作QR分解,还可以按照每次消除主对角下方的一行的方式: 首先消除主对角线下方第2行(只有(2,1)号元素), 然后消除主对角线下方第3行(有(3,1,), (3,2)两个元素), 如此重复直到消除了主对角线下方第n行元素(该行除主对角元素以外的所有元素)。 容易看出,按这种次序变换, 在消去第i+1行下三角部分的元素时,不会影响前面的第1,2,\dots,i行中下三角部分已经变成零的元素, 但是会改变第1,2,\dots, i行的上三角部分(包含对角线)的元素。 在回归计算中采用这种变换次序可以适应不断获得新的观测需要更新回归结果的情况(参见 Monahan (2001) §5.8)。
借助于QR分解可以很容易地计算帽子矩阵X(X^T X)^{-1} X^T的主对角线元素h_i的值, 并由此计算多种回归诊断统计量,比如外学生化残差。 回归的一般线性约束的假设检验也可以借助于正交分解的方法计算。 参见(Monahan 2001) §5.9和§5.10。 消去变换(sweep)是解线性方程组时完全消元法(把系数矩阵通过初等变换化为单位阵)的变种, 在原来系数矩阵的存储空间中保存得到的逆矩阵, 为回归变量选择的计算提供了很多便利, 参见 (高惠璇 1995) §5.6 和 (Monahan 2001) §5.12。
在R软件中,
用qr()
函数可以计算矩阵的QR分解。
Julia程序示例见 Julia的Givens变换示例 。
习题
习题4
对线性模型 (30.1), 设X列满秩, 编写用QR分解求\hat{\boldsymbol\beta}、SSE、 回归残差\boldsymbol e = \boldsymbol y - X \boldsymbol \beta和(X^T X)^{-1}的R程序。