1. 程式人生 > >【矩陣分解】 QR, Householder變換

【矩陣分解】 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元素的變換啦。