数值分析中的QR分解及其代码实现

QR分解

若 A∈Cm×k A ∈ C m × k A\in C^{m×k} 是一个列满秩的矩阵，rank(A) = k，则矩阵A 可以分解为 A=QR A = Q R A = QR ,

Q∈Cm×k Q ∈ C m × k Q \in C^{m×k} ，Q 的列向量为A 的列空间的标准正交基， R∈Ck×k R ∈ C k × k R\in C^{k×k} ，是一个可逆的上三角矩阵，

A 的列向量线性无关， A=(α1,α2,...αk) A = ( α 1 , α 2 , . . . α k ) A= (\alpha_1, \alpha_2, ...\alpha_k) ，将这k个列向量进行Schmidt正交化，得到A 的列向量空间的标准正交基，

正交化： β1=α1 β 1 = α 1 \beta_1 = \alpha_1

β2=α2−(α2,β1)(β1,β1)β1 β 2 = α 2 − ( α 2 , β 1 ) ( β 1 , β 1 ) β 1 \beta_2 = \alpha_2 -\frac{(\alpha_2, \beta_1)}{(\beta_1, \beta_1)} \beta_1

βk=αk−∑k−1i=1(αk,βi)(βi,βi)βi β k = α k − ∑ i = 1 k − 1 ( α k , β i ) ( β i , β i ) β i \beta_k = \alpha_k - \sum_{i=1}^{k-1} {\frac{(\alpha_k, \beta_i)}{(\beta_i, \beta_i)} \beta_i}

标准化： ϵi=βi∥βi∥ ϵ i = β i ‖ β i ‖ \epsilon_i = \frac{\beta_i}{\Vert \beta_i \Vert}

把正交化方法和标准化方法结合在一起：
β1=α1 β 1 = α 1 \beta_1 = \alpha_1 , ϵ1=β1∥β1∥ ϵ 1 = β 1 ‖ β 1 ‖ \epsilon_1 = \frac {\beta_1}{\Vert \beta_1 \Vert}

β2=α2−(α2,ϵ1) β 2 = α 2 − ( α 2 , ϵ 1 ) \beta_2 = \alpha_2 - (\alpha_2, \epsilon_1) , ϵ2=β2∥β2∥ ϵ 2 = β 2 ‖ β 2 ‖ \epsilon_2 = \frac{\beta_2}{\Vert \beta_2 \Vert}

…..

βn=αn−∑n−1i=1(αn,ϵi)ϵi β n = α n − ∑ i = 1 n − 1 ( α n , ϵ i ) ϵ i \beta_n = \alpha_n - \sum_{i = 1} ^{n-1} (\alpha_n , \epsilon_i) \epsilon_i , ϵn=βn∥βn∥ ϵ n = β n ‖ β n ‖ \epsilon_n = \frac {\beta_n}{\Vert \beta_n \Vert}

得到 Q(ϵ1,ϵ2...ϵk) Q ( ϵ 1 , ϵ 2 . . . ϵ k ) Q(\epsilon_1,\epsilon_2 ...\epsilon_k) ，由 A=QR A = Q R A = QR 可知， R=Q−1A=QTA R = Q − 1 A = Q T A R = Q^{-1} A = Q^T A ，得到R 矩阵为上三角矩阵，对角线主元素 rii=∥βi∥>0 r i i = ‖ β i ‖ > 0 r_{ii} = \Vert \beta_i \Vert > 0 ，正定矩阵，

当 A 为普通矩阵时， A=QR A = Q R A = QR ，为QR分解；

当 A∈Cn×n A ∈ C n × n A \in C^{n×n} 为可逆方阵时， A=UR A = U R A = UR ，为UR分解。

QR方法是求一般矩阵全部特征值的最有效并广泛应用的方法之一，在应用中，先把一般的矩阵经过正交相似变换化成上Hessenberg矩阵（拟上三角矩阵），再应用QR方法求其特征值和特征向量。

Householder矩阵

Householder矩阵为反射变换或镜像变换的变换矩阵，将向量x映射为关于“与单位向量u正交的n-1维子空间”对称的向量y的镜像变换定义如下：

y=x−2u(uTx)=(I−2uuT)x=Hx y = x − 2 u ( u T x ) = ( I − 2 u u T ) x = H x y = x- 2 u(u^T x) =(I-2 u u^T) x = H x

设单位向量 u∈Rn u ∈ R n u \in R^n ，称 H=I−2uuT H = I − 2 u u T H = I -2 u u^T ，为Householder矩阵（初等反射矩阵），

性质：对称矩阵，正交矩阵，对合矩阵，

Hessenberg矩阵

若是矩阵 A∈Rn×n A ∈ R n × n A \in R^{n×n} 的次对角线下的元素均为零，即 i >j+1 时， aij=0 a i j = 0 a_{i j} = 0 ，这样的矩阵称为 Hessenberg矩阵（拟上三角矩阵）

