【矩陣分解】 QR, Householder變換
阿新 • • 發佈:2018-11-12
使用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元素的變換啦。