# Direct Inversion via LU Decomposition

- Motivation
- Mechanics of solving a least squares problem
- First perform LU decomposition
- Now let's applied what we've learned to an interpolation example
- Compare results

# Motivation

- this tutorial is also available as a colab notebook

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.

```
cd 'drive/My Drive/Colab Notebooks/Solving Least Squares Problems'
```

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 decomposition

https://en.wikipedia.org/wiki/LU_decomposition

It 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/0898713617

LU decomposition is essentially Gaussian elimination to find the upper matrix while recording your operations in the lower matrix

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

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:

- solve for $y$ in $Ly = b$
- then solve for $x$ after solving for $y$ in $Ux = y$

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

```
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
```

Finally let's check against numpy's built in least squares solver...

```
print(np.linalg.lstsq(A, b, rcond = None)[0])
```

Yay we get the same answer as numpy's least squares solver!

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

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

```
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()
```

```
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()
```

Move on to the next notebook on applying Vandermonde interpolation to a beam hardening problem in CT imaging: https://colab.research.google.com/drive/1AMKoOjs7y6l18d70RcbHMLU9ZFxkzRpO