前面已经说过LU，Cholesky和QR分解，这次介绍的是实Schur分解。对这个分解的定义是任意一个矩阵A，可有如下形式的分解：
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
这段代码将矩阵A转换成上hessenberg矩阵。
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足以解决我大部分的任务，慢慢来吧。