6.4. QR decomposition

QR decomposition is another decomposition method that factorises a matrix in the product of an orthogonal matrix \(Q\) and an upper triangular matrix \(R\) such that

\[ \begin{align*} A = QR. \end{align*} \]

QR decomposition is particularly useful for solving least-squares problems and computing eigenvalues. Unlike LU decomposition, QR decomposition can be applied to rectangular matrices and is numerically more stable.

Definition 6.3 (Orthogonal vectors)

A set of vectors \(\{ \mathbf{v}_1 ,\mathbf{v}_2 ,\mathbf{v}_3 ,\dots \}\) is said to be orthogonal if \(\mathbf{v}_i \cdot \mathbf{v}_j =0\) for \(i\not= j\). Furthermore the set is said to be orthonormal if \(\mathbf{v}_i\) are unit vectors.

Definition 6.4 (Orthogonal matrix)

An orthogonal matrix is a matrix whose columns form an orthonormal set. Equivalently

\[ \begin{align*} A^\mathsf{T} A = AA^\mathsf{T} = I. \end{align*} \]

For example consider the matrix

\[\begin{split} \begin{align*} A= \begin{pmatrix} 0.8 & -0.6\\ 0.6 & 0.8 \end{pmatrix}. \end{align*} \end{split}\]

This is an orthogonal matrix since

\[\begin{split} \begin{align*} A^\mathsf{T} A=\begin{pmatrix} 0.8 & 0.6\\ -0.6 & 0.8 \end{pmatrix} \begin{pmatrix} 0.8 & -0.6\\ 0.6 & 0.8 \end{pmatrix} = \begin{pmatrix} 1 & 0\\ 0 & 1 \end{pmatrix}=I. \end{align*} \end{split}\]

---

6.4.1. QR decomposition using the Gram-Schmidt process

Consider the \(3 \times 3\) matrix \(A\) represented as the concatenation of the column vectors \(\mathbf{a}_1, \mathbf{a}_2, \mathbf{a}_3\)

