In-Class Assignment: Solve Linear Systems of Equations using QR Decomposition¶
Image From: https://en.wikipedia.org/wiki/Hydra
In order to successfully complete this assignment you need to participate both individually and in groups during class. If you attend class in-person then have one of the instructors check your notebook and sign you out before leaving class. If you are attending asyncronously, turn in your assignment using D2L no later than 11:59pm on the day of class. See links at the end of this document for access to the class timeline for your section.
Image From: https://en.wikipedia.org/wiki/Hydra
</p>
import numpy as np
import sympy as sym
In this assignment, we try to solve the linear systems $Ax = b$ in three different categories.
When we have the same number of equations as unknowns, we have the following system $$Ax= b$$ with a squre matrix $A$. In this section, we assume that the matrix $A$ has full rank. That is the system has an unique solution. We talked about many ways to solve this system of equations. Some examples are:
In this assignment, we will show that we can solve it by QR decomposion.
Consider the following system of equations:
$$\begin{bmatrix}5&-2&2 \\ 4 & -3 &4 \\ 4& -6 &7 \end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}=\begin{bmatrix}1\\2\\3\end{bmatrix}$$A = np.matrix([[5, -2, 2], [4, -3, 4], [4,-6,7]])
b = np.matrix([[1],[2],[3]])
display(sym.Matrix(A))
display(sym.Matrix(b))
Back substitution. Let's first implement the back substitution in Python. The back substitution function back_subst solves the system
$$R x = b$$
where $R$ is an upper-triangular matrix.
When we solve for $x$, we start with $x_n$: $$x_n = {b_n\over R_{n,n}}$$ Then we solve for $x_{n-1}$ as $$x_{n-1} = {b_{n-1}-R_{n-1,n}x_n\over R_{n-1,n-1}}$$ $$x_{n-2} = {b_{n-2}-R_{n-2,n-1}x_{n-1}-R_{n-2,n}x_n\over R_{n-2,n-2}}$$ Then we can find $x_{n-2},\cdots,x_1$. We can solve for all components of $x$ in the reserved order. So this is call back substitution.
✅ **DO THIS:** Complete the following code for back substitution.
def back_subst(R,b):
n = R.shape[0]; x = np.zeros(n);
for i in reversed(range(n)):
x[i] = b[i]
for j in range(i+1,n):
## Your code starts ##
x[i] = # Complet this line of code.
## Your code ends ##
x[i] = x[i]/R[i,i]
return np.matrix(x).T
File "<ipython-input-1-fbd4caaeb949>", line 7 x[i] = # Complet this line of code. ^ SyntaxError: invalid syntax
✅ **DO THIS:** When we solve for $Ax=b$ with QR decomposition. We have the following steps:
Use these steps to solve $Ax=b$ with the given $A$ and $b$. Compare the result with np.linalg.solve.
## Your code starts ##
x =
## Your code ends ##
print(type(x)) # x is a column vector in the np.matrix type
np.allclose(x, np.linalg.solve(A,b))
File "<ipython-input-1-3ebbe42bac76>", line 2 x = ^ SyntaxError: invalid syntax
When we have more equations than unknowns, we have the overdetermined system $$Ax\approx b$$ In this assignment, we assume that the matrix $A$ has full column rank. Therefore, the system may not be feasible, i.e., we can not find $x$ such that $Ax=b$. Then we want to find the $x$ to minimize $\|Ax-b\|^2$, which is the least squares problem. We mentioned in previous assignments that we can solve this least squares problem by finding the left inverse of the matrix $A$. That is $$x=(A^\top A)^{-1}A^\top b$$
In this assignment, we show that we can solve it by QR decomposion.
Consider the following system of equations:
$$\begin{bmatrix}5&-2&2 \\ 4 & -3 &4 \\ 4& -6 &7 \\ 6 & 3 & -3\end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}=\begin{bmatrix}1\\2\\3\\-1\end{bmatrix}$$✅ **DO THIS:** We solve the least squares problem in the following steps
back_subst function we defined before to solve $Rx = Q^\top b$A = np.matrix([[5, -2, 2], [4, -3, 4], [4,-6,7], [6,3,-3]])
b = np.matrix([[1],[2],[3],[-1]])
display(sym.Matrix(A))
display(sym.Matrix(b))
## Your code starts ##
x =
## Your code ends ##
print(type(x)) # x is a column vector in the np.matrix type
print(x)
File "<ipython-input-1-df614842c29f>", line 2 x = ^ SyntaxError: invalid syntax
We can not use the np.linalg.solve because the matrix $A$ is not a square matrix (you can try if you do not believe it). However, we can use the np.linalg.lstsq function to find the least squares solution to minimize $\|Ax-b\|^2$. The next cell compares the result from lstsq and our result from the QR decomposition. (If everything is correct, you will expect a True result.)
np.allclose(x, np.linalg.lstsq(A,b)[0])
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-1-58aca6ede5bb> in <module> ----> 1 np.allclose(x, np.linalg.lstsq(A,b)[0]) NameError: name 'x' is not defined
✅ **DO THIS:** Explain why we can use the QR decomposition to solve the least squares problem.
Put Your Answer Here
When we have more unknowns than equations, we have the underdeterminated system $$Ax= b$$ In this assignment, we assume that the matrix $A$ has full row rank. This system has infinite many solution, i.e., we can not find an unique $x$ such that $Ax=b$. Then we want to find the smallest $x$ (by smallest, we mean the smallest $\|x\|^2$) such that $Ax=b$, which is also the least squares problem.
In this assignment, we show that we can also solve it by QR decomposion.
Consider the following system of equations:
$$\begin{bmatrix}5&-2&2 & 1 \\ 4 & -3 &4 &2 \\ 4& -6 &7 & 4\end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\\x_4\end{bmatrix}=\begin{bmatrix}1\\2\\3\end{bmatrix}$$✅ **DO THIS:** We solve the least squares problem in the following steps
forward_subst function defined below to solve $x = Q (R^\top)^{-1}b$A = np.matrix([[5, -2, 2, 1], [4, -3, 4, 2], [4,-6,7, 4]])
b = np.matrix([[1],[2],[3]])
display(sym.Matrix(A))
display(sym.Matrix(b))
def forward_subst(L,b): # This function solves $L x= b$ when $L$ is the lower-trigular matrix
n = L.shape[0]; x = np.zeros(n);
for i in range(n):
x[i] = b[i]
for j in range(i):
x[i] = x[i] - L[i,j]*x[j]
x[i] = x[i]/L[i,i]
return np.matrix(x).T
## Your code starts ##
x =
## Your code ends ##
print(type(x)) # x is a column vector in the np.matrix type
print(x)
File "<ipython-input-1-df614842c29f>", line 2 x = ^ SyntaxError: invalid syntax
We can not use the np.linalg.solve because the matrix $A$ is not a square matrix. However, we can use the np.linalg.lstsq function to find the least squares solution to minimize $\|Ax-b\|^2$ with underdeterminated systems. The next cell compares the result from lstsq and our result from the QR decomposition. (If everything is correct, you will expect a True result.)
np.allclose(x, np.linalg.lstsq(A,b)[0])
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-1-58aca6ede5bb> in <module> ----> 1 np.allclose(x, np.linalg.lstsq(A,b)[0]) NameError: name 'x' is not defined
✅ **DO THIS:** Explain why we can use the QR decomposition to solve the least squares problem. (HINT: you may need the orhogonal decomposition to the four fundamental spaces of $Q$.)
Put Your Answer Here
If you attend class in-person then have one of the instructors check your notebook and sign you out before leaving class. If you are attending remote, turn in your assignment using D2L.
Written by Dr. Dirk Colbry, Michigan State University
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.