a11a21a12a22a32......a33..a1n−1a2n−1.....ann−1a1na2na3n..ann a 1 1 a 1 2 . . . a 1 n − 1 a 1 n a 2 1 a 2 2 . . . a 2 n − 1 a 2 n a 3 2 a 3 3 . . . a 3 n . . . . . . a n n − 1 a n n \begin{matrix} a_{1 1} & a_{1 2} & ...& a_{1 n-1} & a_{1 n} \\ a_{2 1 }& a_{2 2} & ... & a_{2 n-1}& a_{2 n} \\ & a_{3 2 } & a_{3 3}& ...& a_{3 n} \\ & \\ &&& a_{n n-1} & a_{n n} \end{matrix}

对任意矩阵A ，总存在正交阵Q 使得 Q^{-1}A Q 为上Hessenberg 矩阵。Householder矩阵即为这样的一个正交矩阵。

Givens矩阵

一般的，在n维空间 Rn R n R^n 中，令 e1,e2,....en e 1 , e 2 , . . . . e n e_1 , e_2 , ....e_n 为标准正交基，于是在平面 [ei,ej] [ e i , e j ] [ e_i, e_j] 中的旋转变换定义如下：

Ti,j=1....1c−s1....1sc1....1 T i , j = 1 . . . . 1 c s 1 . . . . 1 − s c 1 . . . . 1 T_{i ,j } =\begin{matrix}1 \\ & .. \\ && .. \\ &&& 1 \\ &&&& c &&&&& s\\ &&&&&1 \\&&&&&& ..\\ &&&&&&& .. \\&&&& & &&&1 \\ &&&& -s &&&&& c \\ &&&&&&&&&&1 \\&&&&&&&&&&&.. \\&&&&&&&&&&&&..\\ &&&&&&&&&&&&&1 \end{matrix}

其中c , s 满足 c2+s2=1 c 2 + s 2 = 1 c^2 +s^2 = 1 ，行数分别为i , j ，也可记为 Ti,j=Ti,j(c,s) T i , j = T i , j ( c , s ) T_{i, j} = T_{i ,j}(c, s) ，由Givens矩阵确定的线性变换称为Givens变换（初等旋转变换）。存在角度 θ θ \theta ，使得 c=cosθ，s=sinθ c = c o s θ ， s = s i n θ c = cos \theta ， s = sin \theta

任何n阶实非奇异矩阵 A=ai,j A = a i , j A = { a_{i,j} }可通过左连乘 Givens 矩阵化为上三角矩阵 TA=Y T A = Y T A = Y ；Y 为上三角矩阵

QR分解求特征值和特征向量的工程实现方法：

1）通过正交相似变换将矩阵A 化为Hessenberg矩阵 B ，其变换矩阵可用Householder矩阵H ，H矩阵均为对称正交矩阵。

A(2)=H1AH1 A ( 2 ) = H 1 A H 1 A^{(2)} = H_1 A H_1, A(3)=H2A(2)H2 A ( 3 ) = H 2 A ( 2 ) H 2 A^{(3)} = H_2 A^{(2)} H_2 …… A(n−1)=Hn−2...H2H1AH1A2...Hn−2 A ( n − 1 ) = H n − 2 . . . H 2 H 1 A H 1 A 2 . . . H n − 2 A^{(n-1)} =H_{n-2} ...H_2 H_1 A H_1 A_2 ... H_{n-2}

保持A 序列变换为相似变换，所以右乘 Hi H i H_i ，之所以采用Householder矩阵变换而不采用Givens矩阵，是因为可以减少运算次数。

记 P=H1H2...Hn−2 P = H 1 H 2 . . . H n − 2 P = H_1 H_2...H_{n-2} ，所以有 A(n−1)=PTAP A ( n − 1 ) = P T A P A^{(n-1)} = P^T A P ，因而 A(n−1) A ( n − 1 ) A^{(n-1)} 与A 相似。

A(n−1) A ( n − 1 ) A^{(n-1)} 为上Hessenberg 矩阵，记为矩阵 B ，

2）将上Hessenberg 矩阵B 进行QR 分解，运用 Givens 旋转变换方法将矩阵变换为上三角矩阵 R ，Givens矩阵 Ti,j T i , j T_{i,j} 连乘即为标准正交矩阵 Q ，令 B1=B=Q1R1 B 1 = B = Q 1 R 1 B_1 =B = Q_1 R_1 ，QR分解得到： B1=Q1R1 B 1 = Q 1 R 1 B_1 = Q_1 R_1 ，

让上式右端逆序相乘，矩阵相乘令： B2=R1Q1 B 2 = R 1 Q 1 B_2 = R_1 Q_1 ，

又对 B2 B 2 B_2 ，进行QR分解，得： B3=Q2R2 B 3 = Q 2 R 2 B_3 = Q_2 R_2 ，

