Regularization¶

Introductory Theory¶

Regularization is a technique for reducing the complexity of a model in supervised learning. The processing of performing regularization can be viewed with multiple interpretations:

Complexity Penalty

In the ERM framework, regularization can be viewed as the penalty or cost for adding extra complexity to a model. Usually, the complexity is a function of the weights learned for each feature. For now, we'll consider the class of linear models: $f(X) = W\cdot X+t$. Remember in ERM we seek a function $f^{opt}$ such that:

$f^{opt}= \underset{f \in \mathbb{F}} {\mathrm{argmin}} \frac{1}{n} \sum\limits_{i=1}^n \mathbb{L}(f(x_i),y_i)$

where $\mathbb{L}(f(x),y)$ is some arbitrary loss-function. The complexity of linear functions $f$ is set by the $L^p$-norm of the vector $W$, which is defined as $\|W\|_p =(|W_1|^p+|W_2|^p \dots |W_m|^p)^{\frac{1}{p}}$. The two norms we'll deal with the most in this class are $\|W\|_1$ and $\|W\|_2$. We want to add a term to the loss function that penalizes it as the norm of $W$ increases. We do this by adding a constraint to the above function:

$f^{opt}= \underset{f \in \mathbb{F}} {\mathrm{argmin}} \frac{1}{n} \sum\limits_{i=1}^n \mathbb{L}(f(x_i),y_i)$ subject to $R(W)\leq t$

