QR分解
从矩阵分解的⾓度来看，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
基于householder vector计算QR分解的code如下：
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
这段Matlab代码为了确保数值计算的稳定性，做了⼀些优化。在计算第i个向量的时候，同时将i+1到n个向量在i个向量⽅向上的分量减去，对应的统计数据如下：
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和特征值的计算都有极其重要的基础作⽤。

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。