U*A*U’ = B;其中B是拟上三角矩阵，拟上三角矩阵的定义是在矩阵的对角线上存在2x2大小的矩阵，而且矩阵U是正交矩阵，因为矩阵A的特征值和B的特征值相同。而且A的特征值出现在B的对角线上。计算特征值分解和SVD都依靠这个算法做最基本的处理，然后根据不同的任务有不同的处理。
计算schur分解的方法是是QR算法，这个算法的原理相当的简单，可以用如下的伪代码表示：
for i = 1 …
A(i-1)= QR
A(i) = R*Q
end
这段代码所做的变化类似于A(i) = R*Q = (Q’)*Q*R*Q = (Q’)*A(i-1)*Q;因此这段代码的基本思想就是使用正交矩阵Q不停的对矩阵A做相似变化。在这样的变化中将矩阵A的下半三角矩阵中的数全部消去。但是在实际中使用这样的算法是不现实的，因为每一次QR分解都需要大量的计算，同时完全的矩阵相乘R*Q也需要大量的计算。对这种方法的改进是首先将矩阵A化为Hessenberg型，然后对Hessenbert计算QR分解，对应的code如下:
function [H, U] = zhess(A)
%for any matrix A, turn it into a upper hessenberg matrix by orthogonal
%transformation
[m, n] = size(A);
if m ~= n
error('support square matrix only')
end
H = A;
U = eye(n);
for k=1:n-2
%compute the householder matrix
[v, beta] = zhouse(H(k+1:end, k));
temp_U = eye(n);
temp_U(k+1:n,k+1:n) = eye(n-k) - beta*v*(v');
H = temp_U*H;
U = U * temp_U;
%fprintf('after %d iteration\n', k);
%disp(H);
end
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足以解决我大部分的任务，慢慢来吧。
没有评论:
发表评论