从矩阵分解的角度来看，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和特征值的计算都有极其重要的基础作用。