2019-01-18

# Linear Regression - behind the scenes

## Which problem do we want to solve?

We are given some data points ${(x_1,y_1), (x_2,y_2), (x_3,y_3), (x_n,y_n)}$ which we call our training set. Then we get a new data point $x_{new}$ for which we want to predict the corresponding $y_{new}$ value in a way that it follows the same pattern as in our training set. Thus our task is to find a model which can capture/learn these underlying patterns (aka relationships) in our training set and apply them to new points to make accurate predictions.

Some models types can capture/learn only a specific type of pattern (e.g. Linear Regression), others are more powerful and can learn any type of pattern (e.g. Neural Networks). Unfortunately the more powerful models will not only learn useful patterns but also the noise in our data. Usually it’s best to pick the simplest model type which is yet able to learn the most important patterns in your data because the simpler the model is the less chance there is that it learns noise.

## Implicit model assumption for Linear Regression

Linear Regression is a simple model type which assumes that the relationship between $x$ and $y$ is linear. Because of the added bias term the relationship it is strictly speaking affine and not linear. But you can apply a neat trick to make it linear again.

For scalar valued $x$:
$y = w_1*x_1 + bias$

For vector valued $x$:
$y = w_1*x_1 + w_2*x_2 + w_3*x_3 + bias = w^{T}x + bias$

Trick: affine => linear
$y = w_1*x_1 + w_2*x_2 + w_3*x_3 + w_4*(+1) = w^{T}x$
By adding a new entry $+1$ to every data point $x$ we can absorb the bias term ($w_4*1 = bias$) and make the relationship between $x$ and $y$ linear.

## What if the assumption doesn’t hold?

Here are three examples of how the relationship between x and y might look like. The relationship can either be exactly linear, roughly linear or not linear at all.

We can fit a linear regression model in any case but the predictions won’t make much sense in the case when the data doesn’t follow a linear relationship.

## How can we find the best model paramter w?

1. You define a function which tells you how good your current estimate of the best paramters $w$ is. This function is usually called error function or loss function.
2. You minimize the error function by changing your paramter values $w$. When people say they fit a model they refer to this step of minimizng an error function. For some model types and error functions there exist a closed form solution for $w$ which minimizes the error function for others you need to apply other (often iterative) minimization methods, e.g. for logistic regression.

### Error meassure 1: least square error

The vectors and matrices have the following dimensions.
$w$: Dx1
$x_i$: Dx1
$y_i$: 1x1
$y$: Nx1
$X$: NxD

In order to minizme the error function we take the gradient with respect to $w$ set it to zero and solve for $w$.

$$There exists no closed form solution for the least absolute error because of the absolute value in the equation. That’s why an iterative approach like Gradient Descent is required which won’t be introduced at this point. For the first plot (above) the least square error was used to fit the model, for the plot below the least absolute error was used. The results are pretty much the same. ## Basis functions - how to deal with nonlinear data In case the relationship between $x$ and $y$ is not linear you have two options. You can either switch to a completely new type of model which is able to capture nonlinear relationships or you can transform your data into a new space where the relationship between your transformed data x and y is linear again. The function which transforms your data is generally called a basis function. Example 1: Let’s assume we have a quadratic relationship between $x$ and $y$ in the form of $y = 1.5x^2$. If we transform our data to a new space (again just a 1D space) with the new x axis $x_{transformed} = x^2$ then the relationship between $y$ and $x_{transformed}$ will be linear again $y = 1.5x_{transformed}$ as you can see in the plot below. Example 2: Usually you don’t know to which space you should transform your data in order to make the raltionship linear again. To overcome this issue you can transform your data to a very high dimensional space and let your model decide which transformations to use. For example assume the true relationship between $y$ and $x$ has the form of $y = 1.5x^3$ but you don’t know it in advance. Now you transform your data to a high dimensional space with the new axes $x_{transformed 1} = x^2$, $x_{transformed 2} = x^3$, $x_{transformed 3} = x^4$, $x_{transformed 4} = x^5$, $x_{transformed 5} = x^6$. This works well in theory as your model can always decide not to use a transformation by setting the corresponding model parameter to $0$. Example 3: Using this technique you can fit your model to data which only approximately follows a polynomial relationship like for example data which was generated in the follwing way: $y = sin(x)$. If we use the same transformation technique up to an exponent of 5 we get the following model when we fit our model with a least square error. You can see in the plot below that our model captures the data pretty well and it allows us to make good predictions. Example 4: So far it seems best to transform the data to an as high dimensional space as possible and let your model figure out which transformations to choose. But there is a drawback. Assume that your data was generated in the following way $y = x^2 + noise$. As before let’s transform our data up to an exponent of 40 and see what happens. In red you see the fitted model which uses the transformed high dimensional data and in green our desired model. This phenomen is called overfitting. Our model is to complex and uses the extra power (high dimensional features e.g. $x^{40}$) to follow not only the major underlying quadratic pattern but also the noise. In conclusion, you can say that our model is not able to figure our itself which transformations to use. When we provide higher dimensional features to our model it will not ignore the ones which are not necessary to describe the major underlying pattern but rather use them to also capture noise in our data. We somehow need to find the sweet spot of complexity where our model is powerful enugh to capture the important patterns but not powerful enough to also capture the noise. There is one technique called Regularization which can help us with that. ## Regularization Regularization describes methods which prevent our model from overfitting by reducing the complexity of your model. This is a quite a general definition but provides you with the right intuition for future more complex models where regularization can occure in many different ways and not just by tampering the error function which we will do in a second for linear regression. ### L2 Regularization for Linear Regression$$

