# Fitting Models to Data with Ax = b

## 1. Distilling our data into a model
We are often faced with the task of fitting a model to data that has been aquired (from experiment or computation). The usual situation is that we have more data points than are needed to determine the coefficients of the model - an **over-determined system** - and we seek some kind of 'best fit' to the data. 

## 2. Ax = b


If our model is a linear combination of terms, we can express it compactly using the matrix equation $\mathbf{Ax} = \mathbf{b}$, where:
 - $\mathbf{A}$ is the **model** matrix. Each row of A is a different data point, each column is a different model term. For example, if we seek a straight line fit to a series of points $(x,y)$ then the first column of A would be a vector of $1$s and the second would be a vector of $x$ values.
 - $\mathbf{x}$, our solution vector, is the vector of **loadings** (coefficients) for each term in our model that best matches our data points.
 - $\mathbf{b}$, our **outcomes** vector. This is made up of each $y$ value from our set of $(x,y)$ data points.
 
In an over-determined system, we have more data points (rows of $\mathbf{A}$) than we have coefficients (entries in $\mathbf{x}$).

## 3. Solving a uniquely determined system using the inverse of A
As an illustration, we start with two data points and we seek a line of best fit! We have two data points and two coefficients for our model, so we have a uniquely determined system. In this case, $\mathbf{A}$ is a square matrix and we can solve $\mathbf{Ax}=\mathbf{b}$ (to find our coefficients vector, $\mathbf{x}$) by inverting $\mathbf{A}$.


```python
import numpy as np 
import matplotlib.pyplot as plt
from sklearn.linear_model import Lasso
%matplotlib inline
```


```python
slope = 3.                             # Our data points lie on a line defined by slope,
intercept = 2.                         # and intercept
x_data = np.array([0,1])
y_data = intercept + slope*x_data
plt.scatter(x_data, y_data, c='r');
```


```python
A = np.ones([2,2])                     # Assemble the A matrix. Start by making all entries 1
A[:,1] = x_data                        # then replace the second column with our independent variable, x_data
coeffs = np.linalg.inv(A) @ y_data     # We call our solution vector x (in Ax=b) coeffs, and find it by inverting A
coeffs
```


```python
y_fit = coeffs[0] + coeffs[1]*x_data   # Generate our model line using coeffs
plt.scatter(x_data, y_data, c='r')
plt.plot(x_data, y_fit);               # Plot the line to check the fit
```

## 4. Solving an over-determined system using the pseudo-inverse of A
We now construct the more common case where we have many more data points than coefficients.


```python
npts = 100                            # Number of data points
x_data = np.linspace(0,1,npts)
y_data = intercept + slope*x_data + 0.2*np.random.normal(size=npts)  # Add noise to the data
plt.scatter(x_data,y_data, c='r', alpha=0.5);
```


```python
A = np.ones([npts,2])                 # We constuct A in the same way as before - our model (linear fit) has not changed
A[:,1] = x_data
coeffs = np.linalg.pinv(A) @ y_data   # We now use the pseudo-inverse of A
coeffs
```


```python
y_fit = coeffs[0] + coeffs[1]*x_data
plt.scatter(x_data, y_data, c='r',alpha = 0.5)
plt.plot(x_data, y_fit, 'b', linewidth = 2);
```

## 5. Relationship between the pseudo-inverse and the SVD
If we obtain the SVD of $\mathbf{A}$ we may write,

\begin{align}
\mathbf{b} &= \mathbf{A}\mathbf{x} \\
&= \mathbf{U}\mathbf{\Sigma}\mathbf{V^T}\:\mathbf{x}\quad\text{,}
\end{align}

and so

\begin{align}
\mathbf{x} &= (\mathbf{U}\mathbf{\Sigma}\mathbf{V^T})^{-1}\:\mathbf{b} \\
&= \mathbf{V}\mathbf{\Sigma^{-1}}\mathbf{U^T}\:\mathbf{b}\quad\text{.}
\end{align}

Since $\mathbf{\Sigma}$ is a diagonal matrix, the inverse is obtained by taking the reciprocal of all the non-zero diagonal entries.

We can see that the 'pseudo-inverse' is $\mathbf{V}\mathbf{\Sigma^{-1}}\mathbf{U^T}$.


```python
U,S,VT = np.linalg.svd(A)              # Obtain the SVD of A
Sinv = np.zeros([2,npts])              # Construct the inverse of S by
Sinv[:2,:2]=np.diag(np.reciprocal(S))  # taking the reciprocal of the entries on the diagonal 
Ainv = VT.T @ Sinv @ U.T               # Invert the SVD
coeffs = Ainv @ y_data                 # Find the coeffs vector 
coeffs
```

## 6. What is the pseudo-inverse doing?
When we use the pseudo-inverse, we are obtaining a 'least squares' fit. Formally, this is written as,

\begin{equation}
\hat{x} = \underset{x}{\text{argmin}}\: ||\mathbf{Ax}-\mathbf{b} ||_2
\end{equation}

where the $\ell_2$ norm means,

\begin{equation}
||\mathbf{a}||_2 = \sqrt{\sum_{i=1}^n\:a_i^2}
\end{equation}


## 7. Over-fitting and sparsity
A property of the $\ell_2$ norm is that it leads to coefficients in all of the model terms and this may not be a desirable property. 

As an example, we fit a higher order polynomial to our set of points: 


```python
npoly = 10                           # The order of the polynomial we are fitting
A = np.zeros([npts, npoly])          # Form our model matrix, A
for n in range(npoly):
    A[:,n] = x_data**n               # Column n of A is our x_data vector raised to the power n
coeffs = np.linalg.pinv(A) @ y_data  # Use the pseudo-inverse to obtain a least squares fit
coeffs
```