反复进行这种迭代运算，得到一个矩阵序列 Bk B k { B_k }，

由QR分解 Bk=QkRk B k = Q k R k B_k = Q_k R_k ，得到 Rk=QTkBk R k = Q k T B k R_k = Q_k^{T} B_k

又由矩阵相乘 Bk+1=RkQk B k + 1 = R k Q k B_{k+1} = R_k Q_k ，得到 Bk+1=QTkBkQk B k + 1 = Q k T B k Q k B_{k+1} = Q_k^{T} B_k Q_k , k=1,2,... k = 1 , 2 , . . . k = 1,2, ...
因而 Bk+1 B k + 1 B_{k+1} 与 Bk B k B_k 相似，所有的矩阵 Bk B k B_k 都与矩阵 B B B 相似，因而也与原矩阵A 相似，它们都有相同的特征值。在迭代过程中，迭代序列Bk" role="presentation" style="position: relative;">BkBk{ B_k } 收敛与上三角矩阵，上三角矩阵的主对角线元素就是矩阵 A 的特征值，矩阵 Q 即为特征值对应的特征向量。

matlab实现QR分解代码如下：

1.构造Householder正交相似变换矩阵

function H = househ(x)
%功能：对于向量x，构造Householder变换矩阵H，
%     使得Hx = (*, 0, 0,......, 0)',其中|*|=norm(x, 2)
%输入：列向量x
%输出：Householder变换矩阵Hn = length(x);               % x为矩阵A的对角线下部的列向量 x = A(k+1: n, k)I = eye(n);sn = sign(x(1));            %判断x向量第一个元素的符号if(sn == 0)sn=1;endz = x(2:n);if(norm(z,inf) == 0)H = I;return;endsigma = -sn*norm(x);                 %  sigma = -sign(x1) ||x||infu = x;u(1) = u(1)-sigma;                 %u1 = u1 -sigma* 1; u向量中第一个元素值lo = sigma*(sigma-x(1));          %lo 为转换比例因子，防止计算过程中的溢出和误差累积H = I-u*u'/lo;

u(1) 为 x(1) 与 norm(x) 即向量 x 的范数之间的差值，最终得到Householder变换矩阵；

2.得到上Hessenberg矩阵

function A = hessen(A)
%功能：用Householder变换化矩阵A为上Hessenberg型
%输入：n阶实方阵A
%输出：A的Hessenberg型
[n, n] =size(A);
for k =1: n-2x = A(k+1: n, k);  % x 取为矩阵A 每列对角线下的元素组成的向量H =househ(x);      % 得到其Householder变换矩阵，A(k+1:n, k:n) = H*A(k+1:n, k:n);     A(1:n, k+1:n) = A(1:n, k+1:n)*H;  % A^(n-1) = HAH,经过 n-2次的变换,矩阵A变为Hessen矩阵
end

3.基于Givens旋转变换的QR分解

function A = qrtran(m, A)
%功能：对A的左上角的m阶对角块做QR变换；
%      先用Givens变换作QR分解 A=QR，再做相似变换A：= Q'AQ = RQ
%输入：n阶Hessenberg型矩阵A，其中A(m+1，m) = 0, m>2
%输出：变换后的Hessenberg型矩阵A
Q = eye(m);
for i =1:m-1xi = A(i, i);xk = A(i+1, i);if(xk ~=0)d = sqrt(xi^2 +xk^2);c = xi/d;s = xk/d;J = [c s; -s c];A(i:i+1, i:m) = J*A(i:i+1, i:m);
% J为每次的2阶Givens矩阵，J*A的两行矩阵，得到上三角矩阵R的两行Q(1:m, i:i+1) = Q(1:m, i:i+1)*J';
%将所有的Givens矩阵连乘起来，得到正交单位矩阵Qend
end
A(1:m, 1:m) = A(1:m, 1:m)*Q;
% A(n+1) = R(n)*Q(n)，矩阵相乘得到下一个迭代矩阵A

正定矩阵

广义定义：设M是n阶方阵，如果对任何非零向量z，都有 zTMz>0 z T M z > 0 z^TMz > 0，其中 zT z T z^T 表示 z 的转置，就称M为正定矩阵。

狭义定义：一个n阶的实对称矩阵M是正定的的条件是当且仅当对于所有的非零实系数向量 z ，都有 zTMz z T M z z^T M z> 0。其中 zT z T z^T表示 z 的转置。

等价命题 :

对于n阶实对称矩阵A，下列条件是等价的：

（1）A是正定矩阵；

（2）A的一切顺序主子式均为正；

（3）A的一切主子式均为正；

（4）A的特征值均为正；

（5）存在实可逆矩阵C，使 A=CTC A = C T C A=C^TC ；

（6）存在秩为n的m×n实矩阵B，使 A=BTB A = B T B A=B^TB ；

（7）存在主对角线元素全为正的实三角矩阵R，使 A=RTR A = R T R A=R^TR