第四章 （三）

1.Householder QR分解

Householder变换可以实现任意$m\times n$矩阵$A$的QR分解，其原理是使用变维向量的Householder变换，使得该向量除第一个元素外，其他元素皆变成0。

根据Householder变换一节的分析，欲使一个$p$维向量$x=[x_1,x_2,\dots ,x_p{]}^T$的第1个元素后面的所有元索变为0，则$p$维的Householder向量应取
$$\omega =\frac{x-\beta e_1}{\sqrt{\overline{(\beta )}(\beta -x_1)}}\text{\hspace{1em}(1)}$$
式中
$$\overline{(\beta )}=-|x_1|||x||,\beta =-\frac{x_1}{|x_1|}||x||\text{\hspace{1em}(2)}$$
假定$m\times n$知阵$A$的列分块形式为
$$A_{m\times n}=[a_1,a_2,\dots ,a_n]$$
首先令$x=a_1=[a_{11},a_{21},\dots ,a_{m1}{]}^T$，并取$p=m$，则按照式(1)和式(2)，可以计算得到$u_1={\omega }_m$。此时，
$$H_1=I-u_1{u}_1^T\to A_1=H_{1}A=[a_1^{(1)},a_2^{(1)},\dots ,a_n^{(1)})]$$
变换后，矩阵$A_1$的第1列$a_1^{(1)}$的第一个元素等于$(a_{11}^2+a_{21}^2+\dots +a_{m1}^2{)}^{1/2}$，而该列的其他元素全部为0。
第二步针对矩阵$A_1$的第2列$a_2^{(1)}$，令$p=m-1$和
$$x=[a_{22}^{(1)},a_{32}^{(1)},\dots ,a_{m2}^{(1)}{]}^T$$
又可按照式(1)和式(2)求出(m-1)维向量${\omega }_{m-1}$。 此时，取 u 2 = [ 0 ω m − 1 ] u_2=

u2​=[0ωm−1​​]又可得到
$$H_2=I-u_2{u}_2^T\to A_2=H_2{A}_1=H_2{H}_{1}A=[a_1^{(1)},a_2^{(1)},\dots ,a_n^{(1)}]$$
变换后，矩阵$A_2$的第1列与$A_1$的第1列相同，而第2列$a_1^{(1)}$的第一个元素等于$a_{12}^{(1)}$第二个元素等于$[|a_{22}^{(1)}{|}^2+|a_{32}^{(1)}{|}^2+\dots +|a_{m2}^{(1)}{|}^2{]}^{1/2}$,而该列的其他元素全部为0。
类似地，又可针对矩阵$A_2$的第3列设计Householder变换矩阵$H_3$，使得$A_2$的第一、二个元素保持不变，其他元素组成的$m-2$维向量$x=[a_{33}^{(1)}{|}^2+|a_{43}^{(1)}{|}^2+\dots +|a_{m3}^{(1)}{|}^2]$变换为除第一个元素外的全部元素变为0。

假定矩阵$A$经过$k-1$次Householder变换后，已变成$A^{(k-1)}$，即
$$A^{(k-1)}=H_{k-1}{A}^{(k-2)}=H_{k-1}\dots H_{1}A$$
$$=[a_1^{(k-1)}+a_2^{(k-1)}+\dots +a_n^{(k-1)}],k=2,3,\dots$$
并且其前$k-1$列具有以下变换结果:
$$a_j^{(k-1)}=[a_{1j}^{(k-1)},...,a_{jj}^{(k-1)},0,...,0{]}^T，j=1,2,\dots ,k-1$$
因此，第$k$次Householder变换的目的就是保持前$k-1$列不变，实现$A^{(}k-1)$的第$k$列的下述变换:
H ˉ k [ a k , k ( k − 1 ) a k + 1 , k ( k − 1 ) ⋮ a m , k ( k − 1 ) ] = [ a k , k ( k ) 0 ⋮ 0 ] \bar{H}_k

= Hˉk​⎣ ⎡​ak,k(k−1)​ak+1,k(k−1)​⋮am,k(k−1)​​⎦ ⎤​=⎣ ⎡​ak,k(k)​0⋮0​⎦ ⎤​
这相当于对矩阵$A^{(k-1)}$进行Householder变换$H_k{A}^{(k-1)}$时取
H k = [ I k − 1 0 0 H ~ k ] H_k=

[Ik−100H~k]" role="presentation" style="position: relative;">[Ik−100H~k]$$\left[\begin{array}{cc}{I}_{k-1}& 0\\ 0& {\tilde{H}}_k\end{array}\right]$$

Hk​=[Ik−1​0​0H~k​​]
$n$次Householder变换后，即可实现QR分解。

例
已知线性方程组
$x_1+2x_2=5$
$2x_1+3x_2=-2$
$6x_1+7x_2=1$
用Householder变换求上述方程组的最小二乘解。

