Using MATLAB, the implementation process of QR was written based on householder change
 
1, Householder change algorithm
 
function [H, V, beta] = Householder (x)
% x:inout param. x is a vector which the size is n*1
% v and beta:is PA Ram 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 (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 decomposition
 
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 (ABS (Q2) < 1e-10) = 0;
R2 (ABS (R2) < 1e-10) = 0;
TOC
[Q1, R1] = QR (A);%matlab internal algorithm (very very fast)
q1*r1-a
q2*r2-a erroq
= Q1 + Q2
ErroR = R1 + R2
 
Mainly write the implementation of the QR process, the algorithm is easy to understand. But for yourself matrix data forms are also converted.
 
For example: sparse matrices do not compute the transformation of 0 elements.