<jupyter_start><jupyter_text>Direct Inversion via LU Decomposition- toc: true - hide: false - badges: true- comments: true- categories: [machine learning]  Motivation- this tutorial is also available as a [colab notebook](https://colab.research.google.com/drive/1SJtGSP_1XcjLh0pyF_RwmakUMdDuNcxt?usp=sharing)Continuing with the Vandermonde interpolation example this tutorial goes step by step on how computers typically implement matrix inversions for the purpose of solving least squares problems.This typically takes the form of a matrix factorization of the matrix to be inverted followed by forward and back substitution. The matrix factorization demoed here is LU decomposition. We compare our hand-coded solvers with Numpy library solvers to compare results.<jupyter_code>cd 'drive/My Drive/Colab Notebooks/Solving Least Squares Problems'<jupyter_output>[Errno 2] No such file or directory: 'drive/My Drive/Colab Notebooks/Solving Least Squares Problems'
/content<jupyter_text>Mechanics of solving a least squares problem Introduction to LU decomposition  Our goal is to explore how to solve systems of linear equations following Trefethan and fastai's numerical linear algebrea course: https://www.youtube.com/watch?v=O2x5KPJr5ag&list=PLtmWHNX-gukIc92m1K0P6bIOnZb-mg0hY&index=5.  First perform LU decompositionhttps://en.wikipedia.org/wiki/LU_decompositionIt is also based on Trefethen ch. 20 on LU decompositions and ch. 17 on Backsubstitution.https://www.amazon.com/Numerical-Linear-Algebra-Lloyd-Trefethen/dp/0898713617LU decomposition is essentially Gaussian elimination to find the upper matrix while recording your operations in the lower matrix<jupyter_code>import numpy as np
def LU(A):
    U = np.copy(A)
    m, n = A.shape
    L = np.eye(n)
    for k in range(n-1):
        for j in range(k+1,n):
            L[j,k] = U[j,k]/U[k,k]
            U[j,k:n] -= L[j,k] * U[k,k:n]
    return L, U
A = np.array([[1,2,-3],[2,-1,1],[1,4,-2]]).astype(np.float)
print(A)
L, U = LU(A)
print("L = ")
print(L)
print("U = ")
print(U)
print("LU = ")
print(L@U)
b = np.array([[1],[1],[9]])
print(b)<jupyter_output>[[1]
 [1]
 [9]]
[[1]
 [1]
 [9]]<jupyter_text>Here L and U must be square. Once we have our matrix factorizations of $A$, then we can solve for $x$ in a 2 step process:1. solve for $y$ in $Ly = b$2. then solve for $x$ after solving for $y$ in $Ux = y$<jupyter_code>def forwardsub(L, b):
    y = np.zeros(len(L))
    y[0] = b[0]
    for i in range(1,len(y)):
        y[i] = b[i] - np.dot(L[i,:],y)
    return y
def backsub(U, y):
    N = len(U)
    x = np.zeros([N,1])
    x[-1] = y[-1]/U[-1,-1] 
    for i in range(N-2, -1, -1):
        x[i] = (y[i] - np.dot(U[i,:],x))/float(U[i,i])
    return x
y = forwardsub(L, b)
x = backsub(U, y)
print(x)
bhat = (A@x)
np.allclose(b,bhat)<jupyter_output><empty_output><jupyter_text>Yay we got back the same solution!  Put into function and compare against NumpyWe can combine all 3 of these functions in a least squares solver..<jupyter_code>def LU_solve(A, b):
    L, U = LU(A)
    y = forwardsub(L, b)
    x = backsub(U, y)
    return x
x = LU_solve(A, b)
print(x)
A@x<jupyter_output><empty_output><jupyter_text>Finally let's check against numpy's built in least squares solver...<jupyter_code>print(np.linalg.lstsq(A, b, rcond = None)[0])<jupyter_output>[[1.]
 [3.]
 [2.]]
[[1.]
 [3.]
 [2.]]<jupyter_text>Yay we get the same answer as numpy's least squares solver!  Now let's applied what we've learned to an interpolation example Our goal is to interpolate the function below<jupyter_code>import matplotlib.pyplot as plt
nSamples = 10
t = np.linspace(0, 10, nSamples)
y = np.exp(-5*(t-3)**2)+0.1*np.random.randn(nSamples)

plt.plot(t,y,'*')
plt.show()
deg=6
A = np.vander(y, N=deg, increasing=True)<jupyter_output><empty_output><jupyter_text>Here we solve the interpolation problem by solving for $x$ in the given equation $y = Ax$, found by taking the pseudo inverse:\begin{equation}(A^*A)^{-1}A^*y = x\end{equation}This works by left multiply $y=Ax$ by $A^*$ to get $A^*y = A^*Ax$ and then noticing that $A^*A$ is a square matrix, which thus can be inverted using the LU decomposition solve we introduced above(If $[A] = [m \times n]$ then $[A*A] = [n \times m] \times [m \times n] = [n \times n]$) Above is how things work mathematically, to actually do this in code we need to look at the problem like this...Our goal is to estimate $y$ by using the Vandermonde matrix and an appropriate $x$ vector such that \begin{equation}y \approx \hat{y} = Ax\end{equation}however to find the appropriate $x$ we must perform least squares fitting....\begin{equation}x = \min ||y-\hat{y}||_2^2 = \min || y - Ax ||_2^2\end{equation}Anyway to solve equation 1 we actually back up a step before finding the inverse\begin{equation}A^*y = A^*A x\end{equation}which we can rewrite as $Y = A^*y$ and $A^*A = M$ to get \begin{equation}Mx = Y\end{equation}and invert to find x exactly like we did above<jupyter_code>Astar = np.conjugate(A.transpose())
print(np.shape(Astar))
print(np.shape(Astar@A))
x = LU_solve(Astar@A, Astar@y)
plt.plot(t,y,'*',t,A@x)
plt.show()<jupyter_output><empty_output><jupyter_text>Compare resultsFinally we put it all into a function a show that this is a solution for Least squares of *rectangular* matrices<jupyter_code>def LUsolveRect(A, y):
    Astar = np.conjugate(A.transpose())
    x = LU_solve(Astar@A, Astar@y)
    return x
x = LUsolveRect(A, y)
print(x)
plt.plot(t,y,'*',t,A@x)
plt.show()
x = np.linalg.lstsq(A, y, rcond=None)[0]
print(x)
plt.plot(t,y,'*',t,A@x)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13927227101631118, 'fraction_numerical': 0.01577343609965944, 'mean_word_length': 4.58, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 30, 'https://': 4, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '358240', 'license': 'Apache-2.0', 'n_tokens_mistral': 2150, 'n_tokens_neox': 2032, 'n_words': 639, 'path': '_notebooks/2021-01-28-LU Solving.ipynb', 'repo_name': 'bnel1201/brandon-j-nelson', 'size': '5579'}