L2 regularization adds a new term to the error function which prevents the parameters w to become too large. It uses the squared euclidean distance/squared L2 norm to do so. Let’s again minimize the error function with respect to $w$.

### L1 Regularization for Linear Regression

$$(Xw-y)^T(Xw-y) + \lambda\lVert w\rVert_1$$

L1 regularization also adds a new term to the error function which pretty much does the same thing as the l2 regularization but uses the L1 norm (not squared) instead.

Let’s look at an example similiar to the previous one. Our data was generated in the following way: $y = 8x^2 + 20sin(4x) + noise$. I played around with the coefficients to make the plots look nice. Here are the results.

Firstly, you can see that a model (red line) with regularization becomes smoother compared to a model without regularization. Secondly, you can see that L1 regularization leads to an even “smoother” model than L2 regularization. Why is this so? The L2 norm doesn’t penalize all errors in the same way. For example when your model is far off $w^Tx >> y$ you will get a huge error because of the quadratic term in the L2 norm. However when you are only slighly off $w^Tx ~ y$ then your error will be almost zero again because of the quadratic term. This is different in the L1 norm where you use the absolute value of the distance. That’s why L1 regularization can push parameters w down to zero whereas L2 regularization will make them small but not zero.

The following picture provides another intuitive way to grasp what regularization for linear regression does. Essentially regularization offers a way to shrink the complexity of our model without changing its structure.

## Behind the scenes: probabilistic view

Now we come to the fun part where we will explore the connection between a regularized non probabilistic model and a probabibilistic model with certain assumptions.

Our probabilistic model: $y_i = w^Tx_i + N(\epsilon|0,\sigma^2)$

As the noise is the only probabilistic term it follows directly: $p(y_i|x_i,w,\sigma^2) = N(\epsilon|w^Tx_i,\sigma^2)$

It’s a bit confusing at first to read those formulas. In our case we interpret $p(y_i|x_i,w,\sigma^2)$ not as a function of $y_i$ but as a function over $w$ and $\sigma^2$. This interpretation of the formula is generally called a likelihood function. The integral over $w$ and $\sigma^2$ doesn’t sum up to 1. That’s why it is not a probability distribution. If we would assume $w$ and $\sigma^2$ are given/fixed and we look at it as a function of $y_i$ then it is a proper probability distribution and the integral over $y_i$ does sum of to 1.

## How do we get w and sigma?

We are now interested in finding the unknown parameters $w$ and $\sigma^2$. In order to do so, we have three options.

### 1. Maximum Likelihood Estimation (MLE)

 We assume that our data is identical and independently distributed (aka i.i.d) meaning every data point was generated using the same distribution and all data points are independent of each other. With this assumption we can compute the likelihood function $$p(y,x w,\sigma^2) = \prod_{i=1}^{N} p(y_i,x_i ,w,\sigma^2)$. Now we can simply maximize this function in order to find the parameter$w$and$\sigma$$ which give us the most likely description for our observed data points.

Our task: $w*, \sigma^{2*} = \mathop{\arg\,\max}\limits_{w, \sigma^2} \prod_{i=1}^{N} p(y_i,x_i|,w,\sigma^2) = \mathop{\arg\,\max}\limits_{w, \sigma^2} L(w,\sigma^2)$

