双位移即Francis位移，双位移隐式QR迭代算法的思想是

1. 应用上Hessenberg算法将A矩阵转换为对应的上Hessenberg矩阵H；
2. 对H矩阵应用隐式QR迭代计算，直到H矩阵的对角线收敛到稳定值，即A对应的特征值。

第一步对应的matlab代码如下

% [H,Q] = hessred(A)
%
% Compute the Hessenberg decomposition H = Q'*A*Q using
% Householder transformations.
%
function [H,Q] = hessred(A)

  n = length(A);
  Q = eye(n);      % Orthogonal transform so far
  H = A;           % Transformed matrix so far

  for j = 1:n-2

    % -- Find W = I-2vv' to put zeros below H(j+1,j)
    u    = H(j+1:end,j);
    u(1) = u(1) + sign(u(1))*norm(u);
    v    = u/norm(u);

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

  end

end

第二步对应的matlab代码如下

% [H] = hessqr_francis(H)
%
% Compute a (double) implicit Hessenberg QR step with Francis shift.
% Compare to hessqr_basic.
%
function [H] = hessqr_francis(H)
  % Implicit QR step using a Francis double shift
  % (there should really be some re-scalings for floating point)

  % Compute double-shift poly and initial column of H^2 + b*H + c*I
  [b,c]   = francis_poly(H);
  C1      = H(1:3,1:2)*H(1:2,1);
  C1(1:2) = C1(1:2) + b*H(1:2,1);
  C1(1)   = C1(1)   + c;

  % Apply a similarity associated with the first step of QR on C
  % v      = C1;
  % v(1)   = v(1) + sign(v(1))*norm(v);
    
  v        = house(C1);
  H(1:3,:) = H(1:3,:)-2*v*(v'*H(1:3,:));
  H(:,1:3) = H(:,1:3)-(H(:,1:3)*(2*v))*v';

  % Do "bulge chasing" to return to Hessenberg form
  n = length(H);
  for j = 1:n-2
    k = min(j+3,n);

    % -- Find W = I-2vv' to put zeros below H(j+1,j), H := WHW'
    v          = house(H(j+1:k,j));
    H(j+1:k,:) = H(j+1:k,:)-2*v*(v'*H(j+1:k,:));
    H(:,j+1:k) = H(:,j+1:k)-(H(:,j+1:k)*(2*v))*v';
    H(k,j)     = 0;

  end

end

其中隐式Q定理中用到的Q，计算函数为francis_poly()

% [b,c] = francis_poly(H)
%
% Compute b, c s.t. z^2 + b*z + c = (z-sigma)(z-conj(sigma))
% where sigma is the Francis double shift for H.
%
function [b,c] = francis_poly(H)

  % Get shifts via trailing submatrix
  HH    = H(end-1:end,end-1:end);
  trHH  = HH(1,1)+HH(2,2);
  detHH = HH(1,1)*HH(2,2)-HH(1,2)*HH(2,1);

  if trHH^2 > 4*detHH   % Real eigenvalues

    % Use the one closer to H(n,n)
    lHH(1) = (trHH + sqrt(trHH^2-4*detHH))/2;
    lHH(2) = (trHH - sqrt(trHH^2-4*detHH))/2;
    if abs(lHH(1)-H(end,end)) < abs(lHH(2)-H(end,end))
      lHH(2) = lHH(1);
    else
      lHH(1) = lHH(2);
    end

    % z^2 + bz + c = (z-sigma_1)(z-sigma_2)
    b = -lHH(1)-lHH(2);
    c =  lHH(1)*lHH(2);

  else

    % In the complex case, we want the char poly for HH
    b = -trHH;
    c = detHH;

  end

想要了解代码的更多细节，请参考潘建瑜老师的《矩阵计算讲义》2016年4月19日版，4.5节 带位移的隐式QR迭代。

下面是笔者在琢磨过程中得出的两个性质：

‍

性质1
对于Householder矩阵H，有如下性质：HH=I; HH^(-1)=I; HH^T=I; 得出H=H^(-1)=H^T;

性质2
任何向量x对应的Householder矩阵的第一列与其平行。

hessqr_francis()中v        = house(C1);就是基于性质2得出的。