基于Givens变换的QR分解

01.function [Q,R]=qrgv(A)

02.% 基于Givens变换，将方阵A分解为A=QR，其中Q为正交矩阵，R为上三角阵

03.%

04.% 参数说明

05.% A：需要进行QR分解的方阵

06.% Q：分解得到的正交矩阵

07.% R：分解得到的上三角阵

08.%

09.% 实例说明

10.% A=[-12 3 3;3 1 -2;3 -2 7];

11.% [Q,R]=qr(A) % 调用MATLAB自带的QR分解函数进行验证

12.% [q,r]=qrgv(A) % 调用本函数进行QR分解

13.% q*r-A % 验证 A=QR

14.% q'*q % 验证q的正交性

15.% norm(q) % 验证q的标准化，即二范数等于1

16.%

17.% 线性代数基础知识

18.% 1.B=P*A*inv(P)，称A与B相似，相似矩阵具有相同的特征值

19.% 2.Q*Q'=I，称Q为正交矩阵，正交矩阵的乘积仍为正交矩阵

20.%

21.% by dynamic of Matlab技术论坛

22.% see also http://www.matlabsky.com

23.% contact me matlabsky@gmail.com

24.% 2010-01-17 22:51:18

25.%

26.n=size(A,1);

27.R=A;

28.Q=eye(n);

29.for i=1:n-1

30.    for j=2:n-i+1

31.        x=R(i:n,i);

32.        rt=givens(x,1,j);

33.        r=blkdiag(eye(i-1),rt);

34.        Q=Q*r';

35.        R=r*R;

36.    end

37.end

38.

39.function [R,y]=givens(x,i,j)

40.% 求解标准正交的Given变换矩阵R，使用Rx=y，其中y(j)=0，y(i)=sqrt(x(i)^2+x(j)^2)

41.%

42.% 参数说明

43.% x：需要进行Givens变换的列向量

44.% i：变为sqrt(x(i)^2+x(j)^2)的元素下标

45.% j：变为0的元素的下标

46.% R：Givens变换矩阵

47.% y：Givens变换结果

48.%

49.% 实例说明

50.% x=[1 3 5 9 6]'; % 将3等效到9上

51.% [R,y]=givens(x,4,2) % 注意3的下标为2，9的下标为4

52.% R*x-y % 验证Rx=y

53.% R'*R % 验证正交性

54.% norm(R) % 验证标准性，就是范数为1

55.%

56.% 关于Givens变换说明

57.% 1.Givens矩阵是标准正交矩阵，也叫平面旋转矩阵，它是通过坐标旋转的原理将元素j的数值等效到元素i上

58.% 2.Givens变换每次只能将一个元素变为0，而Householder变换则一次可以将任意个元素变为0

59.% 3.Givens变换常用于将矩阵A变为对角阵

60.%

61.xi=x(i);

62.xj=x(j);

63.r=sqrt(xi^2+xj^2);

64.cost=xi/r;

65.sint=xj/r;

66.R=eye(length(x));

67.R(i,i)=cost;

68.R(i,j)=sint;

69.R(j,i)=-sint;

70.R(j,j)=cost;

71.y=x(:);

72.y([i,j])=[r,0];
 

文章出处：<a target=_blank href="http://www.matlabsky.com/thread-4850-1-1.html">http://www.matlabsky.com/thread-4850-1-1.html</a>