matlab实现上Hessenberg阵的QR算法

实现Hessenberg阵的QR算法可以按照以下步骤进行：

1. 初始化Hessenberg阵H，令Q为单位阵，i=1.

2. 当i<=n-2时，进行以下步骤：

a. 计算H的第i列到第n列的Householder矩阵，得到阵H'=Q'HQ.

b. 计算H'的第(i+1)行第i列元素的Givens旋转矩阵，令H''=RH'，其中R为Givens旋转矩阵.

c. 计算H''的第i行到第n行的Householder矩阵，得到阵H'''=Q''H''Q''，其中Q''为第二步所得到的Q与R的乘积.

d. 令H=H'''，Q=Q''.

e. i=i+1.

3. 对于H的最后两列进行双重隐式QR迭代，直到收敛.

Matlab代码如下：

function [Q,H] = hessqr(H)
% Hessenberg阵的QR分解
% 输入：H为Hessenberg阵
% 输出：Q为正交阵，H为上Hessenberg阵的上三角阵

n = size(H,1);
Q = eye(n);
i = 1;

while i <= n-2
% 计算第i列到第n列的Householder矩阵
[v,beta] = house(H(i+1:n,i));
H(i+1:n,i:n) = H(i+1:n,i:n) - beta*v*(v'*H(i+1:n,i:n));
H(:,i+1:n) = H(:,i+1:n) - (H(:,i+1:n)*v)*beta*v';
% 计算第(i+1)行第i列元素的Givens旋转矩阵
[c,s] = givens(H(i,i),H(i+1,i));
H(i:i+1,i:n) = [c -s; s c]'*H(i:i+1,i:n);
H(:,i:i+1) = H(:,i:i+1)*[c -s; s c];
% 计算第i行到第n行的Householder矩阵
[v,beta] = house(H(i,i+1:n));
H(i:n,i+1:n) = H(i:n,i+1:n) - beta*(H(i:n,i+1:n)*v)*v';
H(i+1:n,i) = v;
H(:,i+1:n) = H(:,i+1:n) - (H(:,i+1:n)*v)*beta*v';
% 更新Q
Q(:,i:i+1) = Q(:,i:i+1)*[c -s; s c];
Q(:,i+1:n) = Q(:,i+1:n)-beta*(Q(:,i+1:n)*v)*v';
i = i+1;
end

% 双重隐式QR迭代
while n > 1
% 找到H的右下角2x2子阵
if n==2
mu = H(1,1);
a = H(1,2);
b = H(2,1);
c = H(2,2);
else
mu = H(n-1,n-1);
a = H(n-1,n);
b = 0;
c = H(n,n);
end
% 双重隐式QR迭代
[c,s] = givens(mu-a,c);
Qn = eye(n);
for k = 1:n-1
Qk = eye(n);
[c,s] = givens(H(k,k)-mu,H(k+1,k));
Qk(k:k+1,k:k+1) = [c -s; s c]';
H = Qk'*H*Qk;
Qn = Qn*Qk;
end
H(n-1:n,n-1:n) = [c -s; s c]'*H(n-1:n,n-1:n);
Q = Q*Qn;
n = n-1;
end

end

% Householder矩阵
function [v,beta] = house(x)
sigma = x(2:end)'*x(2:end);
v = [1; x(2:end)];
if sigma == 0 && x(1)>=0
beta = 0;
else
if x(1) <= 0
v(1) = -norm(x);
else
v(1) = norm(x);
end
beta = 2/(sigma+v(1)^2);
v = v/v(1);
end
end

% Givens旋转矩阵
function [c,s] = givens(a,b)
if b == 0
c = sign(a);
s = 0;
elseif abs(b) > abs(a)
t = -a/b;
s = 1/sqrt(1+t^2);
c = s*t;
else
t = -b/a;
c = 1/sqrt(1+t^2);
s = c*t;
end
end