41.3: Overdetermined Systems

Last updated
Save as PDF

Page ID: 70550

\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \) \( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)\(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\) \(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\)\(\newcommand{\AA}{\unicode[.8,0]{x212B}}\)

import numpy as np
import sympy as sym

When we have more equations than unknowns, we have the overdetermined system \(Ax=b\). In this assignment, we assume that the matrix \(A\) has full column rank. Therefore, the system may not be feasible, i.e., we cannot 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} \nonumber \]

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 =  
## Your code ends ##
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.