```python
import numpy as np
import matplotlib.pyplot as plt
import math
from numpy.linalg import inv
from numpy import matmul as mm
```

- Dataset $\mathcal{D}=\{x_i,y_i\}^N_{i=1}$ of N pairs of inputs $x_i$ and targets $y_i$. This data can be measurements in an experiment.
- Goal: predict target $y_*$ associated to any arbitrary input $x_*$. This is known as a **regression** task in machine learning.

### Generate Dataset
$$y_i = (x_i+1)^3+\epsilon_i$$
where $\epsilon_i \sim \mathcal{N}(0,1)$


```python
x = np.linspace(-4,2,20)
y = (x+1)**3+10*np.random.normal(0,1,20)
```


```python
plt.plot(x,y,'b.')
```

### Model of the data
Assume that the dataset can be fit with a $M^{th}$ order polynomial (**polynomial regression**),

$$f_w(x) = w_0+w_1x+w_2x^2+w_3x^3+...+w_Mx^M=\sum^M_{j=1}w_j\phi_j(x)$$
The $w_j$ are the weights of the polynomial, the parameters of the model, and $\phi_j(x)$ is a basis function of our linear in the parameters model.

### Fitting model parameters via the least squares approach
- Measure the quality of the fit to the training data.
- For each train point, measure the squared error $e^2_i = (y_i-f(x_i))^2$.
- Find the parameters that minimize the sum of squared errors:

$$E(\mathbf{w})=\sum^N_{i=1}e^2_i=\|\mathbf{e}\|^2 = \mathbf{e}^\top\mathbf{e}=(\mathbf{y}-\mathbf{f})^\top(\mathbf{y}-\mathbf{f})$$
where $\mathbf{y} = [y_1,...y_N]^\top$ is a vector that stacks the N training targets, $\mathbf{f}=[f_\mathbf{W}(x_1),...,f_\mathbf{w}(x_N)]^\top$ stacks the prediction evaluated at the N training inputs.

Therefore, 

\begin{align}
\mathbf{y}&=\mathbf{f}+\mathbf{e}=\mathbf{\Phi w}+\mathbf{e} \\
\begin{pmatrix} y_1\\y_2\\\vdots\\y_N\end{pmatrix}&=\begin{pmatrix}1&x_1&x_1^2&...&x_1^M\\1&x_2&x_2^2&...&x_2^M\\\vdots&\vdots&\vdots&\cdots&\vdots\\1&x_N&x_N^2&...&x_N^M\end{pmatrix}\begin{pmatrix}w_0\\w_1\\\vdots\\w_M\end{pmatrix} + \mathbf{e}
\end{align}


The sum of squared errors is a convex function of $\mathbf{w}$: 

$$E(\mathbf{w})=(\mathbf{y}-\mathbf{\Phi w})^\top(\mathbf{y}-\mathbf{\Phi w})$$
To minimize the errors, find the weight vector $\mathbf{\hat{w}}$ that sets the gradient with respect to the weights to zero, 

$$\frac{\partial E(\mathbf{w})}{\partial \mathbf{w}}=-2\mathbf{\Phi}^\top(\mathbf{y}-\mathbf{\Phi w})=2\mathbf{\Phi^\top\Phi w}-2\mathbf{\Phi}^\top \mathbf{y}=0$$

The weight vector is 

$$\mathbf{\hat{w}}=(\mathbf{\Phi^\top\Phi})^{-1}\mathbf{\Phi^\top y}$$


```python
def polynomialFit(x,y,order=3):
    for i in range(order+1):
        if i == 0:
            Phi = x**i
        else:
            Phi = np.vstack((Phi,x**i))
    Phi = Phi.T
    if order == 0:
        Phi = Phi.reshape(-1,1)
    w = mm(mm(inv(mm(Phi.T,Phi)),Phi.T),y)
    f = mm(Phi,w)
    dif = y-f
    err = mm(dif.T,dif)
    return f,err,w,Phi
```


```python
f,err,w,Phi = polynomialFit(x,y)
```


```python
plt.plot(x,y,'b.')
plt.plot(x,f,'r-')
print(w) # ideal: 1,3,3,1
print(err)
```

### M-th order Polynomial