You might think that we will run into issues when we try to optimize over two variables simultaneously. Generally that can be an issue but not in this case as $w^*$ doesn’t depend on $\sigma^2$.

We made use of a simple trick. We took the log of the whole expression knowing that it doesn’t change the position of the maximum. In order to find the maximum, we take the gradient with respect to $w$, set it to zero and solve for $w$.

You can see that $w_{MLE}$ is independent of $\sigma^2$ and also that the result is the same as for the non probabilistic least square case with no regularization.

Let’s go on to derive the solution for $\sigma_{MLE}^2$. Same procedure, we take the gradient with respect to $\sigma^2$ , set it to 0 and solve for $\sigma^2$.

This is a satisfying result because it is just finding the sample average of the squared deviations between what $w_{MLE}$ predicts and what the training data actually are. It feels exactly like what happens when you compute the maximum likelihood estimate of the variance of a univariate Gaussian distribution.

Let’s quickly recap. The solution for $w$ with the probabilistic approach (+ our assumptions) and taking the Maximum Likelihood Estimate leads exactly to the same solution as for the non probabilistic model. Which means that by using the least squared error function we implictely assume that any noise we have has a mean of zero and is symmetric around the mean. We can’t assume that the noise strictly follows a normal distribution as we would get the same result when the noise would be e.g. uniformly distributed with a mean of 0. But at least we got some better intuition about what’s going on behind the scenes. Let’s go a step further and see if we can also find a corresponding probabilisitc assumption for the regularization term.

NOT QUITE SURE YET

### 2. Maximum A Posteriori Estimation (MAP)

Bayes formula: $\text{A posteriori }= \frac{\text{likelihood * prior}}{\text{evidence}}$

Translated to our model: $p(w, \sigma^2|x,y) = \frac{p(x,y|w, \sigma^2) * p(w,\sigma^2)}{p(x,y)}$

The evidence term doesn’t depend on $w$ or $\sigma^2$ and thus doesn’t affect the maximum of the posterior distribution $p(w, \sigma^2|x,y)$ . So in order to find the maximum of the posterior distribution we only need the likelihood function $p(x,y|w, \sigma^2)$ and the prior distribution $p(w,\sigma^2)$. We already know how to get the likelihood function from the previous step but how do we get the prior?

Eventually we need to make some hard assumptions about $w$ and $\sigma^2$. Let’s discuss our options. We could assume that both $w$ and $\sigma^2$ are indepedent random variables with follow a certain distribution which we can freely choose. For example $p(w, \sigma^2) = p(w) * p(\sigma^2) = N(w|a,b)*N(\sigma^2|c,d)$. Not much progress made, now we need to come up with good values for $a,b,c,d$. We could again assume they follow a certain assumptions to take uncertainty into account but you see that we are trapped in a loop. At some point we have a make a hard assumption and set a parameter to a fixed value. It’s good practice too pick a reasonable fixed value when you can’t make good asusmptions about the distribution of a parameter anymore. In our case we will assume $\sigma^2$ is not a random variable but a fixed scalar and $w$ is gaussian distributed with a variance of $\frac{1}{\lambda}$ thus $p(w) = N(\sigma^2|0,\frac{1}{\lambda})$. Our prior simplifies to $p(w, \sigma^2) = N(\sigma^2|0,\frac{1}{\lambda})$.

No we can ccalculate the maximum a poserior estimation for our only left random variable $w$.

Now we set the gradient with respect to $w$ to zero and solve for $w$.

Here we get another insight. If you go back to the solution for the optimal $w$ when we used a non probabilistic model with L2 regularization, you see that the result is exactly the same as for the probabilistic model if we assume $\sigma^2 =1$. Apparently using a gaussian prior corresponds to an added L2 reguarization term to our error function. If you play around a bit with different prior distributions you can find out that a laplacian prior corresponds to L1 regularization.

When I first learnt about these correspondences I was pretty amazed!

### 3. Fully Bayesian Analysis

So far we only did a so called point estimates for $w$. This can work out nicely but we are losing information about the uncertainty of our estimate for $w$w. If we want to keep the information about the uncertainty we need to compute or estimate the full posterior distribution. I won’t go into more detail here. We will discuss more about this in the next post “Probabilistic Inference - behind the scenes”.