(仅供参考)矩阵QR分解MATLAB自编程序

3矩阵的QR分解
3.1 基于Schmidt正交化的QR分解
3.1.1原理
Schmidt正交化方法是矩阵的QR分解最常用的方法.主要依据下面的两个结论:
结论1设A是n阶实非奇异矩阵,则存在正交矩阵Q和实非奇异上三角矩阵R使A有QR分解;且除去相差一个对角元素的绝对值（模）全等于1的对角矩阵因子外,分解是唯一的.
结论 2 设A是 m × n实矩阵,且其n个列向量线性无关,则 A有分解A =QR,其中Q是m ×n实矩阵,且满足Q H TQ = E,R是n阶实非奇异上三角矩阵该分解除去相差一个对角元素的绝对值（模）全等于1的对角矩阵因子外是唯一的.
用Schmidt正交化分解方法对矩阵进行QR分解时，所论矩阵必须是列满秩矩阵。

为方便与后续基于Householder变换的QR分解法对比，这里以方阵为例.
3.1.2算法
1写出矩阵的列向量;
2列向量按照Schmidt正交化正交;
3得出矩阵的Q′,R′;
4对Q′的列向量单位化得到Q，R′的每行乘Q′每列的模得R.
3.1.3流程图
3.1.4程序
function [X,Q,R] = QRDecomsch(A,b)
%方阵的QR的Gram-Schmidt正交化分解法，并用于求解AX=b方程组[m,n]=size(A);
if m~=n
disp('不满足QR分解要求');
end
Q=zeros(m,n);
X=zeros(n,1);
R=zeros(n);
for k=1:n
R(k,k)=norm(A(:,k));
if R(k,k)==0
break;
end
Q(:,k)=A(:,k)/R(k,k);
for j=k+1:n
R(k,j)=Q(:,k)'*A(:,j);
A(:,j)=A(:,j)-R(k,j)*Q(:,k);
end
end
if nargin==2
b=Q'* b;
X(n)=b(n)/R(n,n);
for i=n-1:-1:1
X(i)=(b(i)-sum(R(i,i+1:n).*X(i+1:n)'))/R(i,i);
end
else
X=[];
end
3.2 基于Householder变换的QR分解
3.2.1原理
设A为任一n阶方阵，则必存在n阶酉矩阵Q和n阶上三角阵R，使得A=QR
设w∈C n是一个单位向量，令2H
=-
H Iωω
则称H 是一个Householder 矩阵或Householder 变换。

则对于任意的C n 存在Householder 矩阵H ，使得Hx=-au 。

其中 3.2.2算法
第一步，将矩阵A 按列分块写成A=(α1，α2，…，αn ）.如果 α1≠0，则可得，存在n 阶householder 矩阵H 1使得
H 1α1=-a 1e 1，|a 1|=||α1||，e 1∈C n
于是有H 1A=（H 1α1，H 1α2，…，H 1αn ）= 11*0
n a A --⎛⎫
⎪⎝⎭
如果α1=0，则直接进行下一步，此时相当于取H 1=I n ，而a 1=0. 第二步，将矩阵A n-1按列分块写成A n-1=（αi ,α2，… ,αn-1）。

如果α1≠0，则可得，存在n-1阶householder 矩阵H ’2使得H ’2α2=-a 2e 1，| a 2 |=||α2||，e 1∈C n
于是有H ’2 A n-1=（H ’2α2，…，H ’2αn-1）=2
2*0
n a A --⎛⎫
⎪⎝⎭
此时，令H 2=2100'T H ⎛⎫
⎪⎝⎭
则H 2是n 阶Householder 矩阵，且使
H 2H 1A=1
22
**0
*0n a a A --⎛⎫ ⎪- ⎪ ⎪⎝
⎭
如果α1=0，则直接进行下一步。

第三步，对n-2阶矩阵继续进行类似的变换，如此下去，之多在第n-1步，我们可以找到Householder 矩阵H 1，H 2，…，H n-1使得
a x =
H n-1…H2H1A=
1
2
***
0** 00*
000'
nn a
a
α
-⎛⎫ ⎪
-
⎪ ⎪
⋅⋅⋅
⎪⎝⎭
令Q= H n-1…H2H1，则Q是酉矩阵之积，从而必有酉矩阵并且A=QR
3.2.3流程图
3.2.4程序
function [ X,Q,R ] = QRDecomhouse(A,b)
%用Householder变换将方阵A分解为正交Q与上三角矩阵R的乘积，并用于求解AX=b方程组
[n,n]=size(A);
E=eye(n);
X=zeros(n,1);
R=zeros(n);
P1=E;
for k=1:n-1
%构造w,使Pk=I-2ww'
s=-sign(A(k,k))* norm(A(k:n,k));
R(k,k)=-s;
if k==1
w=[A(1,1)+s,A(2:n,k)']';
else
w=[zeros(1,k-1),A(k,k)+s,A(k+1:n,k)']';
R(1:k-1,k)=A(1:k-1,k);
end
if norm(w)~=0
w=w/norm(w);
end
P=E-2*w*w';
A=P*A;
P1=P*P1;
R(1:n,n)=A(1:n,n);
end
Q=P1';
if nargin==2
b=P1*b;
X(n)=b(n)/R(n,n);
for i=n-1:-1:1
X(i)=(b(i)-sum(R(i,i+1:n).*X(i+1:n)'))/R(i,i);
end
else
X=[];
end
3.3 运行验证与分析
A=
1 2
1 0
1−1
2 1
1−1
−1 1
1 2
−3 1
，b=[1 0 1 1]T，求A得QR分解，并解方程组
A x=b.
3.3.1 MATLAB自带QR分解函数计算
在MATLAB中调用：[Q0,R0]=qr(A)，运行结果如下：>> clear all
>> A=[1,2,1,-1;1,0,2,1;1,-1,1,2;-1,1,-3,1];
>> b=[1 0 1 1]';
>> X0=inv(A)*b
>> [Q0,R0]=qr(A)
Q0 =
-0.5000 0.8165 -0.0577 -0.2828
-0.5000 -0.0000 0.1732 0.8485
-0.5000 -0.4082 -0.7506 -0.1414
0.5000 0.4082 -0.6351 0.4243
R0 =
-2.0000 0 -3.5000 -0.5000
0 2.4495 -0.8165 -1.2247
0 0 1.4434 -1.9053
0 0 0 1.2728
X0 =
2.0000
-0.0000
-1.0000
3.3.2 Schmidt正交化计算结果
>> clear all
>> A=[1,2,1,-1;1,0,2,1;1,-1,1,2;-1,1,-3,1];
>> b=[1 0 1 1]';
>> [X1,Q1,R1] = QRDecomsch(A,b)
X1 =
2.0000
0.0000
-1.0000
0.0000
Q1 =
0.5000 0.8165 -0.0577 -0.2828
0.5000 0 0.1732 0.8485
0.5000 -0.4082 -0.7506 -0.1414
-0.5000 0.4082 -0.6351 0.4243 R1 =
2.0000 0
3.5000 0.5000
0 2.4495 -0.8165 -1.2247
0 0 1.4434 -1.9053
0 0 0 1.2728 3.3.3 Householder变换运行结果
>> clear all
>> A=[1,2,1,-1;1,0,2,1;1,-1,1,2;-1,1,-3,1];
>> b=[1 0 1 1]';
>> [X2,Q2,R2] = QRDecomhouse(A,b)
X2 =
2.0000
0.0000
-1.0000
0.0000
Q2 =
0.5000 0.8165 0.0577 -0.2828
0.5000 0.0000 -0.1732 0.8485
0.5000 -0.4082 0.7506 -0.1414
-0.5000 0.4082 0.6351 0.4243 R2 =
2.0000 0
3.5000 0.5000
0 2.4495 -0.8165 -1.2247
0 0 -1.4434 1.9053
0 0 0 1.2728 3.3.4 数据分析
由1、2、3可知，Gram-Schmidt正交化法与Householder变换法所得Q、R以及方程组解与MATLAB自带函数结果一致，可见程序运行结果正确。