实Schur分解

function [NH] = zhessqr(H)

% perform QR algorithm on upper hessenberg matrix

% firstly, we need to verity this is a hessberg matrix

[m, n] = size(H);

if m ~= n

    error(‘error, support square matrix only‘)

end

NH = H;

c = zeros(1, n-1);

s = zeros(1, n-1);

for k=1:n-1

    %compute gives rotation at first

    [c(k), s(k)] = zgivens(NH(k, k), NH(k+1, k));

    p = [c(k) s(k); -s(k) c(k)];

    NH(k:k+1, k:n) = (p‘)*NH(k:k+1, k:n);

    %fprintf(‘after %d iteration\n‘, k);

    %disp(NH);

end

for k=1:n-1

    p = [c(k) s(k); -s(k) c(k)];

    NH(1:k+1, k:k+1) = NH(1:k+1,k:k+1)*p;

end

这段代码计算Hessenberg型矩阵A的一个QR step，在这里使用givens旋转来获得矩阵的QR分解提高了效率。

在QR算法中最主要的步骤就是QR step，首先做QR分解，然后R*Q，为了加快算法收敛的速度，使用了基于位移的QR算法，基本的伪代码如下：

for i=1 ...

    A - a*I = QR

    R*Q + a*I = A

end

使用这个方法可以加快QR算法的收敛速度。在这个算法的基础上有隐式双位移算法，称为Francis QR。首先给出显式的双位移算法：

for i=1:infinite

    H(0) - a1*I = Q(1)*R(1)

    H(1) = R(1)*Q(1) + a1*I

    H(1) - a2*I = Q(2)*R(2)

    H(2) = R(2)*Q(2) + a2*I

end

略加推倒可以获得如下的公式：

>> H(2) = (Z’)*H*Z; M = Z*R; M = (H-a1*I)(H-a2*I);

Francis QR可以在不显式的构造矩阵M的情况下，完成H2=Z’ * H * Z;

Francis QR算法可以使用依赖于隐式Q定理，对应的Matlab代码如下：

function [H, U] = zfrancisqr(A)

%compute one of the step by implicitly shifted QR step

[m, n] = size(A);

if m ~= n

    error(‘support square matrix only‘)

end

m = n-1;

s = A(m, m) + A(n, n);

t = A(m, m)*A(n, n) - A(m, n)*A(n, m);

x = A(1,1)*A(1,1) + A(1,2)*A(2,1) - s*A(1,1) + t;

y = A(2,1)*(A(1,1) + A(2,2) - s);

z = A(2,1)*A(3,2);

for k=0:n-3

    [v, beta] = zhouse([x y z]‘);

    q = max([1 k]);

    %orthogonal transformation

    ot = (eye(3) - beta*v*(v‘));

    A(k+1:k+3,q:n) = ot*A(k+1:k+3, q:n);

    r = min([k+4 n]);

    A(1:r, k+1:k+3) = A(1:r, k+1:k+3)*(ot‘);

    x = A(k+2, k+1);

    y = A(k+3, k+1);

    if k < n-3

        z = A(k+4, k+1);

    end

end

[v, beta] = zhouse([x y]‘);

ot = eye(2) - beta*v*(v‘);

A(n-1:n, n-2:n) = ot*A(n-1:n, n-2:n);

A(1:n, n-1:n) = A(1:n, n-1:n)*(ot‘);

H = A;

U = eye(n);

隐式Q定理的基本内容如下：对于矩阵A，存在两个不同的相似变换Q’*A*Q = H, V’*A*V=G，H和G是上Hessenberg矩阵，如果Q和V的第一列相同，那么这两个不同的相似变换就是等价的。因此Francis QR的第一步就是计算矩阵M的第一列，然后使用householder reflector将之变成e1，然后将变换后的矩阵转换成上Hessenberg矩阵，这个时候就完成了一步Francis QR。这个方式之所以可以使用隐式Q定理是因为第一个householder reflector是针对M的第一列计算的，而且后来的householder reflector的第一列都是e1，因为最后计算出的变换矩阵的第一列和直接在M上计算QR分解是相同的。

这就是QR算法涉及的主要内容，事实上QR算法的研究很多，但是这里这是给出基本的，而且没有给出完成的计算程序，是因为我现在还不能完全理解整个过程。下面对于Spectral Decomposition和Singular Value Decompositon介绍也要搁置一段时间，第一是因为两个算法很复杂，需要一段时间来理解；第二个原因是因为现在没有很强的需求去研究到这样的细节。目前依靠LAPACK和Matlab足以解决我大部分的任务，慢慢来吧。