\[\begin{split} \begin{pmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{pmatrix} = \begin{pmatrix} \uparrow & \uparrow & \uparrow \\ \mathbf{a}_1 & \mathbf{a}_2 & \mathbf{a}_3 \\ \downarrow & \downarrow & \downarrow \end{pmatrix} \end{split}\]

where \(\mathbf{a}_j = (a_{1j}, a_{2j}, a_{3j})^\mathsf{T}\). To calculate the orthogonal \(Q\) matrix we need a set of vectors \(\{\mathbf{q}_1, \mathbf{q}_2, \mathbf{q}_3\}\) that span the same space as \(\{\mathbf{a}_1, \mathbf{a}_2, \mathbf{a}_3\}\) but are orthogonal to each other. One way to do this is using the Gram-Schmidt process. Consider the diagram shown in Fig. 6.1 that shows the first two non-orthogonal vectors \(\mathbf{a}_1\) and \(\mathbf{a}_2\).

Fig. 6.1 The Gram-Schmidt process

We begin by letting \(\mathbf{u}_1 = \mathbf{a}_1\) and seek to find a vector \(\mathbf{u}_2\) that is orthogonal to \(\mathbf{u}_1\) and in the same span as \(\mathbf{a}_2\). To do this we subtract the vector projection of \(\mathbf{a}_2\) onto \(\mathbf{u}_1\), i.e.,

\[ \mathbf{u}_2 = \mathbf{a}_2 - \left( \frac{\mathbf{u}_1 \cdot \mathbf{a}_1}{\mathbf{u}_1 \cdot \mathbf{u}_1} \right) \mathbf{u}_1.\]

We can simplify things by normalising \(\mathbf{u}_1\) such that \(\mathbf{q}_1 = \dfrac{\mathbf{u}_1}{\|\mathbf{u}_1\|}\) so we have

(6.4)\[ \mathbf{u}_2 = \mathbf{a}_2 - (\mathbf{q}_1 \cdot \mathbf{a}_2) \mathbf{q}_1, \]

which we also normalise to give \(\mathbf{q}_2 = \dfrac{\mathbf{u}_2}{\|\mathbf{u}_2\|}\). For the next vector \(\mathbf{a}_3\) we want the vector \(\mathbf{u}_3\) to be orthogonal to both \(\mathbf{q}_1\) and \(\mathbf{q}_2\) so we subtract both the vector projections of \(\mathbf{a}_3\) onto \(\mathbf{q}_1\) and \(\mathbf{q}_2\)

(6.5)\[ \mathbf{u}_3 = \mathbf{a}_3 - (\mathbf{q}_1 \cdot \mathbf{a}_3) \mathbf{q}_1 - (\mathbf{q}_2 \cdot \mathbf{a}_3) \mathbf{q}_2, \]

which is normalised to give \(\mathbf{q}_3\). Using \(\mathbf{u}_i = (\mathbf{q}_i \cdot \mathbf{a}_i) \mathbf{q}_i\) and rearranging (6.4) and (6.5) to make \(\mathbf{a}_i\) the subject and

\[\begin{split} \begin{align*} \mathbf{a}_1 &= (\mathbf{q}_1 \cdot \mathbf{a}_1) \mathbf{q}_1, \\ \mathbf{a}_2 &= (\mathbf{q}_1 \cdot \mathbf{a}_2) \mathbf{q}_1 + (\mathbf{q}_2 \cdot \mathbf{a}_2) \mathbf{q}_2, \\ \mathbf{a}_3 &= (\mathbf{q}_1 \cdot \mathbf{a}_3) \mathbf{q}_1 + (\mathbf{q}_2 \cdot \mathbf{a}_3) \mathbf{q}_2 + (\mathbf{q}_3 \cdot \mathbf{a}_3) \mathbf{q}_3, \\ & \vdots \end{align*} \end{split}\]

The columns of the \(Q\) matrix are formed using the vectors \(\mathbf{q}_i\)

\[\begin{split} Q = \begin{pmatrix} \uparrow & \uparrow & & \uparrow \\ \mathbf{q}_1 & \mathbf{q}_2 & \cdots & \mathbf{q}_n \\ \downarrow & \downarrow & & \downarrow \end{pmatrix}, \end{split}\]

and the elements of the \(R\) matrix are formed from the dot products of the \(\mathbf{q}_i\) and \(\mathbf{a}_j\) vectors such that

\[\begin{split} R = \begin{pmatrix} (\mathbf{q}_1 \cdot \mathbf{a}_1) & (\mathbf{q}_1 \cdot \mathbf{a}_2) & \cdots & (\mathbf{q}_1 \cdot \mathbf{a}_n) \\ 0 & (\mathbf{q}_2 \cdot \mathbf{a}_2) & \cdots & (\mathbf{q}_2 \cdot \mathbf{a}_n) \\ \vdots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & (\mathbf{q}_n \cdot \mathbf{a}_n) \end{pmatrix}. \end{split}\]

6.4.1.1. Positive diagonal elements of \(R\)

It has become standard convention to impose a condition that the diagonal elements of \(R\) must be positive so that we get a unique factorisation. Consider \(A = QR\) then multiplying a column of \(Q\) by \(-1\) and the corresponding row of \(R\) by \(-1\) doesn’t change the product \(QR\). So if any of the diagonal elements of \(R\) are negative we multiply the corresponding column of \(Q\) and row of \(R\) by \(-1\).

Algorithm 6.6 (QR decomposition using the Gram-Schmidt process)

Inputs: An \(m \times n\) matrix \(A\).

Outputs: An \(m \times n\) orthogonal matrix \(Q\) and and \(n \times n\) upper triangular matrix \(R\) such that \(A = QR\).

• For \(j = 1, \ldots, n\) do

  • \(r_{ij} \gets \mathbf{q}_i \cdot \mathbf{a}_j, \qquad i = 1, \ldots, j - 1\)

  • \(\mathbf{u}_j \gets \mathbf{a}_j - \displaystyle\sum_{i=1}^{j-1} r_{ij} \mathbf{q}_i\)

  • \(r_{jj} \gets \| \mathbf{u}_j \|\)

  • \(\mathbf{q}_j \gets \dfrac{\mathbf{u}_j}{r_{jj}}\)

• \(Q \gets (\mathbf{q}_1, \mathbf{q}_2, \ldots, \mathbf{q}_n)\)

• For \(i = 1, \ldots, n\) do

  • If \(r_{ii} < 0\) then

    • \(R(i,: ) \gets -R(i,:)\)

    • \(Q(:,i) \gets -Q(:,i)\)

• Return \(Q\) and \(R\)

Example 6.8

Calculate the QR decomposition of the following matrix using the Gram-Schmidt process

\[\begin{split} \begin{align*} A = \begin{pmatrix} -1 & -1 & 1 \\ 1 & 3 & 3 \\ -1 & -1 & 5 \\ 1 & 3 & 7 \end{pmatrix}. \end{align*} \end{split}\]

---

Solution

\[\begin{split} \begin{align*} j &= 1: & \mathbf{u}_1 &= \mathbf{a}_1 = \begin{pmatrix} -1 \\ 1 \\ -1 \\ 1 \end{pmatrix}, \\ && r_{11} &= \| \mathbf{u}_1 \| = 2, \\ && \mathbf{q}_1 &= \frac{\mathbf{u}_1}{r_{11}} = \frac{1}{2} \begin{pmatrix} -1 \\ 1 \\ 1 \\ 1 \end{pmatrix} = \begin{pmatrix} -\frac{1}{2} \\ \frac{1}{2} \\ -\frac{1}{2} \\ \frac{1}{2} \end{pmatrix}. \\ j &= 2: & r_{12} &= \mathbf{q}_1 \cdot \mathbf{a}_2 = \begin{pmatrix} -\frac{1}{2} \\ \frac{1}{2} \\ -\frac{1}{2} \\ \frac{1}{2} \end{pmatrix} \cdot \begin{pmatrix} -1 \\ 3 \\ -1 \\ 3 \end{pmatrix} = 4, \\ && \mathbf{u}_2 &= \mathbf{a}_2 - r_{12} \mathbf{q}_1 = \begin{pmatrix} -1 \\ 3 \\ -1 \\ 3 \end{pmatrix} - 4 \begin{pmatrix} -\frac{1}{2} \\ \frac{1}{2} \\ -\frac{1}{2} \\ \frac{1}{2} \end{pmatrix} = \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \end{pmatrix}, \\ && r_{22} &= \| \mathbf{u}_2 \| = 2, \\ && \mathbf{q}_2 &= \frac{\mathbf{u}_2}{r_{22}} = \frac{1}{2} \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \end{pmatrix} = \begin{pmatrix} \frac{1}{2} \\ \frac{1}{2} \\ \frac{1}{2} \\ \frac{1}{2} \end{pmatrix}. \\ j &= 3: & r_{13} &= \mathbf{q}_1 \cdot \mathbf{a}_3 = \begin{pmatrix} -\frac{1}{2} \\ \frac{1}{2} \\ -\frac{1}{2} \\ \frac{1}{2} \end{pmatrix} \cdot \begin{pmatrix} 1 \\ 3\\ 5 \\ 7 \end{pmatrix} = 2, \\ && r_{23} &= \mathbf{q}_2 \cdot \mathbf{a}_3 = \begin{pmatrix} \frac{1}{2} \\ \frac{1}{2} \\ \frac{1}{2} \\ \frac{1}{2} \end{pmatrix} \cdot \begin{pmatrix} 1 \\ 3 \\ 5 \\ 7 \end{pmatrix} = 8, \\ && \mathbf{u}_3 &= \mathbf{a}_3 - r_{13} \mathbf{q}_1 - r_{23} \mathbf{q}_2 = \begin{pmatrix} 1 \\ 3 \\ 5 \\ 7 \end{pmatrix} - 2 \begin{pmatrix} -\frac{1}{2} \\ \frac{1}{2} \\ -\frac{1}{2} \\ \frac{1}{2} \end{pmatrix} - 8 \begin{pmatrix} \frac{1}{2} \\ \frac{1}{2} \\ \frac{1}{2} \\ \frac{1}{2} \end{pmatrix} = \begin{pmatrix} -2 \\ -2 \\ 2 \\ 2 \end{pmatrix}, \\ && r_{33} &= \| \mathbf{u}_3 \| = 4, \\ && \mathbf{q}_3 &= \frac{\mathbf{u}_3}{r_{33}} = \frac{1}{4} \begin{pmatrix} -2 \\ -2 \\ 2 \\ 2 \end{pmatrix} = \begin{pmatrix} -\frac{1}{2} \\ -\frac{1}{2} \\ \frac{1}{2} \\ \frac{1}{2} \end{pmatrix}. \end{align*} \end{split}\]

The diagonal elements of \(R\) are all positive so we don’t need to make any adjustments. Thus the \(Q\) and \(R\) matrices are given by

\[\begin{split} \begin{align*} Q &= \begin{pmatrix} -\frac{1}{2} & \frac{1}{2} & -\frac{1}{2} \\ \frac{1}{2} & \frac{1}{2} & -\frac{1}{2} \\ -\frac{1}{2} & \frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & \frac{1}{2} & \frac{1}{2} \end{pmatrix}, & R &= \begin{pmatrix} 2 & 4 & 2 \\ 0 & 2 & 8 \\ 0 & 0 & 4 \end{pmatrix}. \end{align*} \end{split}\]

We can check that this is correct by verifying that \(A = QR\) and that \(Q^\mathsf{T}Q = I\).

\[\begin{split} \begin{align*} \begin{pmatrix} -\frac{1}{2} & \frac{1}{2} & -\frac{1}{2} \\ \frac{1}{2} & \frac{1}{2} & -\frac{1}{2} \\ -\frac{1}{2} & \frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & \frac{1}{2} & \frac{1}{2} \end{pmatrix} \begin{pmatrix} 2 & 4 & 2 \\ 0 & 2 & 8 \\ 0 & 0 & 4 \end{pmatrix} &= \begin{pmatrix} -1 & -1 & 1 \\ 1 & 3 & 3 \\ -1 & -1 & 5 \\ 1 & 3 & 7 \end{pmatrix} = A, \\ \begin{pmatrix} -\frac{1}{2} & \frac{1}{2} & -\frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & \frac{1}{2} & \frac{1}{2} & \frac{1}{2} \\ -\frac{1}{2} & -\frac{1}{2} & \frac{1}{2} & \frac{1}{2} \end{pmatrix} \begin{pmatrix} -\frac{1}{2} & \frac{1}{2} & -\frac{1}{2} \\ \frac{1}{2} & \frac{1}{2} & -\frac{1}{2} \\ -\frac{1}{2} & \frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & \frac{1}{2} & \frac{1}{2} \end{pmatrix} &= \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} = I. \end{align*} \end{split}\]

