12.3. QR Decomposition

The QR decomposition is used to solve systems of linear equations (\(\m{A}\bm{x} = \m{Q\, R}\bm{x} = \bm{b}\)) and in the eigendecomposition algorithm described in Iterative QR Algorithm, The MATLAB command, [Q, R] = qr(A) finds the factors \(\bf{Q}\) and \(\bf{R}\). As shown in equation (12.4), \(\bf{Q}\) is unitary, and \(\bf{R}\) is upper triangular. The qrFactor function is a simple algorithm implementation. It has an economy option to remove rows of zeros from \(\bf{R}\) and unused columns from \(\bf{Q}\) for over-determined matrices. The qrFactor function makes use of the householder function to compute the Householder reflection matrices that, with multiplication, convert \(\bf{R}\) to upper triangular form. Then, the \(\bf{Q}\) matrix is built as an orthogonal product of transposed reflection matrices. The code is based on pseudo-code in Demmel’s Applied Numerical Linear Algebra text [DEMMEL97] and Golub and Van Loan’s Matrix Computations text [GOLUB13].

(12.4)\[\begin{split}\m{A} = \m{Q\, R} = \begin{bmatrix} \vertbar{} & \vertbar{} & {} & \vertbar{} \\ \bm{q}_1 & \bm{q}_2 & \cdots{} & \bm{q}_n \\ \vertbar{} & \vertbar{} & {} & \vertbar{} \end{bmatrix} \begin{bmatrix} r_{11} & r_{12} & \cdots{} & r_{1n} \\ 0 & r_{22} & \cdots{} & r_{2n} \\ \vdots{} & \vdots{} & \ddots{} & \vdots{} \\ 0 & 0 & \cdots{} & r_{nn} \end{bmatrix}\end{split}\]

function [Q, R] = qrFactor(A, econ)
% QRFACTOR - Simple QR factorization using Householder reflections.
% [Q, R] = QRFACTOR(A), yields Q and R factors such that A = Q*R.
%
% A - input real or complex m-by-n matrix, (m >= n).
% econ - (optional) enter 0 (zero) to get the economy QR,
%        where Q is m-by-n and R is n-by-n.
% Q - Unitary basis vectors of A (m-by-m)
% R - Upper triangular matrix (m-by-n)
%
% See also: QR, EIGQR, HOUSEHOLDER
%
    [m, n] = size(A); R = A; Q = eye(m);
    for k = 1:min((m-1), n)
        H = householder(R(k:m, k), m, k-1);
        R = H*R;
        R(k+1:m,k) = zeros(m-k,1); % no round-off error
        Q = Q*H';
        % For real numbers, H = H', but not in the complex case
    end
    % The econ option is useful for m > n matrices.
    if exist('econ','var') && econ == 0 && m > n
        Q = Q(:,1:n); R = R(1:n, :);
    end
end