【矩阵分解】 QR， Householder变换

使用matlab，基于householder变化写了QR的实现过程

1、Householder变化算法

function [ H, v, beta ] = householder( x )
% x : inout param. x is a vector which size is n*1
% v and beta : is param which construct H matrix
% H is hoseholder Matrix. H = I - beta*v*v' 

%derive parameter v and beta
v = zeros(size(x));
beta = zeros(size(x));

%Houserholder Algorithm begin
x_len = length(x);
x_max = max(abs(x));
x = x./x_max;
zgama = x(2:end)'*x(2:end);
v(1) = 1;
v(2:end) = x(2:end);
if zgama == 0
    beta = 0;
else
    alpha = sqrt( x(1)^2 + zgama);
    if x(1) <=0
        v(1) = x(1) - alpha;
    else
        v(1) = -zgama./( x(1) + alpha);
    end
    beta = 2*v(1)^2./( zgama + v(1)^2 );
    v = v./v(1);
end
%beta = 2./(v'*v);
H = eye(x_len,x_len) - beta*v*v';
end

2、QR分解

A = rand(300,20);
tic
A_hang = size(A,1);
A_lie = size(A,2);
H = cell( A_lie, 1 );
Hs = eye(A_hang, A_hang);
R2 = A;
Q2 = eye(A_hang,A_hang);
for i = 1:A_lie
    [Hi, vi, betai] = householder(R2(i:end,i));
    %H{i} = blkdiag(eye(i-1), Hi);
    R2(i:end,i:end) = Hi*R2(i:end,i:end);
    Q2 = Q2*blkdiag(eye(i-1), Hi);
end
Q2(find(abs(Q2) < 1e-10)) = 0;
R2(find(abs(R2) < 1e-10)) = 0;
toc
[Q1, R1] = qr(A); %matlab内部算法（非常非常的快）
Q1*R1 - A
Q2*R2 - A
erroQ = Q1 + Q2
erroR = R1 + R2

主要写了QR的实现过程，算法很容易明白。但是针对你自己矩阵数据形式还要转换。

例如：稀疏矩阵就不要计算0元素的变换啦。