6.4.2. QR decomposition using Householder reflections

Another method we can use to calculate the QR decomposition of a matrix is by the use of Householder reflections.

Definition 6.5 (Householder matrix)

The Householder matrix is an orthogonal reflection matrix of the form

(6.6)\[ H = I - 2 \frac{\mathbf{v} \mathbf{v}^\mathsf{T}}{\mathbf{v}^\mathsf{T} \mathbf{v}}, \]

where \(\mathbf{v}\) is a non-zero column vector. The matrix \(H\) is symmetric and orthogonal.

Fig. 6.2 The vector \(\mathbf{x}\) is reflected about the hyperplane normal to the vector \(\mathbf{v}\) so that it is parallel to the basis vector \(\mathbf{e}_1\).

Consider the diagram in Fig. 6.2 where the vector \(\mathbf{x}\) is reflected about the hyperplane so that it is parallel to the first basis vector of the standard basis \(\mathbf{e}_1 = (1, 0, \ldots, 0)^\mathsf{T}\). The Householder matrix that achieves this reflection is written in terms of the vector \(\mathbf{v}\) which is normal to the reflecting hyperplane

\[ \mathbf{v} = \mathbf{x} - \|\mathbf{x}\|\mathbf{e}_1. \]

If the magnitude of \(\mathbf{x}\) is close to being parallel to \(\mathbf{e}_1\) then \(\mathbf{v}\) will be close to \(\mathbf{0}\) and the calculation of (6.6) will involve the division by a number close to zero. To overcome this we use a Householder transformation to reflect \(\mathbf{x}\) so that it is parallel to \(-\mathbf{e}_1\) instead of \(+\mathbf{e}_1\).