```python
errlist = []
plt.figure(figsize=(20,20))
for i in range(21):
    plt.subplot(7,3,i+1)
    f,err,w,Phi = polynomialFit(x,y,i)
    errlist.append(err)
    plt.plot(x,y,'b.')
    plt.plot(x,f,'r-')
    plt.title('Order '+str(i)+': '+str(err))
```


```python
plt.plot(np.arange(16),errlist[:16])
```

#### The fitting becomes very unstable after the order of 15. This may be due to the inverse instability. This can be resolved via LU decomposition, Cholesky decomposition or QR decomposition.

### LU decomposition


```python
import scipy
from scipy.linalg import lu_factor,lu_solve
```


```python
def polynomialFitLU(x,y,order=3):
    for i in range(order+1):
        if i == 0:
            Phi = x**i
        else:
            Phi = np.vstack((Phi,x**i))
    Phi = Phi.T
    if order == 0:
        Phi = Phi.reshape(-1,1)
    lu,piv = lu_factor(mm(Phi.T,Phi))
    tmp = lu_solve((lu,piv),Phi.T)
    w = mm(tmp,y)
    f = mm(Phi,w)
    dif = y-f
    err = mm(dif.T,dif)
    return f,err,w,Phi
```


```python
errlistLU = []
plt.figure(figsize=(20,20))
for i in range(21):
    plt.subplot(7,3,i+1)
    f,err,w,Phi = polynomialFitLU(x,y,i)
    errlistLU.append(err)
    plt.plot(x,y,'b.')
    plt.plot(x,f,'r-')
    plt.title('Order '+str(i)+': '+str(err))
```


```python
plt.plot(np.arange(21),errlistLU)
```

### Cholesky decomposition


```python
from scipy.linalg import cho_factor,cho_solve
```


```python
def polynomialFitChol(x,y,order=3):
    for i in range(order+1):
        if i == 0:
            Phi = x**i
        else:
            Phi = np.vstack((Phi,x**i))
    Phi = Phi.T
    if order == 0:
        Phi = Phi.reshape(-1,1)
    c,low = cho_factor(mm(Phi.T,Phi))
    tmp = cho_solve((c,low),Phi.T)
    w = mm(tmp,y)
    f = mm(Phi,w)
    dif = y-f
    err = mm(dif.T,dif)
    return f,err,w,Phi
```


```python
errlistChol = []
plt.figure(figsize=(20,20))
for i in range(21):
    plt.subplot(7,3,i+1)
    f,err,w,Phi = polynomialFitLU(x,y,i)
    errlistChol.append(err)
    plt.plot(x,y,'b.')
    plt.plot(x,f,'r-')
    plt.title('Order '+str(i)+': '+str(err))
```


```python
plt.plot(np.arange(21),errlistChol)
```

### Comparison between inverse, LU decomposition and Cholesky decomposition


```python
plt.plot(np.arange(21),errlist)
plt.plot(np.arange(21),errlistLU)
plt.plot(np.arange(21),errlistChol)
```


```python
plt.plot(np.arange(21),errlistLU)
plt.plot(np.arange(21),errlistChol)
```


```python

```

7ab295bd291dc9cf6f1fa2a3cc07e5cb0e8e3672
ipynb
Jupyter Notebook
Probabilistic Machine Learning/Introduction to Probabilistic Machine Learning/.ipynb_checkpoints/Linear in the parameters regression-checkpoint.ipynb
zcemycl/ProbabilisticPerspectiveMachineLearning
8291bc6cb935c5b5f9a88f7b436e6e42716c21ae
MIT
2019-11-20T10:20:29.000Z
2021-11-09T11:15:23.000Z
Probabilistic Machine Learning/Introduction to Probabilistic Machine Learning/.ipynb_checkpoints/Linear in the parameters regression-checkpoint.ipynb
zcemycl/ProbabilisticPerspectiveMachineLearning
8291bc6cb935c5b5f9a88f7b436e6e42716c21ae
MIT
Probabilistic Machine Learning/Introduction to Probabilistic Machine Learning/.ipynb_checkpoints/Linear in the parameters regression-checkpoint.ipynb
zcemycl/ProbabilisticPerspectiveMachineLearning
8291bc6cb935c5b5f9a88f7b436e6e42716c21ae
MIT
2020-05-27T03:56:38.000Z
2021-05-02T13:15:42.000Z
Qwen/Qwen-72B
1. YES
2. YES