QR Decomposition

The QR decomposition (or factorization) is an algorithm that converts a given matrix ${\bf A}$ into a product ${\bf A}={\bf Q}{\bf R}$ of an orthogonal matrix ${\bf Q}$ and a right or upper triangular matrix ${\bf R}$ with $r_{ij}=0\;(i>j)$. In the following we consider two methods for the QR decomposition.

• The Householder transformation:

  We first construct a Householder matrix ${\bf Q}_1$ based on the first column vector ${\bf a}_1=[a_{11},\;a_{21},\;\cdots,a_{N1}]^T$ of ${\bf A}=[{\bf a}_1,\cdots,{\bf a}_N]$, by which ${\bf a}_1$ will be converted into a vector with its last $N-1$ elements equal to zero:

  $${\bf Q}_1{\bf A}=\left[\begin{array}{c\vert ccc}a'_{11}&a'_{12}&\cdots&a'_{1N}\\ \hline 0& & & \\ \vdots & &{\bf A}'& \\ 0& & & \end{array}\right]$$ (64)

  We then construct second transformation matrix ${\bf Q}_2$

  $${\bf Q}_2=\left[\begin{array}{c\vert ccc} 1&0&\cdots&0\\ \hline 0& & & \\ \vdots & &{\bf Q}'_2 & \\ 0& & & \end{array}\right]$$ (65)

  where ${\bf Q}'_2$ is an $N-1$ dimensional Householder matrix to eliminate the last $N-2$ elements of the first column of ${\bf A}'$:

  $${\bf Q}_2{\bf A}'={\bf Q}_2{\bf Q}_1{\bf A}= \left[\begin{array}{... ...\hline 0&0& & & \\ \vdots & \vdots & &{\bf A}''& \\ 0&0& & & \end{array}\right]$$ (66)

  This process is repeated with the general kth transformation matrix

  $${\bf Q}_k=\left[\begin{array}{cc}{\bf I}_{k-1}&{\bf0}\\ {\bf0}&{\bf Q}'_k\end{array}\right]\;\;\;\;\;\;(k<N)$$ (67)

  where ${\bf Q}'_k$ is an $N-k+1$ dimensional Householder matrix to eliminate the last $N-k$ elements of the first column of the sub-matrix ${\bf A}^{(k-1)}$ of the same dimension obtained from the previous steps. After $N-1$ such iterations ${\bf A}$ becomes an upper triangular matrix:

  $${\bf Q}_{N-1}\cdots{\bf Q}_1{\bf A}={\bf R}$$ (68)

  Define ${\bf Q}=({\bf Q}_{N-1}\cdots{\bf Q}_1)^{-1}$ and pre-multiply it on both sides we get the QR decomposition of ${\bf A}$:

  $${\bf A}=({\bf Q}_{N-1}\cdots{\bf Q}_1)^{-1}{\bf R}={\bf Q}{\bf R}$$ (69)

  where ${\bf Q}$ defined aobve is an orthogonal matrix:
  

  $${\bf Q} = ({\bf Q}_{N-1}\cdots{\bf Q}_1)^{-1} ={\bf Q}_1^{-1}\cdots{\bf Q}_{N-1}^{-1}={\bf Q}_1^T\cdots{\bf Q}_{N-1}^T$$

  $$= ({\bf Q}_{N-1}\cdots{\bf Q}_1)^T=\left({\bf Q}^{-1}\right)^T$$ (70)

  
  

  To find the complexity of this algorithm, we note that the Householder transformation of a column vector is $O(N)$, which is needed for all columns of a submatrix in each of the iterations. Therefore the total complexity is $O(N^3)$. However, if ${\bf A}$ is a Hessenberg matrix (or more specially a tridiagonal matrix), then the Householder transformation of each column vector is constant $O(1)$ (for the only subdiagonal component), and the complexity becomes $O(N^2)$.

  Note that this method cannot convert ${\bf A}$ into a diagonal matrix, by performing Householder transformations on the rows (from the right) as well as on the columns (from the left), similar to the tridiagonalization of a symmetric matrix previously considered. This is because while post-multiplying a Householder matrix to convert a row, the column already converted by pre-multiplying a Householder matrix will be changed again.

• The Gram-Schmidt process:

  The QR decomposition of ${\bf A}$ can also be obtained by converting the $N$ column vectors in ${\bf A}=[{\bf a}_1,\cdots,{\bf a}_N]$, assumed to be independent, into a set of $N$ orthonormal vectors $\{{\bf q}_1,\cdots,{\bf q}_N\}$, which form the columns of the orthogonal matrix ${\bf Q}=[{\bf q}_1,\cdots,{\bf q}_N]$.

  We first find ${\bf q}_1={\bf a}_1/\vert\vert{\bf a}_1\vert\vert$ as the normalized version of the first column of ${\bf A}$, and then find each of the remaining orthogonal vector ${\bf q}_j$ ( $j=2,\cdots,N$), with all previous ${\bf q}_1,\cdots,{\bf q}_{j-1}$ already normalized:

  $${\bf q}'_j={\bf a}_j-\sum_{i=1}^{j-1}{\bf p}_{{\bf q}_i}({\bf a}_... ...{\bf q}_i ={\bf a}_j-\sum_{i=1}^{j-1}\left({\bf a}_j^T{\bf q}_i\right){\bf q}_i$$ (71)

  followed by normalization ${\bf q}_j={\bf q}'_j/\vert\vert{\bf q}'_j\vert\vert$ so that ${\bf q}_j^T{\bf q}_j=1$. Now ${\bf a}_j$ can be represented as a linear combination of $j$ vectors $\{{\bf q}_1,\cdots,{\bf q}_j\}$ as:

  $${\bf a}_j={\bf q}'_j+\sum_{i=1}^{j-1}({\bf a}_j^T{\bf q}_i){\bf q... ...rt{\bf q}_j'\vert\vert{\bf q}_j+\sum_{i=1}^{j-1}({\bf a}_j^T{\bf q}_i){\bf q}_i$$ (72)

  We define $r_{jj}=\vert\vert{\bf q}'_j\vert\vert$ and $r_{ij}={\bf a}_j^T{\bf q}_i\; (i=1,\cdots,j-1)$, and rewrite the above as

  $${\bf a}_j=r_{jj}{\bf q}_j+\sum_{i=1}^{j-1}r_{ij}{\bf q}_i =[{\bf ... ...j}\\ 0\\ \vdots\\ 0\end{array}\right] ={\bf Q}{\bf r}_j\;\;\;\;\;(j=1,\cdots,N)$$ (73)

  where ${\bf r}_j=[r_{1j},\cdots,r_{(j-1)j},r_{jj},0,\cdots,0]^T$ is the ith column of ${\bf R}$, and ${\bf q}_{j+1},\cdots,{\bf q}_N$ are yet to be generated (but they are multiplied by zero and play no role here). Combining the $N$ equations above in matrix form, we get

  $${\bf A}=[{\bf a}_1,\cdots,{\bf a}_N]={\bf Q}[{\bf r}_1,\cdots,{\bf r}_N] ={\bf Q}{\bf R}$$ (74)

  where ${\bf Q}$ is orthogonal and ${\bf R}$ is upper triangular.

  The complexity for computing ${\bf q}_j$ and ${\bf r}_j$ is $O(N^2)$, and the total complexity of this algorithm is $O(N^3)$.

The QR decomposition finds many applications, such as solving a linear system ${\bf Ax}={\bf b}$:

$${\bf A}{\bf x}={\bf Q}{\bf R}{\bf x}={\bf Q}{\bf y}={\bf b}$$ (75)

where ${\bf y}={\bf R}{\bf x}$ can be obtained as:

$${\bf y}={\bf Q}^{-1}{\bf b}={\bf Q}^T{\bf b}$$ (76)

and then ${\bf x}$ can be obtained by solving ${\bf R}{\bf x}={\bf y}$ with an upper triangular coefficient matrix ${\bf R}$ by back-substitution.

More importantly, the QR decomposition is the essential part of the QR algorithm for solving the eigenvalue problem of a general matrix, to be considered in the following section.

The Matlab code listed below carries out the QR decomposition by both the Householder transformation and the Gram-Schmidt method:

function [Q A]=QR_Householder(A)  % QR decomposition by Householder transformation
    [m n]=size(A);
    Q=eye(m);
    for i=1:n-1                   % for all n columns of A
        e=zeros(n+1-i,1);   
        e(1)=1;     
        a=A(i:n,i);               % first column of lower-right block of A
        u=a-norm(a)*e;
        v=u/norm(u);
        P=eye(n);
        P(i:n,i:n)=P(i:n,i:n)-2*v*v';   % Householder matrix
        A=P*A;
        Q=Q*P;
    end
end


function [Q R]=QR_GramSchmidt(A)  % QR decomposition by Gram-Schmidt method
    [m n]=size(A);
    Q=zeros(m,n);
    R=zeros(n);
    R(1,1)=norm(A(:,1));
    Q(:,1)=A(:,1)/R(1,1);
    for j=2:n
        Q(:,j)=A(:,j);
        for i=1:j-1
            Q(:,j)=Q(:,j)-A(:,j)'*Q(:,i)*Q(:,i);
            R(i,j)=A(:,j)'*Q(:,i);
        end
        R(j,j)=norm(Q(:,j));
        Q(:,j)=Q(:,j)/norm(Q(:,j));
    end    
end