Fig. 6.3 The vector \(\mathbf{x}\) is reflected about the dashed line so that it is parallel to the basis vector \(-\mathbf{e}_1\).

To do this we calculate \(\mathbf{v}\) by adding \(\| \mathbf{x} \| \mathbf{e}_1\) to \( \mathbf{x}\) when the sign of \(x_1\) is negative, i.e.,

\[ \mathbf{v} = \mathbf{x} - \operatorname{sign}(x_1) \| \mathbf{x} \| \mathbf{e}_1, \]

where \(\operatorname{sign}(x)\) returns \(+1\) if \(x \geq 0\) and \(-1\) if \(x < 0\). To compute the QR decomposition of a matrix \(A\), we reflect the first column of \(A\), using Householder reflection so that it is parallel to \(\mathbf{e}_1\). This means that \(r_{11}\) is the only non-zero element in the first column of \(R\).

\[\begin{split} \begin{align*} R = H_1 A = \begin{pmatrix} r_{11} & r_{12} & r_{13} & \cdots & r_{1n} \\ 0 & b_{22} & b_{23} & \cdots & b_{2n} \\ 0 & b_{32} & b_{33} & \cdots & b_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & b_{m2} & b_{m3} & \cdots & b_{mn} \end{pmatrix} \end{align*} \end{split}\]

