QR分解

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