### Preparation

* [`numpy`](http://www.numpy.org) is a module for scientific computing (i.e. *vector* and *matrix* processing) with Python.
* [`matplotlib`](https://matplotlib.org) is a module for plotting.


```python
import numpy as np
import matplotlib.pyplot as plt
import seaborn
```


```python
from operator import mul
from IPython.core.interactiveshell import InteractiveShell

InteractiveShell.ast_node_interactivity = "all"
plt.style.use('fivethirtyeight')
```

## An example

Let us look at a hypothetical appartment rent data. We will take a look only at the price per square meter ($m^2$): `price = w0 * 1  + w1 * square_meters`. We call `price` a response variable; in literature usually denoted as $y$. We call `square_meters` a predictor variable. We call `w0` and `w1` weights - it is basically *intercept* and *slope* of a line. `w0` has a special meaning and is called a `bias`. In our example we can interpret it as base fixed price you would pay for 0 sqaure meters of living area.

Here is a graph that visualizes a price with 200 CU (currency units) base price and 10 CU increase per square meter: `price = 200 + 10 * square_meters`. We also add some *observations*, i.e., some data that we read in a finctional announce or whereever. These data points have some *noise* - factors that we cannot account for. For example some renters may have good or bad negotiation skills that will make the price vary around the line.  


```python
square_meters = np.linspace(0, 150)
price = 200 + 10 * square_meters

# a random number between 0 and 1 will be multiplied by 150 to
# kind of represent a spectrum of square meters.
number_of_points = 20
random_square_meters_data = np.random.sample(size = number_of_points) * 150
price_with_noise = 200 + 10 * random_square_meters_data + np.random.normal(size = number_of_points, scale = 50)

# plot data
_ = plt.plot(square_meters, price, color = 'grey')
_ = plt.plot(random_square_meters_data, price_with_noise, 'o', color = 'blue')
_ = plt.xlabel('Square meters')
_ = plt.ylabel('Price in CU')
plt.ylim(0, 2000)
plt.show()
```

The data here has one dimension: the square meters. But we rarely have only one-dimensional data (also see [problems with many dimensions](https://en.wikipedia.org/wiki/Curse_of_dimensionality)). For example, square meters are actually two variables: length and width. We can also add height as a third dimension that may be an additional influence to the price. So is the city, weather and anything else. We usually abstract from the actual names and simly call dimensions by their numbered variable: dimension $i$ is simply encoded in variable $x_i$. For example $x_1$ is `length`, $x_2$ is `width`, $x_3$ is `height` and so on. *Bias* has usually a special notation: $x_0$. We almost always need to take care for *bias* variable and we will do it down bellow.

We now have a conceptual understanding of dimensions and can jump into coding stuff.

The goal of the following is to explain how basic regression works. 

Imagine we have gathered 20 samples of 'data' where for each price (response) we have `length`, `width`, `height`, and `distance` from downtown (predictors).


The weights influence the price as follows: we have `basic_price` ($w_0$), `length` ($x_1$), `width` ($x_2$), `height` ($x_3$), `distance` ($x_4$), and the last one `square_meters` ($x_5$). Thus the price results from the following: 
$$price = w_0 + w_1 \cdot x_1 + w_2 \cdot x_2 + w_3 \cdot x_3 + w_4 \cdot x_4 + w_5 \cdot x_5$$


```python
weights = (np.random.uniform(size = 6) + 1) * [200, .5, .7, 1.5, -.05, 3.5]
weights
```


```python
n_samples = 20
np.set_printoptions(2, suppress=True)
tmp = np.column_stack((1.5 + np.random.sample(size = n_samples) * 10, # length
                       1.5 + np.random.sample(size = n_samples) * 10, # width
                       np.round(2   + np.random.sample(size = n_samples) * 2, decimals=1), # height
                       np.random.sample(size = n_samples) * 3000)) # distance from center in meters
samples_data_X = np.column_stack((tmp, tmp[:,0]*tmp[:,1])) # square meters
del(tmp)
samples_data_X
```


```python
samples_data_price_Y = np.dot(np.column_stack((np.ones(samples_data_X.shape[0]), samples_data_X)), weights)
# we add some noise to the response that adds unpredictable variance to the price.
samples_data_price_Y = samples_data_price_Y + np.random.normal(scale = 2.5, size = np.prod(samples_data_price_Y.shape))
samples_data_price_Y
```


Let's say, we want to find *weights* of a *linear model* to predict the price of an offer (the six $w_i$ from the equation above). Meaning we got data from a makler that send us data for `length=10`, `width=10`, `height=3`. Furthermore we looked at some map data and got the `distance=3044` meters from downtown. The makler didn't tell us the price yet - but lets see if we can estimate the price by ourselfes. In such a way we will go into negotiations and can back up our claims with data.

Having gathered the data, we compute the following:

$$\mathbf{\hat{w}} = (\mathbf{X}^T \mathbf{X})^{-1}\mathbf{X}^T \mathbf{y} $$


```python
bold_X = np.column_stack((np.ones(samples_data_X.shape[0]), samples_data_X))
estimated_weights = np.dot(np.dot(np.linalg.inv(np.dot(bold_X.T, bold_X)), bold_X.T), samples_data_price_Y)
estimated_weights
weights
```

Looks like the estimated weight from noisy data are not so far off. Of course we wouldn't know the true weights in reality. And so, we would have to bargain with a makler for a price around:


```python
np.round(np.dot(estimated_weights, [1, 10, 10, 3, 3044, 10*10]), 2) 
```

### Here be dragons

**What is this matrix vodoo above? And why does it looks like it looks like?**

We will derive the above equation based on amazingly clear [tutorial of Mark L. J. Orr](https://www.cc.gatech.edu/~isbell/tutorials/rbf-intro.pdf). And we will do it in the least mathematical way possible - mostly with **code**.

## Appendix A.1

A *vector* is simply an array of $n$ numbers $\mathbf{x}^T = [x_1, x_2, \ldots, x_n]$. The number $n$ represents the number of *dimensions*. The T-operator $\cdot^T$ simply tells us to treat the vector as a *column* vector - more to that in a second.


```python
x = np.arange(5)
x
```

A matrix is simply a 'collection' of vectors with the same dimension. Most of the time the matrix will have $m$ number of rows and $n$ number of columns. 

**Example**: We want to have two (2) arrays. Each of array contains five (5) numbers, i.e., is a 5-dimensional vector. The resulting matrix has $m=2$ rows and $n=5$ columns. 


```python
v1 = np.arange(1, 6)
v1
v2 = np.arange(1, 6) * 2
v2
simple_matrix = np.stack([v1, v2])
simple_matrix
```

For computational purpuses we seldomly use $x$ as single vector. Instead we treat it as *row* entry in a matrix, i.e. as a *row vector*. Hence, in matrix notation, a vector has one (1) row and $n$ columns. (This is usually denoted as a *tuple* `(1, n)`)


```python
v1 = np.reshape(np.arange(1,6), (1,5))
v1
```

Remember our T-operator? Now, if we use $\cdot^T$ it will *[T]ranspose* the matrix, i.e. rows will be collums and collumns will be rows.


```python
v1.T
```

The same applies to matrices:


```python
simple_matrix = np.stack([v1, v2])
simple_matrix
simple_matrix.T
```

Which bring us directly to matrix and vector operations. We denote a $k$-th vector with bold $\mathbf{x}_k$ (e.g. vector 1 is $\mathbf{x}_1$). We denote a matrix with a bold capital letter, for example $M$. As discussed, we can access $j$-th vector (i.e. $\mathbf{x}_j$) and $i$-th value (value in dimension $i$) - or in short: the value of $M_{ji}.

We have the following operations: (TODO implement in actual code; compare timings with numpy)

* vector addition 
* vector and scalar multiplication
* vector matrix product
* matrix product 

About matrix product: you can only multiply matrices when the 'inner' shape matches: `(shape_1, inner_shape) x (inner_shape, shape_2)`. The resulting shape of the matrix will be `(shape_1, shape_2)`. Here the T-operator comes in play because it reverses the shape. This allow us to multiply matrix and its transpose form. The resulting shape will have the same number of rows and columns and is termed quadratic matrix. One of the important properties of the quadratic matrix is that it is inversible. The inverse of a matrix describes the same effect for simple number. If you have a number $a$ - the inverse is $\frac{1}{a}$ or $a^{-1}$ which means that $a \cdot a^{-1} = 1$.
The same is true for matrices where $\mathbf{M} \cdot \mathbf{M}^{-1} = \mathbf{I}$. $\mathbf{I}$ is the unity matrix - a matrix where the only non-zero entries are in the diagonal with value 1, or in other terms $I_{ii} = 1$. The identity matrix' purpose is to have an operation which yields the same 'object' when multiplied. 


```python
vector = np.arange(5)
vector
vector + 5 # 5 is 'broadcasted' into the array - i.e. + 5 is repeated for every entry
vector + vector # adding vectors element wise
vector * 3 # multiply every entry with 3
vector * vector # element wise multiplication | THAT IS NOT A VECTOR PRODUCT | it is called Hadaman-Product
np.dot(vector, vector) # that is a vector product; it results in a scalar
mat_1 = np.reshape(np.random.randint(1,100, size = 2*3), (2,3))
mat_2 = np.reshape(np.random.randint(1,100, size = 2*3), (2,3))
mat_1
mat_2
np.dot(mat_1, mat_2.T)
np.dot(mat_1, mat_2.T).shape
np.dot(mat_1.T, mat_2)
np.dot(mat_1.T, mat_2).shape
```

#### Why do we use matrix operations ?

In linear models we usually compute. 
$$y = w_0 \cdot 1 + w_1 \cdot x_1 + w_2 \cdot x_2 + w_3 \cdot x_3 + w_4 \cdot x_4 + w_5 \cdot x_5$$
$$y = w_0 + \sum\limits_{i=1}^n w_i \cdot x_i$$

For efficiency reasons we transform the vector $\mathbf{x} = (x_1, x_2, x_3, x_4, x_5)$ into $\mathbf{x} = (1, x_1, x_2, x_3, x_4, x_5)$


```python
example_weights = np.array([20,4,6])
example_weights
example_x = np.array([1, 4, - 5])
example_x
```

The equation above can basically be implemented like this (don't! it is highly ineffcient!):


```python
def get_y(w, x):
    res = 0
    for i in range(len(x)):
        res = res + w[i]*x[i]
    return res
```


```python
get_y(example_weights, example_x)
```

A *vectorized* version is actually 'better':


```python
sum(example_weights * example_x)
```

For matrix operation we need to treat the vectors as tuples (1 row of 3 dimensions, `(1,3)`). The ultimate operation is the *dot-product* and since we need matching of inner dimensions, we need to prepare the vectors to return one number, i.e. the outer shape is 1. Thus the dot-product between these two matrix-vectors can be written as $\mathbf{w}\cdot\mathbf{x}^T$ for *row-vectors* (and actually in mathematical terms for *column vectors*: $\mathbf{w}^T\mathbf{x}$, but never mind).


```python
example_weights2 = np.reshape(example_weights, (1, 3))
example_x2 = np.reshape(example_x, (1, 3))
np.dot(example_weights2, example_x2.T)
```

The timings speak for themselfes:


```python
%timeit get_y(example_weights, example_x)
%timeit sum(example_weights * example_x)
%timeit np.dot(example_weights2, example_x2.T)
```

And that, just for one single small sample. We usually are dealing with more. How about 10000 samples with 1000 dimensions?


```python
example_weights2 = np.reshape(np.random.sample(size = 1*1000), (1, 1000))
example_x2 = np.reshape(np.random.sample(size = 10000*1000), (10000, 1000))
```


```python
%timeit y1 = [get_y(example_weights2[0], x) for x in example_x2]
%timeit y2 = [sum(example_weights2[0] * x) for x in example_x2]
%timeit y3 = np.dot(example_weights2, example_x2.T)
```

**This is the reason why matrix operations are so important!** (not only because of math)!

## Appendix A.4 The optimal weight vector

Now we are getting into the bottom of things: So given data in $\mathbf{X}$ (m rows of n dimensional vectors) we want to determine weights (n (+1) dimensional vector) of a linear model. As [exemplarly shown](https://www.geogebra.org/m/xC6zq7Zv), the goal is to find such a weight vector that minimizes the '*reconstruction error*' of the estimation represented by a *sum of squared errors* (SSE):
$$E = \frac{1}{2}\sum\limits_{\mathbf{x}_i \in \mathbf{X}, \hat{y}_i} (\hat{y}_i - f(\mathbf{x}_i))^2$$
with $f(\mathbf{x}_i)$ beeing our estimator:
$$f(\mathbf{x}) = \sum\limits_{k=0}^n w_k \cdot x_k$$
(we add $\frac{1}{2}$ to make our equations nicer ;-)

Now, to find a minimum of a function (and we want to find parameters such that the error $E$ is minimal) math tells us to
* differentiate the function with respect to the variables (in our case $w_i$ from $\mathbf{w}$)
* eqaute the result with zero
* solve the equation for the variables

Let's to just that:
$$\frac{\delta E}{\delta w_j} = \frac{\delta \frac{1}{2} \sum\limits_{\mathbf{x}_i \in \mathbf{X}, \hat{y}_i} (\hat{y}_i - f(\mathbf{x}_i))^2}{\delta w_j}$$
which is equally valid to
$$\frac{\delta E}{\delta w_j} = \frac{1}{2} \sum\limits_{\mathbf{x}_i \in \mathbf{X}, \hat{y}_i} \frac{\delta (\hat{y}_i - f(\mathbf{x}_i))^2}{\delta w_j}$$

Let us focus on the term $\frac{\delta (\hat{y}_i - f(\mathbf{x}_i))^2}{\delta w_j}$. A simple application of [chain rule](https://en.wikipedia.org/wiki/Chain_rule) tells us that $(h(g(x)))' =  h'(g(x))\cdot g'(x)$ with $h(x)~~\hat{=}~~x^2$ and $g(x)~~\hat{=}~~\hat{y}_i - f(\mathbf{x}_i)$. The derivate of $h(x)$ is $h'(x) = 2\cdot x$. And the derivate of $g(x)$ is $g'(x) = - x_{ij}$.

Substituting the variables yields:
$$\frac{\delta (\hat{y}_i - f(\mathbf{x}_i))^2}{\delta w_j} = 2 \cdot (\hat{y}_i - f(\mathbf{x}_i)) \cdot {-  x_{ij}}$$

Inserting into the error $E$ we get
$$\frac{\delta E}{\delta w_j} =  \frac{1}{2} \sum\limits_{\mathbf{x}_i \in \mathbf{X}, \hat{y}_i} 2 \cdot (\hat{y}_i - f(\mathbf{x}_i)) \cdot {-  x_{ij}}$$
and equating to zero we get:
\begin{align}
0 & = \sum\limits_{\mathbf{x}_i \in \mathbf{X}, \hat{y}_i} (\hat{y}_i - f(\mathbf{x}_i)) \cdot {-  x_{ij}} \\
0 & = \sum\limits_{\mathbf{x}_i \in \mathbf{X}, \hat{y}_i} - \hat{y}_i \cdot x_{ij} + f(\mathbf{x}_i) \cdot x_{ij} \\
0 & = - \sum\limits_{\mathbf{x}_i \in \mathbf{X}, \hat{y}_i} \hat{y}_i \cdot x_{ij} + \sum\limits_{\mathbf{x}_i \in \mathbf{X}, \hat{y}_i} f(\mathbf{x}_i) \cdot x_{ij} \\
\sum\limits_{\mathbf{x}_i \in \mathbf{X}, \hat{y}_i} \hat{y}_i \cdot x_{ij} & = \sum\limits_{\mathbf{x}_i \in \mathbf{X}, \hat{y}_i} f(\mathbf{x}_i) \cdot x_{ij} \\
\end{align}

Now, as we saw earlier, we can rewrite the sums over samples in *matrix form*:

$$\mathbf{X}^T\cdot \hat{\mathbf{y}} = \mathbf{X}^T\cdot \mathbf{f}$$

Our $\mathbf{f}$ is a vector of linear model estimation of $\mathbf{x}_j$: 
$$f_j = f(\mathbf{x}_j) = \sum\limits_{i=0}^n w_i \cdot x_{i}$$
(with $x_0 = 1$) and hence can be rewritten to $\mathbf{x}_j^T\hat{\mathbf{w}}$ as we saw earlier.

$$
\mathbf{f} = \begin{bmatrix}f_1\\\ldots\\f_m\end{bmatrix} = \begin{bmatrix}\mathbf{x}_1^T\hat{\mathbf{w}}\\\ldots\\\mathbf{x}_m^T\hat{\mathbf{w}}\end{bmatrix} = \mathbf{X}\mathbf{\hat{w}}
$$

Then, from above, we have:
\begin{align}
\mathbf{X}^T \hat{\mathbf{y}} & = \mathbf{X}^T \cdot \mathbf{f} \\
\mathbf{X}^T \hat{\mathbf{y}} & = \mathbf{X}^T \mathbf{X} \cdot \mathbf{\hat{w}}
\end{align}

And solving for $\hat{w}$ (i.e. 'dividing' by $\mathbf{X}^T \mathbf{X}$ or better yet: multiplying by the inverse $(\mathbf{X}^T \mathbf{X})^{-1}$) finally yields:

$$ \mathbf{\hat{w}} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \hat{\mathbf{y}}$$

This is our solution from above. Simply coded as `np.dot(np.dot(np.linalg.inv(np.dot(X.T, X)), X.T), Y)`.

## Runtime information

Numpy, BLAS, and LAPACK.


```python
np.__file__
np.__version__
np.show_config()
```


```bash
%%bash
ldd ~/anaconda3/lib/python3.6/site-packages/numpy/core/multiarray.cpython-36m-x86_64-linux-gnu.so
```

d4f58229c32d7130988043a01f1fea009f77f569
ipynb
Jupyter Notebook
linear-models.ipynb
burkovae/linear-models
e8db6856d8e822ec31ec8eead6b509a7fc38e6fc
MIT
linear-models.ipynb
burkovae/linear-models
e8db6856d8e822ec31ec8eead6b509a7fc38e6fc
MIT
linear-models.ipynb
burkovae/linear-models
e8db6856d8e822ec31ec8eead6b509a7fc38e6fc
MIT
Qwen/Qwen-72B
1. YES
2. YES