Our `coeffs` vector has non-zero values for all of the model terms. This is guaranteed to fit our data with minimal error in an $\ell_2$ norm sense:



```python
y_fit = np.zeros(npts)
for n in range(npoly):
    y_fit = y_fit + coeffs[n]*x_data**n    
plt.scatter(x_data,y_data,c='r',alpha=0.5)
plt.plot(x_data, y_fit, 'b', linewidth=2);
```

But is likely to be unreliable if we extrapolate in $x$:


```python
x_test=np.linspace(-0.2,1.2,npts)          # extend our x rang
y_fit = np.zeros(npts)
for n in range(npoly):
    y_fit = y_fit + coeffs[n]*x_test**n
plt.scatter(x_data,y_data,c='r',alpha=0.5)
plt.plot(x_test, y_fit, 'b', linewidth=2);
```

One way to overcome this is to find our model coefficients by solving a different optimisation problem, 

\begin{equation}
\hat{x} = \underset{x}{\text{argmin}}\: ||\mathbf{Ax}-\mathbf{b} ||_2 + \alpha ||\mathbf{x}||_1 \quad\text{,}
\end{equation}

where the $\ell_1$ norm means,

\begin{equation}
||\mathbf{a}||_1 = \sum_{i=1}^n\:|a_i|\quad\text{.}
\end{equation}

What is happening now is that the optimisation is penalising solutions with large sums of loading coefficients. This type of fit will avoid coefficient vectors which have a non-zero value for each model term. This is called the **Lasso** fit and is already implemented in the scikit-learn library.



```python
lassoreg = Lasso(alpha=0.007, max_iter=1e5)  # set up a Lasso regression
lassoreg.fit(A,y_data)                       # apply it to our model matrix A, and outcomes vector, y_data
```


```python
coeffs=np.copy(lassoreg.coef_)               # form our coeffs vector
coeffs[0]=lassoreg.intercept_
coeffs
```

Lasso regression promotes a *sparse* set of model coefficients.

We can see how our model fits our data, 


```python
y_fit = np.zeros(npts)
for n in range(npoly):
    y_fit = y_fit + coeffs[n]*x_data**n
plt.scatter(x_data, y_data, c='r', alpha=0.5)
plt.plot(x_data, y_fit, 'b', linewidth=2);
```

and also how it can be used to extrapolate.


```python
x_test=np.linspace(-0.2,1.2,npts)
y_fit = np.zeros(npts)
for n in range(npoly):
    y_fit = y_fit + coeffs[n]*x_test**n
plt.scatter(x_data, y_data, c='r', alpha=0.5)
plt.plot(x_test, y_fit, 'b', linewidth=2);
```

In the Lasso regression, we must choose a value of $\alpha$. What happens to the fit as we change $\alpha$?

## 8. A more realistic example
Imagine we measured the loss coefficient $Y_p$ of a blade at a range of incidences $i$. We would obtain a 'loss bucket', perhaps subject to some measurement noise: 


```python
npts=100
i=np.linspace(-10,10,npts)
Yp = 0.01 + i*i/500
```


```python
Yp = Yp * (1+ 0.1*np.random.normal(size=npts))
plt.scatter(i, Yp, c='r', alpha=0.5)
plt.xlabel('i')
plt.ylabel('Yp');
```

We have set $Y_p$ to be a quadratic function of $i$. But if we didn't know this beforehand, we may also be trying to fit $Y_p$ to some other variables that we had measured. We simulate this be adding two columns to our model matrix $\mathbf{A}$ that contain random numbers.


```python
A = 0.01*np.random.normal(size=[npts,5])   # First fill A with random numbers
A[:,0] = 1.                                # Then overwrite first 3 columns with 1, i and i*i
A[:,1] = i
A[:,2] = i*i
coeffs = np.linalg.pinv(A) @ Yp            # We use the pseudo-inverse to perform an l2 regression
coeffs
```

Our `coeffs` vector contains loadings of all the model terms, even the random noise.

It will still fit our data:


```python
fit = A @ coeffs                           # If we are applying our model to the independent variables in the data used to form the model, we can use this matrix multiply to obtain the fit
plt.scatter(i,Yp,c='r',alpha=0.5)
plt.plot(i,fit,'b',linewidth=2)
plt.xlabel('i')
plt.ylabel('Yp');
```

Now try with Lasso:


```python
lassoreg = Lasso(alpha=0.0003, max_iter=1e5)
lassoreg.fit(A,Yp)
coeffs=np.copy(lassoreg.coef_)
coeffs[0]=lassoreg.intercept_
coeffs
```

The Lasso regression promotes a sparse set of coefficients (try varying $\alpha$)


```python
fit = A @ coeffs
plt.scatter(i,Yp,c='r',alpha=0.5)
plt.plot(i,fit,'b',linewidth=2)
plt.xlabel('i')
plt.ylabel('Yp');
```


```python

```

3e684d7543e8020e63b2cae52267e5c094a07e35
ipynb
Jupyter Notebook
03 - Fitting Models to Data.ipynb
grahampullan/datascienceBB
87c1127d9deddd3e4e11f1a785b4f9e93a74fb19
MIT
03 - Fitting Models to Data.ipynb
grahampullan/datascienceBB
87c1127d9deddd3e4e11f1a785b4f9e93a74fb19
MIT
03 - Fitting Models to Data.ipynb
grahampullan/datascienceBB
87c1127d9deddd3e4e11f1a785b4f9e93a74fb19
MIT
Qwen/Qwen-72B
1. YES
2. YES