QR方法是做什么的?它可以求一般矩阵的全部特征值和特征向量。为了使用这种方法我们首先要知道一些理论基础。
Householder变换
(课本上的话虽然严谨,但是枯燥。总是感觉大学课本上的东西是给会的人看的,这里我用自己语言聊一下)什么叫Householder变换呢?就是在1932年,Householder这个人发现了一种变换。他给出了这样的公式:
(1)H是对称矩阵。直接写出来就可以看出来了。
(2)H是正交矩阵, ,可以直接证出来。
(3) ,
(4) 这个性质是我们直接用的: ,存在反射矩阵 ,使得 ,这个时候的 .
OK有了这个性质就可以做我们想做的事了。
你看求A的所有特征值,我们最容易想到的就是把A变成上对角矩阵,这样我们就可以直接看出来,对角矩阵上的对角元素就是特征值。怎么变这里就是用了上面我BB的理论。哎,以为A是矩阵,Householder变换是针对向量的,所以这里就用的约化定理(不想写不写了。)再写一个比较直接的式子。
QR方法的迭代格式
好了上代码:
function [Q,R] = Household(A)
[m,~] = size(A);
Q = eye(m);
for i = 1:m-1
y = zeros(m,1);
I = eye(m);
y(1:i,1) = A(1:i,i);
y(i,1) = -sign(A(i,i))*(sum(A(i:end,i).^2)^(1/2));
x = A(:,i);
u = x - y;
w = u/norm(u,2);
H = I - 2*w*w';
A = H*A;
Q = H*Q;
end
R = A;
Q = Q';
end
function demo_five
A = [5 -3 2;6 -4 4;4 -4 5];
for i =1:100 %这里表示迭代次数。
[Q,R] = Household(A);
A = R*Q;
end
diag(A)
end
算法改进
为了加快收敛速度,我们先利用Householder变换把A变换为Hessenberg矩阵(矩阵的下次对角线下方元素均为零),
代码:
function A_hessenberg = Hessenberg(A)
[m,~] = size(A);
for i=1:m-2
x = A(:,i);
y = zeros(m,1);
y(1:i,1) = A(1:i,i);
y(i+1,1) = -sign(A(i+1,i))*(sum(A(i+1:end,i).^2)^(1/2));
w = (x-y)/norm(x - y);
H = eye(m) - 2*(w*w');
A = H*A*H;
end
A_hessenberg = A;
end
利用旋转G求出Q,(这里的作用是:把Hessenberg矩阵的下次对角线全部化为0这样就求出Q了).
function Q = Givens(A)
[m,~] = size(A);
Q = eye(m);
for i = 1:m-1
R = eye(m);
r = sqrt(A(i,i)^2 + A(i+1,i)^2);
c = A(i,i)/r;
s = A(i+1,i)/r;
R(i,i) = c;
R(i+1,i+1) = c;
R(i,i+1) = s;
R(i+1,i) = -s;
A = R*A;
Q = Q*R';
end
end
主函数
function [lamda,vector] = demo_five_1
A = [5 -3 2;6 -4 4;4 -4 5];
A_h = Hessenberg(A);
for i =1:30
Q = Givens(A_h);
A_h = Q'* A_h *Q;
end
lamda = diag(A_h);
vector = fan_mi_fa(A,lamda);
end
求特征向量这里我用的是反幂法。(这个方法的重点是理解: 的迭代)
function vector = fan_mi_fa(A,lamda)
[m,~] = size(A);
vector = ones(m,m);
for i = 1:m
p = lamda(i);
x = ones(m,1);
for j = 1:50
if max(abs(x)) == max(x)
max_x = max(abs(x));
else
max_x = -max(abs(x));
end
y = x / max_x;
x = inv( A - p*eye(m) ) * y;
end
vector(:,i) = y;
end
end
好了写完了,哎其实中间还有好多过程没有搞懂(希望以后搞懂了,继续添加),但是代码是自己写的绝对靠谱。
我写博客的目的有两个:
(1)把自己的学习资料保存一下,以后用到的话方便查询,之前写在纸上的资料,大多数都被自己搞丢了。
(2)希望自己有所收获,提升。如果能交到朋友更好了。