# 21 In-Class Assignment: Solve Linear Systems of Equations using QR Decomposition¶

Image From: https://en.wikipedia.org/wiki/Hydra

## Agenda for today’s class (80 minutes)¶

import numpy as np
import sympy as sym

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-1-3b02787bb8ab> in <module>
1 import numpy as np
----> 2 import sympy as sym

ModuleNotFoundError: No module named 'sympy'


In this assignment, we try to solve the linear systems $$Ax = b$$ in three different categories.

• $$A$$ is a square matrix. Unique solution when $$A$$ is invertible

• overdetermined (more equations than unknowns): If $$A$$ has full column rank, the system has an unique solution when $$b$$ is in the column space of $$A$$, otherwise no solution.

• underdetermined (more unknowns than equations): If $$A$$ has full row rank, there are infinite many solutions.

## 2. Solve Linear Systems¶

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:

• Jacobian iteration/ Gauss-Seidel iteration

• $$x = A^{-1} b$$

• Gaussian elimination

• LU decomposition

In this assignment, we will show that we can solve it by QR decomposion.

Consider the following system of equations:

$\begin{split}\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}\end{split}$
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):
x[i] =   # Complet this line of code.
x[i] = x[i]/R[i,i]
return np.matrix(x).T


DO THIS: When we solve for $$Ax=b$$ with QR decomposition. We have the following steps:

• Find the QR decomposition of $$A$$

• From $$QRx=b$$, we obtain $$Rx =Q^\top b$$

• Solve for $$x$$ using back substitution.

Use these steps to solve $$Ax=b$$ with the given $$A$$ and $$b$$. Compare the result with np.linalg.solve.

## Your code starts ##
x =
print(type(x))   # x is a column vector in the np.matrix type
np.allclose(x, np.linalg.solve(A,b))


## 3. Overdetermined Systems¶

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{split}\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}\end{split}$

DO THIS: We solve the least squares problem in the following steps

• Find the QR decomposition of the matrix $$A$$ such that $$R$$ is a square upper-triangular matrix. ($$Q$$ is not a square matrix any more)

• Use the 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 =
print(type(x))   # x is a column vector in the np.matrix type
print(x)


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])


DO THIS: Explain why we can use the QR decomposition to solve the least squares problem.

## 4. Underdetermined Systems¶

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{split}\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}\end{split}$

DO THIS: We solve the least squares problem in the following steps

• Find the QR decomposition of the matrix $$A^\top$$ such that $$R$$ is a square upper-triangular matrix.

• Use the 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 =
print(type(x))   # x is a column vector in the np.matrix type
print(x)


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])


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$$.)