QR Decomposition

The QR (orthogonal triangular) decomposition is the most effective and widely used method to find all the eigenvalues of a general matrix. It decomposes the matrix into a normal orthogonal matrix \(Q\) and an upper triangular matrix \(R\), so it is called the QR decomposition. There are two common methods of QR decomposition, one is based on Gram-Schmidt, and the other one is based on Householder reflectors.

QR decomposition based on Gram-Schmidt

Let say a matrix \(A=(a_{1} \quad a_{2} \quad a_{3})\) for example, then \(Ax=(a_{1} \quad a_{2} \quad a_{3})x=b\). To find the \(Q\) and \(R\) for \(A\), we need to orthogonalize and normalize every vector in \(A\) (* for the first vector \(a_1\), only normalization is required).

1. Normalize \(a_1\): \(\quad g_{1}=\frac{a_{1}}{\parallel a_{1} \parallel}, \quad a_{1}=\left\|a_{1}\right\| g_{1}\)

2. Orthogonalize \(a_2\): \(\quad a_{2}^{\prime}=a_{2}-(a_{2}^{\top}g_{1})g_{1}\)

3. Normalize \(a_2\): \(\quad g_{2}=\frac{a_{2}^{\prime}}{\left\|a_{2}^{\prime}\right\|}, \quad a_{2}=\left(a_{2}^{\top} g_{1}\right) g_{1}+\left\|a_{2}^{\prime}\right\| g_{2}\)

4. Orthogonalise \(a_{3}\): \(\quad a_{3}^{\prime}=a_{3}-\left(a_{3}^{\top} g_{1}\right) g_{1}-\left(a_{3}^{\top} g_{2}\right) g_{2}\)

5. Normalize \(a_{3}\): \(\quad g_{3}=\frac{a_{3}{ }^{\prime}}{\left\|a_{3}\right\|}, \quad a_{3}=\left(a_{3}^{\top} g_{1}\right) g_{1}+\left(a_{3}^{\top} g_{2}\right) g_{2}+\left\|a_{3}^{\prime}\right\| g_{3}\)

After the orthogonalization and normalization, we get \(Q=\left(\begin{array}{lll} g_{1} & g_{2} & g_{3} \end{array}\right)\) and \(R=\left(\begin{array}{ccc} \left\|a_{1}\right\| & a_{2}^{\top} g_{1} & a_{3}^{\top} g_{1} \\ 0 & \left\|a_{2}\right\| & a_{3}^{\top} g_{2} \\ 0 & 0 & \left\|a_{3}\right\| \end{array}\right)\)

One problem with Gram-Schmidt is that \(g_{1}, g_{2}, g_{3}\) won't be exactly orthogonal due to error in the computations (round-off).

And as \(g_3\) depends on \(g_1\) and \(g_2\), the problem tends to get worse and worse. The vectors get less and less orthogonal in practice. Therefore, it's better to use Householder reflectors when we apply QR decomposition.

QR decomposition based on Householder reflectors

Householder reflector: \(H=I-2uu^{\top}=I-\frac{2 v v^{\top}}{\|v\|_{2}{ }^{2}}, \quad H^{2}=I\)

For every reflector, it is orthogonal and doesn't change length \(\rightarrow H^{2}=H^{\top}H=I\)

Now consider \(A\) as a \(4*3\) matrix for example, then \(Q_{n}\) \((n=1,2,3)\), all are square matrix with size of \(4*4\). We can find \(Q_{n}\) one by one following the steps below.

Find \(Q_1\) that makes \(Q_{1}A=\left(\begin{array}{ccc}x & x & x \\ 0 & x & x \\ 0 & x & x \\ 0 & x & x\end{array}\right)\):

\(\quad v=a_{1}+\left\|a_{1}\right\|e_{1} \quad\) (\(a_{1}\) is the first column in \(A\)), \(\quad Q_{1}=I - \frac{2 v v^{\top}}{\|v\|_{2}{ }^{2}}\)

