隐式Hessenberg矩阵QR分解常见的是应用Givens变换来保证Hessenberg矩阵的不可约性。直到Hessenberg矩阵的对角线收敛到其对应特征值。

下面的matlab代码应用Householder旋转的隐式QR分解，属于直接求解Hessenberg矩阵的特征值，并没有保证Hessenberg矩阵的不可约性。

% [H] = hessqr_basic(H)
%
% Compute one basic (unshifted) implicit Hessenberg QR step via
% Householder transformations.
%
function H = hessqr_basic(H)

  n = length(H);
  V = zeros(2,n-1);

  % Compute the QR factorization
  for j = 1:n-1

    % -- Find W_j = I-2vv' to put zero into H(j+1,j)
    u      = H(j:j+1,j);
    u(1)   = u(1) + sign(u(1))*norm(u);
    v      = u/norm(u);
    V(:,j) = v;

    % -- H := W_j H
    H(j:j+1,:) = H(j:j+1,:)-2*v*(v'*H(j:j+1,:));

  end

  % Compute RQ
  for j = 1:n-1

    % -- H := WHW', Q := QW
    v = V(:,j);
    H(:,j:j+1) = H(:,j:j+1)-(H(:,j:j+1)*(2*v))*v';

  end

end

下面对代码详解

‍

u      = H(j:j+1,j);
u(1)   = u(1) + sign(u(1))*norm(u);
v      = u/norm(u);

这段代码是Householder旋转中求解Householder向量。

‍

H(j:j+1,:) = H(j:j+1,:)-2*v*(v'*H(j:j+1,:));

这段其实是在求QR分解中的R，注意代码中的H并非Householder旋转中的Householder矩阵。

举个例子，假设矩阵M是6*6的矩阵。Householder旋转进行QR分解的步骤应该是：

‍

M1=H1*M
M2=H2*M1
M3=H3*M2
M4=H4*M3
M5=H5*M4

这样M5就是一个上三角阵。H1~H5表示对应的Householder矩阵。若是QR分解，就有

R=M5=H5*H4*H3*H2*H1*M
Q=H1*H2*H3*H4*H5

所以代码中第一个for循环就是在求解QR分解，最终代码中的H阵就是R。

QR迭代的第一步做QR分解已经完成，那么第二步求R*Q就是第二个for循环完成的事。对应M的QR分解示例。

R*Q=R*H1*H2*H3*H4*H5

确实对应于下面代码

v = V(:,j);
H(:,j:j+1) = H(:,j:j+1)-(H(:,j:j+1)*(2*v))*v';

整个过程处理得很有技巧。