41.4: Underdetermined Systems

import numpy as np
import sympy as sym

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 infinitely many solutions, i.e., we cannot find a 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} \nonumber \]

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 =  
## 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. 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\).)