雑記帳

function [Q,R] = QR_Householder(A)
%===============================================
% ハウスホルダー変換によるQR分解
%===============================================
[n,m]=size(A);

for k=0:m-1
    if k==0
        R=A;
        Q=eye(n);
    end
    C2=R(k+1:n,k+1:m);
    x=C2(:,1);
    y=zeros(length(x),1);
    y(1)=-sign(x(1))*sqrt(x'*x);
    u=(x-y)/sqrt((x-y)'*(x-y));
    if k==0
        H=eye(length(x))-2*u*u';
    else
        H=zeros(n);
        H(1:k,1:k)=eye(k);
        H(k+1:n,k+1:n)=eye(length(x))-2*u*u';
    end
    tR=H*R;
    tQ=H*Q;
    clear R Q;
    R=tR;
    Q=tQ;
    clear tR tQ C2 x y u H;
end

end