解
记增广矩阵
A = [ 1 2 5 2 3 − 2 6 7 1 ] A=

A=⎣ ⎡​126​237​5−21​⎦ ⎤​
针对第1列$[1,2,6{]}^T$，由式(1)和式(2)求得Householder向量$u_1=[1.075,0.290,0.871{]}^T$。经过Householder矩阵$H_1=I-u_1{u}_1^T$变换后，得
H 1 A = [ − 6.415 − 7.826 − 23.008 0 0.592 0.875 0 0.808 1.438 ] H_1A=

[−6.415−7.826−23.00800.5920.87500.8081.438]" role="presentation" style="position: relative;">⎡⎣⎢−6.41500−7.8260.5920.808−23.0080.8751.438⎤⎦⎥$$\left[\begin{array}{ccc}-6.415& -7.826& -23.008\\ 0& 0.592& 0.875\\ 0& 0.808& 1.438\end{array}\right]$$

H1​A=⎣ ⎡​−6.41500​−7.8260.5920.808​−23.0080.8751.438​⎦ ⎤​
针对上述矩阵的第2列，计算Householder向量$u_2=[0,1.161,-0.809{]}^T$。经过House holder变换$H_2=I-u_2{u}_2^T$后，即得QR分解为
H 2 H 1 A = [ − 6.415 − 7.826 − 23.008 0 − 1.105 − 1.851 0 0 − 0.160 ] H_2H_1A=

[−6.415−7.826−23.0080−1.105−1.85100−0.160]" role="presentation" style="position: relative;">⎡⎣⎢−6.41500−7.826−1.1050−23.008−1.851−0.160⎤⎦⎥$$\left[\begin{array}{ccc}-6.415& -7.826& -23.008\\ 0& -1.105& -1.851\\ 0& 0& -0.160\end{array}\right]$$

H2​H1​A=⎣ ⎡​−6.41500​−7.826−1.1050​−23.008−1.851−0.160​⎦ ⎤​
由此得线性方程组
$-6.415x_1-7.826x_2=-23.008$
$-1.105x_2=-1.851$
用回代法解之，得$x_2=1.675$和$x_1=1.543$。将这两个解代入原方程，得三个方程的拟合误差分别为$e_1=0.107,e_2=0.111和e_3=0.017$。

2.带有列旋转的Householder OR分解算法

如果矩阵$A_{m\times n}$的秩小于$n$，则上述Householder QR分解不一定能够产生正交矩阵$Q$。为了解决这个问题，Golub与Van Loan提出使用带有列旋转的Householder OR分解算法。该算法计算分解
Q T A Π = [ R 11 R 12 O O ] Q^TA\Pi=

QTAΠ=[R11​O​R12​O​]
式中，$Q$为正交矩阵;$R_{11}$为$r\times r$上三角矩阵，且非奇异;$r=rank(A)$;矩阵$\Pi$为置换矩阵，矩阵乘积$A\Pi$表示对矩阵A的列互换(称为列旋转)。

列旋转的原则是:在对矩阵$A$的第$k$列进行Householder变换之前，先比较$A$的第$k$列至第$n$列中各列的范数，并把范数最大、且与第$k$列相距最近的列记为$piv(k)$，然后将第$piv(k)$列与第$k$列相互交换位置。这相当于第$k$次列旋转矩阵$\Pi$由$n\times n$单位矩阵交换第$piv(k)$和第$k$列而成。

算法 具有列旋转的HouseholderQR分解算法

for j=1:n
	c(j)=A(1:m,j)^TA(1: m,j) 
end
r=0;T=max{c(1),c(2),...,c(n)}
求满足c(k)=T的最小整数k
whilet T>0
	r=r+1
	piv(r)=k:A(1:m,r)↔A(1:m,k);c(r)↔c(k)
	v(r:m)=house(A(r:m,r))
	A(r:m,r:n)=row.house(A(r:m,r:n),v(r:m)) 
	A(r+1:m,r)=v(r+1:m)
	for i=r+1:n
		c(i)=c(i)- A(r,i)^2 
	end
	if r<n
		T=max{c(r+1),c(r+2),...,c(n)}
		在k=r+1,r+2,...,n中，求满足c(k)=T的最小整数k 
	else
		T=0
	end 
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

上述算法中，$v=house(x)$代表针对向量x的Householder向量函数;而$A=row.house(A,v)$为Householder左乘的数，表示用Householder变换结果$PA=(I-2v{v}^T/v^{T}v)A$重写矩阵$A$。

具有列旋转的HouseholderQR分解在矩阵的奇异值分解与特征值分解中有着重要的应用。
虽然我们分别叙述了修正Gram-Schmidt法和Householder法两种QR分解方法，但是这两种方法在数学上是等价的。除了存在$Q$的第$j$列和$R$的第$j$行的符号可同时反转的自由度以外，修正Gram-Schmidt法和Householder法是相同的。