從矩陣分解的角度來看，LU和Cholesky分解目標在於將矩陣轉化為三角矩陣的乘積,所以在LAPACK種對應的名稱是trf（Triangular Factorization）。QR分解的目的在於將矩陣轉化成正交矩陣和上三角矩陣的乘積，對應的分解公式是A=Q*R。正交矩陣有很多良好的性質，比如矩陣的逆和矩陣的轉置相同，任意一個向量和正交矩陣的乘積不改變向量的2范數等等。QR分解可以用於求解線性方程組，線性擬合。更重要的是QR分解是QR算法的基礎，可以用於各種特征值問題，所以QR分集的應用非常廣泛。與QR分解對應的還有一個LQ分解，LQ分解其實就是transpose(A)的QR分解。下面分別給出QR分解的詳細計算方法。
QR分解的計算分解有兩種不同的思路：第一種是使用正交矩陣Q左乘A，選擇特殊的矩陣Q可以將A種的某些元素消零；第二種方法是逐漸從A矩陣的列中構造出正交矩陣，並且計算出對應的R。首先介紹第一種方法，暫且稱之為消零法吧。
第一種消零法基於householder reflector。定義P=(I - 2*v*transpose(v) ./ (transpose(v)*v))，如果P*x = epison*I(:,1)，那么v稱為householder vector。可以發現找到向量x對應的householder vector可以將除第一個元素以為的n-1個元法全部消零，而且P本身是正交矩陣。為了存儲方便將v(1)z標准化為1存儲，計算code如下：
function[v, beta] =zhouse(x)
%generate householder reflector vector v, so that P=(I - 2*(v*v')/(v'*v)
% Px = epison*eye(:,1), good computation methods
[m, n] = size(x);
if n ~= 1
error('only calculate householder reflector for vector not matrix')
end
delta = (x(2:end)')*x(2:end);
v = zeros(m, 1);
v(1) = 1;
v(2:end) = x(2:end);
if abs(delta) < 1e-8
beta = 0;
else
mu = sqrt( x(1)^2 + delta );
if x(1) <= 0
v(1) = x(1) - mu;
else
v(1) = -1*(delta ./ (x(1) + mu));
end
beta = 2 * (v(1)^2) ./ (delta + v(1)^2);
v = v ./ v(1);
end
function[Q, R] =zhouseqr(A)
%compute the QR decomposition of square matrix A
% A = Q*R; Q is orthogonal matrix and R is upper trigular matrix
[m, n] = size(A);
if m ~= n
error('support squre matrix only')
end
temp = A;
Q = eye(n);
for k=1:n-1
[v, beta] = zhouse(temp(k:end, k));
reflector = eye(m-k+1) - beta*v*(v');
temp(k:m, k:m) = reflector*temp(k:m, k:m);
%construct orthogonal matrix
t = eye(n);
t(k:m, k:m) = reflector;
Q = Q*(t');
%
end
R = temp;
這個計算方法非常直接，對矩陣A的每一列計算出對應的householder reflector。下面給出Matlab代碼和Matlab自帶的qr分解之間的差異，考察包括Q的正交性如何，Q*R乘積與原始矩陣的差異，給出對應的均值和方差。
mean of my norm(q*t(q) - eye(n) : 1.65376e-15
mean of matlab norm(q*t(q) - eye(n) : 1.47759e-15
variance of my norm(q*t(q) - eye(n) : 8.32936e-32
variance of matlab norm(q*t(q) - eye(n) : 9.17174e-32
mean of my norm(q*r - test) : 4.42572e-15
mean of matlab norm(q*r - test) : 3.75022e-15
variance of my norm(q*r - test) : 6.34439e-31
variance of matlab norm(q*r - test) : 7.59865e-31
對這個結果進行分析會發現，自己寫的QR分解的均值總比Matlab自帶的差，但是在一個數量級之內。從方差的角度來分析，可以發現自己寫的QR分解的方差會更好寫。我覺得一種可能的解釋是，自己寫的QR分解穩定的比Matlab自帶的差？但是不管怎么說分解的穩定性比較好就對了。
第二種消零法基於givens rotation。givens rotation是特定的對矩陣中特定的某個元素消零。相比於householder reflector用於在矩陣中引入大量的零元素。計算的基本原理可以從一個2X2的矩陣P，左乘2x1的向量x，並消去第二個元素知道。givens rotation對應的Matlab代碼如下：
function [c, s] = zgivens(a, b)
%compute givens rotation for vector [a b]'
if abs(b) < 1e-8
c = 1; s = 0;
else
if abs(b) > abs(a)
tau = -1*a ./ b ;
s = 1 ./ sqrt(1 + tau.^2);
c = s * tau;
else
tau = -1*b ./ a;
c = 1 ./ sqrt(1 + tau.^2);
s = c * tau;
end
end
使用givens rotation來完成QR分解也非常的直接。對於矩陣的第i列，將i+1行到n行的元素全部消零。由於givens rotation需要的矩陣P是正交矩陣，因為最終的乘積仍然是正交矩陣，完全滿足QR分解的要求，對應的Matlab如下：
function[Q, R] =zgivensqr(A)
%compute QR decomposition based on givens rotation
[m, n] = size(A);
if m ~= n
error('support square matrix')
end
Q = eye(n);
temp = A;
for k=1:n-1
for j = m:-1:k+1
[c, s] = zgivens(temp(j-1,k), temp(j,k));
%update the original matrix
for l=k:n
tau1 = temp(j-1, l);
tau2 = temp(j, l);
temp(j-1, l) = c*tau1 - s*tau2;
temp(j, l) = s*tau1 + c*tau2;
end
%update the orthogonal matrix
for l=1:n
tau1 = Q(j-1, l);
tau2 = Q(j, l);
Q(j-1, l) = c*tau1 - s*tau2;
Q(j, l) = s*tau1 + c*tau2;
end
end
end
R = temp;
將基於givens rotation的QR分解與Matlab自帶的Matlab分解的比較統計數據如下：
Comparison of orthogonality of matrix Q
mean of my norm(q*t(q) - eye(n) : 1.66839e-15
mean of matlab norm(q*t(q) - eye(n) : 1.4748e-15
variance of my norm(q*t(q) - eye(n) : 8.53389e-32
variance of matlab norm(q*t(q) - eye(n) : 9.68473e-32
Comparison of QR decomposition
mean of my norm(q*r - test) : 13.3161
mean of matlab norm(q*r - test) : 3.73495e-15
variance of my norm(q*r - test) : 2.65531
variance of matlab norm(q*r - test) : 7.38973e-31
還是從兩個不同的方法對QR分解做比較，一個是分解獲得的Q的精度和穩定性。從結果可以看出，Q的精度比Matlab的稍差，但是Q的穩定性比Matlab的稍好。無論如何，這樣的結果與基於householder reflector的QR分解相同。但是涉及到QR分解本身，發現分解的穩定性與精度比Matlab的差了很多。雖然對這段代碼做過多次檢查，但是沒有發現出對應的bug，可能最終的計算結果就是這樣吧。
第二種從需要分解的矩陣A中構造出正交矩陣Q。這種計算的方法方法是基於如下的觀察，如果有兩個向量x和y，如果將y投影到x方向的分量減去，那么y剩下的部分必然與x正交。如果有n個向量，從i個向量開始，分別減去第i個向量在1到i-1個向量方向上的投影，那么第i個向量剩下的部分必然與1到i-1個向量正交。將這種算法寫成code如下：
function[Q, R] =zgsqr(A)
%compute QR decomposition based on intuitive way
[m, n] = size(A);
if m ~= n
error('support square matrix only')
end
Q = zeros(n, n);
R = zeros(n, n);
for k = 1:n
R(k, k) = norm(A(:,k), 2);
Q(:, k) = A(:,k) ./ R(k, k);
for j=k+1:n
R(k, j) = Q(:,k)'*A(:,j);
A(:, j) = A(:, j) - Q(:,k)*R(k,j);
end
end
mean of my norm(q*t(q) - eye(n) : 1.64192e-14
mean of matlab norm(q*t(q) - eye(n) : 1.47988e-15
variance of my norm(q*t(q) - eye(n) : 1.02173e-26
variance of matlab norm(q*t(q) - eye(n) : 9.49148e-32
Comparison of QR decomposition
mean of my norm(q*r - test) : 1.12525e-15
mean of matlab norm(q*r - test) : 3.74833e-15
variance of my norm(q*r - test) : 2.92379e-32
variance of matlab norm(q*r - test) : 7.40202e-31
從統計數據來看，分解獲得的Q的精度和穩定性都比Matlab自帶的qr分解要查，但是可以滿足日常的需求。同時QR分解的精度和穩定性卻比Matlab自帶的qr分解更好，不得不說這是讓人感到很驚喜的一個方面。如果不使用類似的技術來改進數值穩定性，那么對應的統計數據如下：
Comparion of orthogonality of matrix Q
mean of my norm(q*t(q) - eye(n) : 1.4272e-14
mean of matlab norm(q*t(q) - eye(n) : 1.48239e-15
variance of my norm(q*t(q) - eye(n) : 5.85718e-27
variance of matlab norm(q*t(q) - eye(n) : 9.28478e-32
Comparison of QR decomposition
mean of my norm(q*r - test) : 1.1205e-15
mean of matlab norm(q*r - test) : 3.76528e-15
variance of my norm(q*r - test) : 2.77064e-32
variance of matlab norm(q*r - test) : 7.48176e-31
貌似在精度和穩定性上沒有太多的差別，可能是在特殊的情況下才表現出來的吧。暫時數據上不支持教科書上的結論，希望有能人志士給出解答。使用如下的code進行計算的：
function[Q, R] =zrawgsqr(A)
%compute the QR decomposition based on raw gram-schmit approach
[m, n] = size(A);
if m ~= n
error('support square matrix only');
end
Q = zeros(n, n);
R = zeros(n, n);
for k=1:n
for j=1:k-1
R(j, k) = Q(:,j)' * A(:, k);
A(:, k) = A(:, k) - R(j, k)*Q(:,j);
end
R(k, k) = norm(A(:,k), 2);
Q(:, k) = A(:, k) ./ R(k, k);
end
關於QR分解的具體實現就在這里給出了，總體來說各個方法有好處有壞處，關鍵看你需要什么類型的功能。下一次將給出schur分解的詳細實現，schur分解對於SVD和特征值的計算都有極其重要的基礎作用。