Detailed Explanation of QR Decomposition by Householder Transformation
Episode 1: QR Decomposition by Classical Gram-Schmidt (CGS)and modified Gram-Schmidt (MGS)
Episode 2: QR Decomposition by Householder Transformation
Episode 3: QR Decomposition by Givens Rotation
In Episode 1, we talked about both classical Gram-Schmidt (CGS)and modified Gram-Schmidt (MGS) processes. In this episode (episode 2), we will go through the Householder Transformation method to compute QR decomposition.
If you are not familiar with QR decomposition, these are some key points that you need to know before we start.
“Q” in QR decomposition represents an orthogonal matrix, which means:
Remarks: for an orthogonal matrix, the dot product of any 2 column vectors in orthogonal matrix must be equal to 0
“R” in QR decomposition represents right (upper) triangular matrix, which means:
Therefore, QR decomposition is sometimes called QU decomposition.
Some Characteristics of QR Decomposition
- QR Decomposition always exists, but may not be unique
- If the determinant of matrix A is not 0 and all diagonal entries of R > 0, then this QR decomposition is unique.
- There are two types of QR decomposition:
a. Full QR Decomposition
b. Reduced QR Decomposition
For simplicity, we will only include full QR decomposition for the square matrix here, whereas matrix A is full rank (i.e. All columns are linearly independent).
QR decomposition is also applicable for m x n (m-by-n) matrixes, where m ≥ n.
In that case, Q will be a m x m (m-by-m) matrix and R will be an m x n (m-by-n) matrix.
4 Methods to perform QR decomposition
- (Classical) Gram-Schmidt Process (CGS)
- Modified Gram-Schmidt Process (MGS)
- Household Transformation
- Givens Rotation
In this article, we will cover Household Transformation.
Set-Up
Like the last episode, we will go through the steps of QR decomposition and implementation of QR decomposition with Python code.
We created a Python class called Matrix, which will be used in our example. First, we create a 3-by-3 Matrix A for our QR decomposition example.
# main.py
# Create the matrix A
A = Matrix(
n_row=3,
n_col=3,
two_d_array=[[2, -1, -2], [-4, 6, 3], [-4, -2, 8]],
)
QR Decomposition by Householder Transformation
We assume you already know what vector projection is. If you do not, you may refer to the last episode.
Householder Transformation
In Householder transformation, similar to Gram-Schmidt, we try to decompose each column vector in A to a set of linear combinations of orthogonal vectors in Q.
Therefore, we will try to map each column to a set of orthogonal vectors.
In order to map each vector to the orthogonal vectors above, we applied householder reflection to each column vector (or the Matrix A)
To reflect vector a to the direction of e_1 (i.e. a’), we reflect vector a about the direction of vector u, while vector w is orthogonal to u.
From the equation above, we can first calculate the vector w. Then we can obtain the Householder Matrix Q.
Remarks:
For this Matrix Q, it decomposes all column vectors in A for the direction of e_1. Then we will move to the next iteration.
In the next iteration, we will only use the sub-matrix (minor) to compute the next orthogonal vector and so on. Let’s look into our example.
Our Example — First Iteration
In the first iteration, we will use e_1 to calculate the vector w first and compute Q_1. Then compute the multiplication of Q_1 and A for the next iteration.
Caution: The red box in the above image is misaligned. It should be bounding the sub-matrix ((4,-2),(-4,3)).
Our Example — Second iteration.
Caution: The red box in the above image is misaligned. It should be bounding the sub-matrix (sqrt(2) / 2).
Our Example — Third iteration.
Our Example — Result
[Step-by-Step Python Code] QR Decomposition by Householder Transformation:
# QR Decomposition by Householder Transformation
def QR_Householder(M: Matrix) -> tuple():
Qs = []
Q_M = Matrix(n_row=M.n_col, n_col=M.n_col)
Q_M.set_diagonal()
curr_M = M
# ===========================================================
for k in range(0, M.n_row):
# Get the sub-matrix, nothing happened then k = 0
sub_matrix_values = []
for i in range(M.n_row):
row = []
for j in range(M.n_col):
if i < k or j < k:
continue
row.append(curr_M.values[i][j])
if len(row) > 0:
sub_matrix_values.append(row)
# Get Sub-Matrix for next iteration
sub_matrix = Matrix(
n_row=M.n_row - k, n_col=M.n_col - k, two_d_array=sub_matrix_values
)
# ==================== Start Calculation ============================
sign_a11 = 1 if sub_matrix.values[0][0] >= 0 else -1
# Get the current column vector
q = [row[0] for row in sub_matrix.values]
e = [1 if j == 0 else 0 for j in range(sub_matrix.n_col)]
# Normalise the column vector
b = q - (-sign_a11) * sub_matrix.length(q) * np.array(e)
q = b / sub_matrix.length(b)
# Calculate matrix Q
q_values = []
for i in range(sub_matrix.n_row):
row = []
for j in range(sub_matrix.n_col):
if i == j:
row.append(1 - 2 * q[i] * q[j])
else:
row.append(-2 * q[i] * q[j])
q_values.append(row)
Q = Matrix(n_row=sub_matrix.n_row, n_col=sub_matrix.n_col, two_d_array=q_values)
# ===========================================================
# Check if need extend, nothing to do if k = 0
if Q.n_row < M.n_row and Q.n_col < M.n_col:
# Extend the Matrix to original size
dim_this = Q.n_row
dim_first = M.n_row
dim_diff = dim_first - dim_this
if dim_diff > 0:
extended_Q_values = []
for i in range(dim_first):
row = []
for j in range(dim_first):
# Add the extended part
if i < dim_diff or j < dim_diff:
if i == j:
row.append(1)
else:
row.append(0)
# Add back the result of the Q
else:
row.append(Q.values[i - dim_diff][j - dim_diff])
extended_Q_values.append(row)
Q = Matrix(n_row=dim_first, n_col=dim_first, two_d_array=extended_Q_values)
# ===========================================================
# Save the Q matrix
Qs.append(Q)
curr_M = Q.multiply(M)
# Calculate Q and R with all the {Q_1 ... Q_k}
final_R = M
final_Q = Matrix(n_row=M.n_row,n_col=M.n_col)
final_Q.set_diagonal()
for matrix in Qs:
final_R = matrix.multiply(final_R)
matrix.transpose()
final_Q = final_Q.multiply(matrix)
return final_Q, final_R# main.py
# ----- Example of Householder Transformation -----
print("----- Householder Transformation -----")
A.print_values(matrix_name="A")
Q, R = QR_Householder(A)
Q.print_values(matrix_name="Q", num_digits=4)
R.print_values(matrix_name="R", num_digits=4)
print("----- Validate Result -----")
new_A = Q.multiply(R)
new_A.print_values(matrix_name="new_A", num_digits=4)Python Output
Done! We successfully implement Householder Transformation for QR decomposition.
The Python code of the above example can be found on GitHub.
If you like this article, give me a like and share it with your friends. You may also leave me a comment/message too.
