%% Householder Transformation to Compute QR Decomposition
% This example shows how to use Householder transformations to compute the QR decomposition of a matrix A = QR, where Q is an orthogonal matrix and R is an upper triangular matrix.
% First, set the random number generator to the default value, and create a 6-by-3 rectangular matrix of random numbers from the standard normal distribution.

rng('default')
A = randn(6,3)

% To create a Householder matrix, use the function [v,beta] = gallery('house',x). This function takes a column vector x, and returns v and \beta such that H = I - \beta v v^T is a Householder matrix. The Householder transformations are used to zero out all but the first element of vector x.
% Compute a Householder matrix P_1 and perform the transformation A_1 = P_1A. The matrix A_1 has only zeros below the diagonal in the first column.

[v1,beta1] = gallery('house',A(:,1));
P1 = eye(6) - beta1*(v1*v1');
A1 = P1*A

% Next, compute a Householder matrix P_2 such that A_2 = P_2A_1 has only zeros below the diagonal in the first and second columns.

[v2,beta2] = gallery('house',A1(2:end,2));
v2 = [0;v2];
P2 = eye(6) - beta2*(v2*v2');
A2 = P2*A1

% Finally, compute a Householder matrix P_3 such that A_3 = P_3A_2 has only zeros in the subdiagonal elements.

[v3,beta3] = gallery('house',A2(3:end,3));
v3 = [0;0;v3];
P3 = eye(6) - beta3*(v3*v3');
R = P3*A2

% The matrix R = P_3P_2P_1A is an upper triangular matrix. Since the Householder matrices are involutory (the matrices are equal to their own inverses), the QR decomposition of A becomes A = QR with Q = P_1P_2P_3.

Q = P1*P2*P3

% Compare this result with the computation using the qr function.

[Qa,Ra] = qr(A)

% Verify that A = QR, within machine precision.

norm(A - Q*R)

%
% Copyright 2019 The MathWorks, Inc.