(仅供参考)矩阵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)']';