Find \(Q_2\) that makes \(Q_{2}(Q_{1}A)=\left(\begin{array}{ccc}x & x & x \\ 0 & x & x \\ 0 & 0 & x \\ 0 & 0 & x\end{array}\right)\):

\(\quad v=a_{2}+\left\|a_{2}\right\|e_{1} \quad\) (\(a_{2}\) is last three elements of the second column in \(Q_{1}A\))

\(\quad Q_{2}^{\prime}=I - \frac{2 v v^{\top}}{\|v\|_{2}{ }^{2}}\), \(\quad Q_{2}=\left(\begin{array}{cccc}1 & 0 & 0 & 0 \\ 0 & Q_{2_{1,1}}^{\prime} & Q_{2_{1,2}}^{\prime} & Q_{2_{1,3}}^{\prime} \\ 0 & Q_{2_{2,1}}^{\prime} & Q_{2_{2,2}}^{\prime} & Q_{2_{2,3}}^{\prime} \\ 0 & Q_{2_{3,1}}^{\prime} & Q_{2_{3,2}}^{\prime} & Q_{2_{3,3}}^{\prime} \end{array}\right)\)

Find \(Q_3\) that makes \(Q_{3}(Q_{2}Q_{1}A)=\left(\begin{array}{ccc}x & x & x \\ 0 & x & x \\ 0 & 0 & x \\ 0 & 0 & 0\end{array}\right)\):

\(\quad v=a_{3}+\left\|a_{3}\right\|e_{1} \quad\) (\(a_{3}\) is last two elements of the third column in \(Q_{2}(Q_{1}A)\))

\(\quad Q_{3}^{\prime}=I - \frac{2 v v^{\top}}{\|v\|_{2}{ }^{2}}\), \(\quad Q_{3}=\left(\begin{array}{cccc}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & Q_{3_{1,1}}^{\prime} & Q_{3_{1,2}}^{\prime} \\ 0 & 0 & Q_{3_{2,1}}^{\prime} & Q_{3_{2,2}}^{\prime} \end{array}\right)\)

\(Q_{3} Q_{2} Q_{1} A=R \rightarrow A=Q_{1}^{\top} Q_{2}^{\top} Q_{3}^{\top} R=Q R\), and we get \(Q=Q_{1}^{\top} Q_{2}^{\top} Q_{3}^{\top}\)

QR Decomposition in Python

import math
import argparse
import numpy as np
from typing import Union

# QR decomposition based on Gram-Schmidt
def gram_schmidt(A):
    cols = A.shape[1]
    Q = np.copy(A)
    R = np.zeros((cols, cols))
    for col in range(cols):
        for i in range(col):
            k =  np.sum(a[:, col] * Q[:, i]) / np.sum( np.square(Q[:, i]) )
            Q[:, col] -= k*Q[:, i]
        Q[:, col] /= np.linalg.norm(Q[:, col])
        for i in range(cols):
            R[col, i] = Q[:, col].dot( A[:, i] )   
    return Q, R  

# QR decomposition based on Householder reflectors  
def householder(alpha: float, x: np.ndarray) -> Union[np.ndarray, int]:
    s = math.pow(np.linalg.norm(x, ord=2), 2)
    v = x
    if s == 0:
        tau = 0
    else:
        t = math.sqrt(alpha * alpha + s)
        v_one = alpha - t if alpha <= 0 else -s / (alpha + t)
        tau = 2 * v_one * v_one / (s + v_one * v_one)
        v /= v_one
    return v, tau

def qr_decomposition(A: np.ndarray, m: int, n: int) -> Union[np.ndarray, np.ndarray]:
    H = []
    R = A
    Q = A
    I = np.eye(m, m)
    for j in range(0, n):
        # Apply Householder transformation.
        x = A[j + 1:m, j]
        v_householder, tau = householder(np.linalg.norm(x), x)
        v = np.zeros((1, m))
        v[0, j] = 1
        v[0, j + 1:m] = v_householder
        res = I - tau * v * np.transpose(v)
        R = np.matmul(res, R)
        H.append(res)
    return Q, R