隐式qr方法matlab

### 隐式QR方法的MATLAB实现

隐式QR方法是一种用于计算矩阵特征值的有效数值算法。该方法的核心在于通过对给定矩阵施加一系列变换，最终将其转换成Schur形式，从而简化特征值求解过程。

#### Householder变换的应用
为了使矩阵逐步接近上三角结构，在隐式QR过程中会频繁使用Householder变换来减少非零元的数量。这种变换能够有效地将一般方阵转化为Hessenberg型矩阵——一种除了主对角线下一层外其余位置均为零的形式[^2]。

function H = householder(A)
    n = length(A);
    for k=1:n-2,
        x=A(k+1:n,k); % Select column to eliminate below diagonal.
        v=zeros(n-k,1);
        if norm(x)==0; continue; end;
        
        v(1)=sign(x(1))*norm(x)+x(1);
        v=v/norm(v)*sqrt(2*norm(x)^2+v'*A(k+1:n,k));
        A(k+1:n,k:n) = A(k+1:n,k:n)-v*(v'*A(k+1:n,k:n))/(v'*v);
        A(:,k+1:n) = A(:,k+1:n)-(A(:,k+1:n)*v)*(v'/2/(v'*v));    
    end
end

#### 上Hessenberg变换的作用
一旦获得了Hessenberg矩阵，则可以通过进一步的操作使其更易于处理。特别是当目标是获得完整的谱信息时，继续执行QR迭代直到达到所需的精度水平是非常重要的。此阶段的目标是在保持相似性的前提下尽可能多地消除次对角线上的条目[^3]。

function [Q,R]=qrstep(H,i,j)
% Perform a single QR step on rows/columns i through j of matrix H.

if (j-i)<1
    Q=eye(j-i+1); R=H(i:j,i:j);
else
    beta=-sign(H(i,i))*(abs(H(i,i))+...
         sqrt(sum(abs(diag(H(i:j,i:i)).^2))));
         
    u=[beta-H(i,i), zeros(1,j-i)]';
    u=u/sqrt(u'*u);

    Q=(eye(j-i+1)-2*u*u');
    R=triu(Q'*H(i:j,i:j));

    % Update the rest of the columns using Givens rotations.
    for p=j+1:size(H,2),
       c=R(j,p)/sqrt(R(j,j)^2+R(j,p)^2);
       s=-R(j,j)/sqrt(R(j,j)^2+R(j,p)^2);
       
       temp=c*R(:,p)+s*R(:,j);
       R(:,j)=c*R(:,j)-s*R(:,p);
       R(:,p)=temp;

       temp=c*H(:,p)+s*H(:,j);
       H(:,j)=c*H(:,j)-s*H(:,p);
       H(:,p)=temp;
    end
    
    % Apply transformation back onto original matrix.
    H(i:j,:)=Q'*H(i:j,:);
end

#### 完整流程概述
整个隐式QR方法涉及多个步骤：首先是利用Householder反射器构建初始的Hessenberg矩阵；接着反复施行带位移的QR分解并累积所有正交因子形成总变换矩阵\(Q\)；随着迭代次数增加，原输入矩阵逐渐趋近于其标准Schur形态，即准上三角状态下的酉相似版本\[T\]。

对于实际编程而言，上述操作通常封装在一个循环内完成，直至收敛条件被满足为止：

function T=eig_qr(A,tol,maxiter)
    m=size(A,1);
    I=speye(m,m);
    shift=trace(A)/m*I;
    B=A-shift;
    [~,B]=householder(B);
    
    iter=0;
    while true && iter<maxiter
        [Q,B]=qrstep(B,1:m-1);
        B=B-diag(diag(B).*I)+diag(shift.*ones(size(diag(B)))) ;
        err=max(max(abs(tril(B,-1))));
        disp(['Iteration ',num2str(iter),' error=', num2str(err)]);
        if err<=tol || isnan(err)|| isinf(err), break ; end
        iter=iter+1;
    end
    T=B+shift;
end