where $R(W)$ can be any $L^p$-norm (mostly we'll choose $\|W\|_1$ or $\|W\|_2$). Using the method of Lagrange multipliers, we can set this up as:

$f^{opt}= \underset{f \in \mathbb{F}} {\mathrm{argmin}} \frac{1}{n} \sum\limits_{i=1}^n \mathbb{L}(f(x_i),y_i)+\lambda (R(W)-t)=\underset{f \in \mathbb{F}} {\mathrm{argmin}} \frac{1}{n} \sum\limits_{i=1}^n \mathbb{L}(f(x_i),y_i)+\lambda R(W)$

The parameter $\lambda$ is the Lagrange multiplier here, and in most cases we choose it using some form of cross-validation or model selection. We can think of it as a lever that controls the amount of penalty we enforce as $R(W)$ increases. The most commonly used forms of regularization are:

$f^{opt}= \underset{f \in \mathbb{F}} {\mathrm{argmin}} \frac{1}{n} \sum\limits_{i=1}^n \mathbb{L}(f(x_i),y_i)+\lambda \sum\limits_{j=1}^m |W_j|$

$f^{opt}= \underset{f \in \mathbb{F}} {\mathrm{argmin}} \frac{1}{n} \sum\limits_{i=1}^n \mathbb{L}(f(x_i),y_i)+\lambda \sum\limits_{j=1}^m W_j^2$

The first one above is commonly called "$L1$" regularization or the "lasso." The second is usually called "$L2$" or "ridge." Both can be thought of as penalties that force the optimal solution into a bounded region, where the bound is defined by some sort of $L^p$ circle with radius $t$.

Bayesian Formulation

Regularization also has a very nice Bayesian interpretation, which we will develop here. Regularization can be thought of as a Maxium a Posteriori estimation problem, which is an extension of Maximum Likelihood estimation. If you recall, we used the MLE method to find optimal parameters for Logistic Regression. This was defined as:

$\hat{\beta}_{MLE} =\underset{\beta} {\mathrm{argmin}} \:L(\beta|X,Y)=\underset{\beta} {\mathrm{argmin}} \:\prod\limits_{i=1}^nP(x_i,y_i|\beta)=\underset{\beta} {\mathrm{argmin}} \:\prod\limits_{i=1}^np_i^{y_i}(1-p_i)^{1-y_i}$

In MAP estimation, we assert a prior belief onto the value of $\beta$, so instead of finding $\beta$ that optimizes the likelihood, we find one that optimizes the posterior distribution, which is defined as:

$P(\beta|X,Y) = \frac{P(X,Y|\beta)*P(\beta)}{P(X,Y)}$

Usually with MAP, the denominator above doesn't factor into the optimization, so we only consider:

$P(\beta|X,Y) \propto P(X,Y|\beta)*P(\beta)=Likelihood*Prior$

For logistic regression, we have already defined the likelihood. Now we can define a prior on $\beta$. For reasons of mostly mathematical convenience, we can assume that $\beta$ comes from one of two priors:

$P(\beta_j) = N(0,\tau_j)=\frac{1}{\sqrt{2\pi \tau_j}}exp(\frac{-\beta_j^2}{2\tau_j}),\: \:j=1,2,\dots d$
or $P(\beta_j) = \frac{\lambda_j}{2}exp(-\lambda_j|\beta_j|), \: \:j=1,2,\dots d$

The first above is the normal distribution, and the second is the Laplace distribution. Again, there might be no reason to assume one prior is better (or more accuracte) than the other. We normally choose based on properties of the MAP, given a specific prior. Also, we defined the above priors as functions of individual parameters $\beta_j$. We usually assume that the prior on each $\beta_j$ is independent, so that $P(\beta)=\prod\limits_{i=1}^d P(\beta_j)$. We can now define a new optimization task, based on the MAP formulation:

$\hat{\beta}_{MAP} =\underset{\beta} {\mathrm{argmax}} \:L(\beta|X,Y)\:P(\beta)$

Using the above priors, we can define two different MAP estimators.

$\hat{\beta}_{MAP}=\underset{\beta} {\mathrm{argmax}} \prod\limits_{i=1}^n p_i^{y_i}(1-p_i)^{1-y_i} \: \prod\limits_{j=1}^d \frac{1}{\sqrt{2\pi \tau_j}} exp(\frac{-\beta_j^2}{2\tau_j})$
or $\hat{\beta}_{MAP}=\underset{\beta} {\mathrm{argmax}} \prod\limits_{i=1}^n p_i^{y_i}(1-p_i)^{1-y_i} \: \prod\limits_{j=1}^d \frac{\lambda_j}{2}exp(-\lambda_j|\beta_j|)$

We don't try to optimize these functions as is. Instead we take the log, and in doing so, we end up with equivalent formulation of $L1$ and $L2$ regularization defined above.

Examples¶

As a starting point, we can visualize how regularization constrains our search space for possible optimum values of $\beta$. We'll generate a high variance data set that represents a bournoulli process where $P(Y|X)$ is governed by the inverse logit. We'll then plot the following:

• Contour plot of the log-likelihood loss as a function of $\beta$
• Constraint boundaries of the following constraint $|\beta|_p
• The path of $\hat{\beta}_{MAP}$ as we vary regularization strength

The next two plots show all of the above with different sized datasets. The first is relatively small, with 50 examples, and the 2nd has 50k examples. We know the true value of the data generating distribution is $\beta=[1, 1]$. The dark green elipsoid shows where the log-loss is at a minimum. Its center is the optimum point of the data. With more data, we can see that this center gets closer to $[1, 1]$.

The circular constraint boundaries basically say that the coordinates of $\hat{\beta}_{MAP}$ have to be within that boundary. The MAP estimator will found the lowest value of the log-loss (i.e., darker color) within the circle.

What we can learn here is that with a very loose regularization weight (bigger circle), we are more likely to choose $\hat{\beta}_{MAP}$ that fits the data well, but is far from the truth (which we know because we created the data generating distribution). With the right regulariztion strength (such as the red circle in the above plot), we can see that our MAP estimate is actually closer to the truth. The problem is less severe with more data, because the dark green elipsoid shifts closer to the true optimum point.

Note that with $L1$ regularization, we can do the same analysis, though the constraint boundaries would look different. Instead of circles they would be straight lines.

The above examples were a bit theoretical, and are mainly presented for explaining the underlying optimization process (and to reinforce the idea of the bias-variance tradeoff going on).

Below we'll look at an actual dataset. First we'll read it in and rescale it. Rescaling is very important because the norm of $\beta$ is going to be a function of the scale of each feature. Since we apply the same regularization factor to each feature, rescaling ensures that it has the same effect per feature.

Next, we'll look at how a feature's weight evolves as we change the regularization weight.

This plot shows how the weights progress with $L2$ regularization. We can see a generally smooth progression that converges to a point as the regularization strength decreases.

Next we'll do the same for $L1$.

It is not so obvious from the above plot that $L1$ holds the feature weight at 0 for some time, but we can see this if we directly compare the paths of individual features.

Next we want to choose our regularization strength. Remember that from a Bayesian point of view, the regularization strength respresents our prior belief on the variance of the feature weights. If we assume a lower prior variance, we are effectively not allowing the weights to deviate too far from zero. Whether we take this point of view or not, the optimal strength is largely an empirical question. It is generally better to test rather than to assume what the optimal strength should be.

The various implementations of learning algorithms differ in how they treat the regularization parameter. Many implementations formulate the loss function as such:

$\underset{f \in \mathbb{F}} {\mathrm{argmin}} \frac{1}{n} \sum\limits_{i=1}^n \mathbb{L}(f(x_i),y_i)+\lambda R(W)$

In the above, the parameter $\lambda$ controls how much of the loss function is determined by the regularization function. Other implementations view it like this:

$\underset{f \in \mathbb{F}} {\mathrm{argmin}} \frac{1}{n} \sum\limits_{i=1}^n C*\mathbb{L}(f(x_i),y_i)+R(W)$

In the latter formulation the parameter controls how much the training error influences the total loss. It doesn't really matter which of the above is used (one is just a scalar multuple of the other), but it is helpful to know which system your implementation uses. I.e., it is useful to know if increasing the regularization parameter increases or decreases the amount of regularization being done. In SKLearn, the latter formulation is used.

Here we take the above example and use cross validation to find an optimal regularization strength.

Now lets plot the cross-validated error as a function of C, with a separate series for the different penalization types

This is interesting. We can see that the AUC for $L1$ is $0.5$ where the regularization is stronger, but for $L2$ the AUC is pretty good in the same region. This is because all of the feature weights in $L1$ are zero for that range. In $L2$ they'll likely be small, but they won't be zero. This means they'll still be able to rank the instances.