We then look to do the same for the sub-matrix of \(R\) formed by omitted its first row and column. We do this by calculating a \((m-1) \times (m-1)\) Householder matrix \(H'\) using the second column of \(R\) starting from the diagonal element, i.e., \(\mathbf{x} = (b_{22}, b_{32}, \ldots, b_{m2})^\mathsf{T}\). This is used to form the \(m \times m\) Householder matrix \(H_2\) by padding out the first row and column of the identity matrix.

\[\begin{split} \begin{align*} H_2 = \begin{pmatrix} 1 & \begin{matrix} 0 & \cdots & 0 \end{matrix} \\ \begin{matrix} 0 \\ \vdots \\ 0 \end{matrix} & H' \end{pmatrix} \end{align*} \end{split}\]

The Householder matrix \(H_2\) is then applied to \(R\) so that the elements of the second column below the diagonal are zeroed out.

\[\begin{split} \begin{align*} R = H_2 R = \begin{pmatrix} r_{11} & r_{12} & r_{13} & \cdots & r_{1n} \\ 0 & r_{22} & r_{23} & \cdots & r_{2n} \\ 0 & 0 & c_{33} &\cdots & c_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & c_{n3} & \cdots & c_{nn} \end{pmatrix} \end{align*} \end{split}\]

We repeat this process until column \(j = m - 1\) when \(R\) will be an upper triangular matrix. If \(A\) has \(n\) columns then

\[R = H_{m-1}\ldots H_2 H_1 A,\]

rearranging gives

\[\begin{align*} H_1^{-1} H_2^{-1} \ldots H_{m-1}^{-1} R = A. \end{align*} \]

The Householder matrices are symmetric and orthogonal so \(H^{-1} = H\) and since we want \(A = QR\) then

\[ \begin{align*} Q = H_1H_2 \ldots H_{m-1}. \end{align*} \]

Since each Householder matrix is orthogonal, their product is also orthogonal. Therefore \(Q\) is automatically orthogonal.

Algorithm 6.7 (QR decomposition using Householder reflections)

Inputs: An \(m \times n\) matrix \(A\).

Outputs: An \(n \times n\) orthogonal matrix \(Q\) and and \(n \times n\) upper triangular matrix \(R\) such that \(A = QR\).

• \(Q \gets I_m\)

• \(R \gets A\)

• for \(j = 1, \ldots, \min(m - 1, n)\)

  • \(\mathbf{x} \gets\) column \(j\) of \(R\) starting at the \(j\)-th row

  • \(\mathbf{v} \gets \mathbf{x} - \operatorname{sign} (x_1) \|\mathbf{x}\|\mathbf{e}_1\)

  • \(H' \gets I_{m-j+1} - 2 \dfrac{\mathbf{v} \mathbf{v} ^\mathsf{T}}{\mathbf{v}^\mathsf{T} \mathbf{v}}\)

  • \(H(j:m, j:m) \gets H'\) with the first \(j-1\) rows and columns of \(H\) set to those of \(I_m\)

  • \(R \gets H R\)

  • \(Q \gets Q H\)

• for \(i = 1, \ldots, n\) do

  • If \(R(i,i) < 0\) then

    • \(R(i,: ) \gets -R(i,:)\)

    • \(Q(:,i) \gets -Q(:,i)\)

• Return \(Q\) and \(R\).

Example 6.9

Calculate the QR decomposition of the following matrix using Householder reflections

\[\begin{split} \begin{align*} A=\begin{pmatrix} 1 & -4 \\ 2 & 3 \\ 2 & 2 \end{pmatrix}. \end{align*} \end{split}\]

---

Solution

Set \(Q = I\) and \(R = A\). We want to zero out the elements below \(R(1,1)\) so \(\mathbf{x} = (1, 2, 2)^\mathsf{T}\) and \(\|\mathbf{x}\| = 3\).

\[\begin{split} \begin{align*} j &= 1: & \mathbf{x} &= \begin{pmatrix} 1 \\ 2 \\ 2 \end{pmatrix}, \\ && \mathbf{v} &= \mathbf{x} - \operatorname{sign}(x_1) \|\mathbf{x}\| \mathbf{e}_1 = \begin{pmatrix} 1 \\ 2 \\ 2 \end{pmatrix} - 3 \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} = \begin{pmatrix} -2 \\ 2 \\ 2 \end{pmatrix}, \\ && H &= I - 2 \frac{\mathbf{v} \mathbf{v}^\mathsf{T}}{\mathbf{v}^\mathsf{T} \mathbf{v}} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} - \frac{2}{12} \begin{pmatrix} 4 & -4 & -4 \\ -4 & 4 & 4 \\ -4 & 4 & 4 \end{pmatrix} = \begin{pmatrix} \frac{1}{3} & \frac{2}{3} & \frac{2}{3} \\ \frac{2}{3} & \frac{1}{3} & -\frac{2}{3} \\ \frac{2}{3} & -\frac{2}{3} & \frac{1}{3} \end{pmatrix}, \\ && R &= H R = \begin{pmatrix} \frac{1}{3} & \frac{2}{3} & \frac{2}{3} \\ \frac{2}{3} & \frac{1}{3} & -\frac{2}{3} \\ \frac{2}{3} & -\frac{2}{3} & \frac{1}{3} \end{pmatrix} \begin{pmatrix} 1 & -4 \\ 2 & 3 \\ 2 & 2 \end{pmatrix} = \begin{pmatrix} 3 & 2 \\ 0 & -3 \\ 0 & -4 \end{pmatrix}, \\ && Q &= QH = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \frac{1}{3} & \frac{2}{3} & \frac{2}{3} \\ \frac{2}{3} & \frac{1}{3} & -\frac{2}{3} \\ \frac{2}{3} & -\frac{2}{3} & \frac{1}{3} \end{pmatrix} = \begin{pmatrix} \frac{1}{3} & \frac{2}{3} & \frac{2}{3} \\ \frac{2}{3} & \frac{1}{3} & -\frac{2}{3} \\ \frac{2}{3} & -\frac{2}{3} & \frac{1}{3} \end{pmatrix}. \\ j &= 2: & \mathbf{x} &= \begin{pmatrix} -3 \\ -4 \end{pmatrix}, \\ && \mathbf{v} &= \mathbf{x} - \operatorname{sign}(x_1) \|\mathbf{x}\| \mathbf{e}_1 = \begin{pmatrix} -3 \\ -4 \end{pmatrix} + 5 \begin{pmatrix} 1 \\ 0 \end{pmatrix} = \begin{pmatrix} 2 \\ -4 \end{pmatrix}, \\ && H' &= I - 2 \frac{\mathbf{v}\mathbf{v}^\mathsf{T}}{\mathbf{v}^\mathsf{T} \mathbf{v}} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} - \frac{2}{20} \begin{pmatrix} 4 & -8 \\ -8 & 16 \end{pmatrix} = \begin{pmatrix} \frac{3}{5} & \frac{4}{5} \\ \frac{4}{5} & -\frac{3}{5} \end{pmatrix}, \\ && H &= \begin{pmatrix} 1 & 0 & 0 \\ 0 & \frac{3}{5} & \frac{4}{5} \\ 0 & \frac{4}{5} & -\frac{3}{5} \end{pmatrix}, \\ && R &= H R = \begin{pmatrix} 1 & 0 & 0 \\ 0 & -\frac{3}{5} & -\frac{4}{5} \\ 0 & -\frac{4}{5} & \frac{3}{5} \end{pmatrix} \begin{pmatrix} 3 & 2 \\ 0 & 3 \\ 0 & -4 \end{pmatrix} = \begin{pmatrix} 3 & 2 \\ 0 & -5 \\ 0 & 0 \end{pmatrix}, \\ && Q &= Q H = \begin{pmatrix} \frac{1}{3} & \frac{2}{3} & \frac{2}{3} \\ \frac{2}{3} & \frac{2}{3} & -\frac{1}{3} \\ \frac{2}{3} & -\frac{1}{3} & \frac{2}{3} \end{pmatrix} \begin{pmatrix} 1 & 0 & 0 \\ 0 & \frac{3}{5} & \frac{4}{5} \\ 0 & \frac{4}{5} & -\frac{3}{5} \end{pmatrix} = \begin{pmatrix} \frac{1}{3} & \frac{14}{15} & \frac{2}{15} \\ \frac{2}{3} & -\frac{1}{3} & \frac{2}{3} \\ \frac{2}{3} & -\frac{2}{15} & -\frac{11}{15} \end{pmatrix}. \end{align*} \end{split}\]

The second diagonal element of \(R\) is negative, so we multiply row 2 of \(R\) and column 2 of \(Q\) by \(-1\)

\[\begin{split} \begin{align*} Q &= \begin{pmatrix} \frac{1}{3} & -\frac{14}{15} & \frac{2}{15} \\ \frac{2}{3} & \frac{1}{3} & \frac{2}{3} \\ \frac{2}{3} & \frac{2}{15} & -\frac{11}{15} \end{pmatrix}, & R &= \begin{pmatrix} 3 & 2 \\ 0 & 5 \\ 0 & 0 \end{pmatrix}. \end{align*} \end{split}\]

We can check whether this is correct by verifying that \(A = QR\) and \(Q^\mathsf{T} Q = I\).

\[\begin{split} \begin{align*} QR &= \begin{pmatrix} \frac{1}{3} & -\frac{14}{15} & \frac{2}{15} \\ \frac{2}{3} & \frac{1}{3} & \frac{2}{3} \\ \frac{2}{3} & \frac{2}{15} & -\frac{11}{15} \end{pmatrix} \begin{pmatrix} 3 & 2 \\ 0 & 5 \\ 0 & 0 \end{pmatrix} = \begin{pmatrix} 1 & -4 \\ 2 & 3 \\ 2 & 2 \end{pmatrix} = A, \\ Q^\mathsf{T} Q &= \begin{pmatrix} \frac{1}{3} & \frac{2}{3} & \frac{2}{3} \\ -\frac{14}{15} & \frac{1}{3} & \frac{2}{15} \\ \frac{2}{15} & \frac{2}{3} & -\frac{11}{15} \end{pmatrix} \begin{pmatrix} \frac{1}{3} & -\frac{14}{15} & \frac{2}{15} \\ \frac{2}{3} & \frac{1}{3} & \frac{2}{3} \\ \frac{2}{3} & \frac{2}{15} & -\frac{11}{15} \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} = I. \end{align*} \end{split}\]

---

6.4.3. Loss of orthogonality

Like any numerical technique, the QR decomposition methods shown here are prone to computational rounding errors. If \(Q\) is perfectly orthogonal then \(Q^\mathsf{T}Q = I\), but if there is bound to be some small error, i.e.,

\[ Q^\mathsf{T}Q = I + E \]

where \(E\) is the orthogonality error matrix. A single value for the loss of orthogonality can be computed using

\[ \| E \|_F = \| Q^\mathsf{T}Q - I \|_F, \]

where \(\| \cdot \|_F\) is the Frobenius norm defined by

\[ \| A \|_F = \sqrt{ \displaystyle\sum_{i=1}^m \displaystyle\sum_{j=1}^n |a_{ij}|^2}. \]

The both QR decomposition methods shown here use an iterative approach where we orthogonalised the columns of the \(A\) matrix in turn based upon the previous iterations. So an error in column \(i\) will have a knock on effect on all subsequent columns. To analyse the methods, the QR decomposition of a random \(200 \times 200\) matrix \(A\) was computed using the Gram-Schmidt process and Householder reflections. The Frobenius norm was calculated for the first \(n\) columns of \(Q\) where \(n = 1 \ldots 200\) giving the loss of orthogonality for those iterations and plotted in below.

Fig. 6.4 Comparisons between the loss of orthogonality using Gram-Schmidt and Householder reflections to compute the QR decomposition of a 200 column matrix.

The loss of orthogonality grows much more slowly for Householder reflections than for classical Gram-Schmidt. For this reason, Householder QR is the method used in most numerical linear algebra software.

6.4.4. Calculating eigenvalues using QR decomposition

Eigenvalues and eigenvectors feature prominently in the study of numerical methods for ODEs. Given a system of ODEs, the eigenvalues of the coefficient matrix provide information about the stability, divergence, oscillatory behavior, and constant solutions of the system.

Definition 6.6 (Eigenvalue)

Let \(A\) be an \(n \times n\) matrix. A scalar \(\lambda\) is an eigenvalue of \(A\) if there exists a non-zero vector \(\mathbf{v}\) such that

(6.7)\[A \mathbf{v} = \lambda \mathbf{v}.\]

The vector \(\mathbf{v}\) is called the eigenvector associated with the eigenvalue \(\lambda\).

Rearranging equation (6.7) we have

\[\begin{split} \begin{align*} A \mathbf{v} &= \lambda \mathbf{v} \\ (A - \lambda I) \mathbf{v} &= \mathbf{0}. \end{align*} \end{split}\]

Since \(\mathbf{v} \ne \mathbf{0}\), the homogeneous system \((A - \lambda I) \mathbf{v} = \mathbf{0}\) has non-trivial solution. This occurs if and only if \(A - \lambda I\) is singular, which means

(6.8)\[ \det(A - \lambda I) = 0.\]

For example, consider the eigenvalues of the matrix

\[\begin{split} A = \begin{pmatrix} 2 & 1 \\ 2 & 3 \end{pmatrix}.\end{split}\]

Using equation (6.8)

\[\begin{split} \begin{align*} \det \begin{pmatrix} 2 - \lambda & 1 \\ 2 & 3 - \lambda \end{pmatrix} &= 0 \\ \lambda^2 - 5 \lambda + 4 &= 0 \\ (\lambda - 1)(\lambda - 4) &= 0 \end{align*} \end{split}\]

so the eigenvalues are \(\lambda_1 = 4\) and \(\lambda_2 = 1\) [1]. The problem with using determinants to calculate eigenvalues is that it is too computationally expensive for larger matrices.

6.4.4.1. The QR algorithm

Although the characteristic polynomial can be used to determine eigenvalues analytically for small matrices, it becomes computationally expensive and numerically unreliable for large matrices. Modern numerical software therefore uses iterative methods such as the QR algorithm.

The QR algorithm is a method of computing the eigenvalues of a square matrix \(A\). Let

\[ A_0 = A, \]

and for \(k = 0, 1, \ldots, \)

\[\begin{split} \begin{align*} A_k &= Q_k R_k,\\ A_{k+1} &= R_k Q_k. \end{align*} \end{split}\]

Since \(Q_k\) is orthogonal then \(Q^{-1} = Q^T\) and

\[ A_{k+1} = R_kQ_k = Q_k^{-1} (Q_kR_k) Q_k = Q^{-1}A_kQ_k = Q_k^\mathrm{T}A_kQ_k. \]

Since

\[ A_{k+1} = Q_k^\mathsf{T} A_k Q_k, \]

each iterate is similar to the previous iterate. Therefore all matrices \(A_0, A_1, A_2, \ldots \) have exactly the same eigenvalues. Under suitable conidtions, \(A_k\) converges to an upper triangular matrix.

\[\begin{split} A_{k \to \infty} = \begin{pmatrix} \lambda_1 & \star & \cdots & \star \\ 0 & \lambda_2 & \ddots & \vdots \\ \vdots & \ddots & \ddots & \star \\ 0 & \cdots & 0 & \lambda_n \end{pmatrix}. \end{split}\]

Theorem 6.4 (QR algorithm)

Let \(A_0 = A\) and suppose

\[ A_k = Q_kR_k, \]

is the QR decomposition of \(A_k\). Define

\[ A_{k+1} = R_kQ_k. \]

Then

\[ A_{k+1} = Q_k^\mathsf{T} A_k Q_k, \]

so \(A_{k+1}\) is similar to \(A_k\) and therefore has the same eigenvalues.

Algorithm 6.8 (The QR algorithm)

Inputs: An \(n \times n\) matrix \(A\) and an accuracy tolerance \(tol\).

Outputs: A vector \((\lambda_1, \lambda_2, \ldots, \lambda_n)\) containing the eigenvalues of \(A\).

• for \(k = 1, 2, \ldots\)

  • Calculate the QR decomposition of \(A\)

  • \(A \gets R Q\)

  • If \(\max(|a_{ij}| \text{ for } i > j) < tol\)

    • Break

• Return \(\operatorname{diag}(A)\)

Example 6.10

Use the QR algorithm to compute the eigenvalues of the matrix

\[\begin{split} A = \begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix}, \end{split}\]

using an accuracy tolerance of \(tol = 10^{-4}\).

---

Solution

Calculate the QR decomposition of \(A_0\)

\[\begin{split} \begin{align*} Q_0 &= \begin{pmatrix} 0.3162 & 0.9487 \\ 0.9487 & -0.3162 \end{pmatrix}, \\ R_0 &= \begin{pmatrix} 3.1623 & 4.4272 \\ 0 & 0.6325 \end{pmatrix} \end{align*} \end{split}\]

and calculate \(A_1 = R_0Q_0\)

\[\begin{split} \begin{align*} A_1 &= \begin{pmatrix} 3.1623 & 4.4272 \\ 0 & 0.6325 \end{pmatrix} \begin{pmatrix} 0.3162 & 0.9487 \\ 0.9487 & -0.3162 \end{pmatrix} \\ &= \begin{pmatrix} 5.2 & 1.6 \\ 0.6 & -0.2 \end{pmatrix} \end{align*} \end{split}\]

Calculate the QR decomposition of \(A_1\)

\[\begin{split} \begin{align*} Q_1 &= \begin{pmatrix} 0.9943 & 0.1146 \\ 0.1146 & -0.9934 \end{pmatrix}, \\ R_1 &= \begin{pmatrix} 5.2345 & 1.5665 \\ 0 & 0.3821 \end{pmatrix} \end{align*} \end{split}\]

and calculate \(A_2 = R_1Q_1\)

\[\begin{split} \begin{align*} A_2 &= \begin{pmatrix} 5.2345 & 1.5665 \\ 0 & 0.3821 \end{pmatrix} \begin{pmatrix} 0.9943 & 0.1146 \\ 0.1146 & -0.9934 \end{pmatrix} \\ &= \begin{pmatrix} 5.3796 & -0.9562 \\ 0.0438 & -0.3796 \end{pmatrix}. \end{align*} \end{split}\]

Continuing the iteration gives the sequence of diagonal estimates shown below.

\(k\)	\(\lambda_1\)	\(\lambda_2\)	\(\max(|a_{ij}|, i > j)\)
0	1.000000	4.000000	3.000000
1	5.200000	-0.200000	0.600000
2	5.379562	-0.379562	0.043796
3	5.371753	-0.371753	0.003026
4	5.372318	-0.372318	0.000210
5	5.372279	-0.372279	0.000015

Using the characteristic polynomial to compute the actual eigenvalues

\[\begin{split} \begin{align*} \det \begin{pmatrix} 1 - \lambda & 2 \\ 3 & 4 - \lambda \end{pmatrix} = \lambda^2 - 5 \lambda - 2, \end{align*} \end{split}\]

so \(\lambda_1 = (5 + \sqrt{33})/2 \approx 5.372281\) and \(\lambda_2 = (5 - \sqrt{33}) / 2 \approx -0.372281\) which shows the QR algorithm has calculated the eigenvalues correct to four decimal places.

---

[1]

Eigenvalues are usually listed in descending order by their absolute values.