前面已經說過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
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足以解決我大部分的任務，慢慢來吧。